[PATCH] UTM zone number is always zero in NITF creation

[rsoconnor]

I created a NITF file with GDALCreate with the ICORDS variable set to "N" for UTM north. The problem is that the zone number is always coming out zero when I read the file back in.

I called the following functions after I wrote to the dataset with RasterIO:
oSRS.SetWellKnownGeogCS ("WGS84");
oSRS.SetUTM (zoneNumber, TRUE);
oSRS.exportToWkt (&pszDstWKT);
GDALSetProjection( pDstDS, pszDstWKT );
GDALSetGeoTransform( pDstDS, dfTransform );

BTW, the zone number comes out fine if I choose the GTiff driver instead.

[Even Rouault]

I think this is a known restriction. See ​http://www.gdal.org/frmt_nitf.html

Extract: On export NITF files are always written as NITF 2.1 with one uncompressed image and no other auxilary layers. Georeferencing can only be written for images using a geographic coordinate system. Coordinates are implicitly treated as WGS84 even if they are actually in a different geographic coordinate system. Pseudo-color tables may be written for 8bit images.

ICORDS: Set to "G" to ensure that space will be reserved for geographic corner coordinates to be set later via SetGeoTransform()

[warmerdam]

In particular, there is no actual implementation of SetProjection() on the NITFDataset(). The ICORDS=N creation options sets up a UTM coordinate system but leaves the zone part at zero (implicitly a central meridian of -183). I tried this program:

#include "gdal_priv.h"
#include "ogr_spatialref.h"

int main()

{
    GDALAllRegister();

    OGRSpatialReference oSRS;
    char *pszDstWKT = NULL;

    oSRS.SetWellKnownGeogCS ("WGS84");
    oSRS.SetUTM (11, TRUE);
    oSRS.exportToWkt (&pszDstWKT);

    printf( "%s", pszDstWKT );

    char *apaszOptions[] = {"ICORDS=N",NULL};
    GDALDriver *poDriver = (GDALDriver *) GDALGetDriverByName( "NITF" );
    GDALDataset *poDS = poDriver->Create( "out.ntf", 5, 5, 1, GDT_Byte, 
                                          apaszOptions );
    double adfTransform[6] = {30000, 10, 0, 1200000, 0, -10 };

    GDALSetProjection( poDS, pszDstWKT );
    GDALSetGeoTransform( poDS, adfTransform );

    GDALClose( poDS );
}

The resulting file has a gdalinfo report of:

Driver: NITF/National Imagery Transmission Format
Files: out.ntf
Size is 5, 5
Coordinate System is:
PROJCS["UTM Zone 0, Northern Hemisphere",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AXIS["Lat",NORTH],
        AXIS["Long",EAST],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-183],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (30000.000000000000000,1200000.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  NITF_FHDR=NITF02.10
...

As Even points out this is more-or-less a known limitation of the driver (though I had to rediscover it for myself by experimenting). I am going to defer this bug report as an enhancement request. It should be practical (for instance) to have SetProjection() reset the zone info for the UTM case.

PS. You can create UTM NITF files using the CreateCopy() method.

[Even Rouault]

I'm attaching a patch that implements SetProjection. Note that the papszCreationOption must contain the 'ICOORDS=' setting that matches the SetProjection string (ICOORDS=G for latlong WGS84, ICOORDS=N for north WGS UTM projection, ICOORDS=S for south WGS UTM projection), so the above test program must be fixed (cf ICOORDS=S and SetUTM(11, TRUE)) The patch also checks in CreateCopy that the ellipsoid is WGS84. It also fixes SetGeoTransform implementation that currently doesn't save the input geotransform in adfGeoTransform member.

[Even Rouault]

This is a duplicate of #753

[warmerdam]

Even,

This looks fine, go ahead and apply in trunk.

I'd like to see a test case in the test suite covering this case. If you aren't comfortable preparing it, let me know and I'll take a crack at it. If that goes well, I wonder if we should risk porting this into 1.5 branch in time for 1.5.1?

[Even Rouault]

Code commited in trunk in r13729 and test added in nitf.py in r13730.