[PATCH] EPSG:2163 incorrectly defined in some instances

[Kyle Shannon]

When applying the EPSG:2163 spatial reference to a file via

$ gdal_translate -a_srs EPSG:2163 in.tif out.tif

The PROJECTION value and the PARAMETER values are not written to the file:

$ gdalinfo -noct out.tif Driver: GTiff/GeoTIFF
Files: out.tif
Size is 4587, 2889
Coordinate System is:
PROJCS["US National Atlas Equal Area",
    GEOGCS["Unspecified datum based upon the Clarke 1866 Authalic Sphere",
        DATUM["Not_specified_based_on_Clarke_1866_Authalic_Sphere",
            SPHEROID["Clarke 1866 Authalic Sphere",6370997,0,
                AUTHORITY["EPSG","7052"]],
            AUTHORITY["EPSG","6052"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4052"]],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","2163"]]
Origin = (-2050500.000000000000000,752500.000000000000000)
Pixel Size = (1000.000000000000000,-1000.000000000000000)
...
Upper Left  (-2050500.000,  752500.000) 
Lower Left  (-2050500.000,-2136500.000) 
Upper Right ( 2536500.000,  752500.000) 
Lower Right ( 2536500.000,-2136500.000) 
Center      (  243000.000, -692000.000) 
...

When I run:

$ testepsg EPSG:2163
Validate Succeeds.
WKT[EPSG:2163] =
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6370997,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",45],
    PARAMETER["longitude_of_center",-100],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1],
    AUTHORITY["EPSG","2163"]]
...
PROJ.4 rendering of [EPSG:2163] = +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs 
...

It properly reports the values. The tif cannot be warped as there is no way to get a proj4 rendering:

$ gdalwarp -t_srs EPSG:4326 out.tif warp.tif
Copying color table from out.tif to new file.
ERROR 1: No PROJ.4 translation for source SRS, coordinate
transformation initialization has failed.
Creating output file that is 4587P x 2889L.
Processing input file out.tif.
ERROR 1: No PROJ.4 translation for source SRS, coordinate
transformation initialization has failed.
0...10...20...30...40...50...60...70...80...90...100 - done.

If I define the source srs as a proj string, it works:

$ gdalwarp -s_srs "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs" -t_srs EPSG:4326 out.tif warp.tif
Copying color table from out.tif to new file.
Creating output file that is 5791P x 2687L.
Processing input file out.tif.
0...10...20...30...40...50...60...70...80...90...100 - done.

I don't know if this a problem with the EPSG csv files, or something else. I would imagine testepsg and gdal_translate use SetFromUserInput() function to populate the OGRSpatialReference.

I believe that this can be reproduced by translating any tif and applying the srs.

[Even Rouault]

I'd suggest the following patch (one piece in ogr_fromepsg.cpp, the other one in geo_normalize.c). The key is that EPSG:2163, 3408, 3409, 3973 and 3974 uses EPSG projection method 1027, that seems to be a variant of 9820 (LAEA) according to EPSG Guidance 7-2 paragraph 1.3.11.1. It needs review from someone more versed in projections than me though to be sure that 1027 can be treated as an alias of 9820.

[Even Rouault]

Patch to OSR and libgeotiff to deal with EPSG projection method 1027

[Kyle Shannon]

I applied the patch and ran the commands again. With CPL_DEBUG set to ON, her is the output of gdal_translate:

$ gdal_translate -a_srs EPSG:2163 in.tif out.tif --config CPL_DEBUG ON
EPSG: No WKT support for projection method 1027.
OGRCT: PROJ >= 4.8.0 features enabled
GDAL: GDALOpen(in.tif, this=0x6bd250) succeeds as GTiff.
Input file size is 4587, 2889
GDAL: GDALDefaultOverviews::OverviewScan()
GDAL: QuietDelete(out.tif) invoking Delete()
GDAL: GDALOpen(out.tif, this=0x6cbad0) succeeds as GTiff.
GDAL: GDALDefaultOverviews::OverviewScan()
GDAL: GDALClose(out.tif, this=0x6cbad0)
0GDAL: GDALOpen(GTIFF_RAW:out.tif, this=0x6cba00) succeeds as GTiff.
GDAL: GDALDatasetCopyWholeRaster(): 4587*2176 swaths, bInterleave=0
...10...20...30...40...50...60...70...80...90...100 - done.
GDAL: GDALClose(out.tif, this=0x6cba00)
GDAL: GDALClose(, this=0x6bd700)
GDAL: GDALClose(in.tif, this=0x6bd250)
GDAL: GDALDeregister_GTiff() called.
kyle@Lucky13:~/Desktop/jolly_gdal$ gdal_translate -a_srs EPSG:1027 in.tif out.tif --config CPL_DEBUG ON
OGRCT: PROJ >= 4.8.0 features enabled
ERROR 6: EPSG PCS/GCS code 1027 not found in EPSG support files.  Is this a valid
EPSG coordinate system?
Failed to process SRS definition: EPSG:1027
GDAL: GDALDeregister_GTiff() called.

EPSG: throws a few lines to the DEBUG handler. And gdalinfo reports the same system:

$ gdalinfo out.tif -noct
Driver: GTiff/GeoTIFF
Files: out.tif
Size is 4587, 2889
Coordinate System is:
PROJCS["US National Atlas Equal Area",
    GEOGCS["Unspecified datum based upon the Clarke 1866 Authalic Sphere",
        DATUM["Not_specified_based_on_Clarke_1866_Authalic_Sphere",
            SPHEROID["Clarke 1866 Authalic Sphere",6370997,0,
                AUTHORITY["EPSG","7052"]],
            AUTHORITY["EPSG","6052"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4052"]],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","2163"]]
Origin = (-2050500.000000000000000,752500.000000000000000)
Pixel Size = (1000.000000000000000,-1000.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-2050500.000,  752500.000) 
Lower Left  (-2050500.000,-2136500.000) 
Upper Right ( 2536500.000,  752500.000) 
Lower Right ( 2536500.000,-2136500.000) 
Center      (  243000.000, -692000.000) 
Band 1 Block=4587x1 Type=Byte, ColorInterp=Palette
  Color Table (RGB with 256 entries)

Is there a problem elsewhere, or does 1027 need to be aliased to 9820 somewhere else? Thanks.

[Even Rouault]

The warning 'EPSG: No WKT support for projection method 1027.' comes from ogr/ogr_fromepsg.cpp and is a clear sign that something went wrong when you've patched your GDAL version or that you are not running the GDAL version you believe... And I'd note that 1027 is an EPSG projection method, not a EPSG CRS code, so you shouldn't try to use it directly. Anyway don't bet too much on the proposed patch. It could be that 9821 is the more appropriate alias. I just don't know...

[Kyle Shannon]

Of course you were right, I had switched my prefix to do some debugging and was using my system install to run tests on this ticket. Sorry about that. It looks like it worked, and qgis picks up the srs, which it did not before. I have heard complaints about that SRS from people who use esri software as well. I guess if a projection person checks this off, I feel good about it. Thanks again Even.

[warmerdam]

I am concerned that 1027 is specifically spherical LAEA and we are losing that distinction, but since EPSG:2163 uses a spherical geographic coordinate system perhaps it isn't important. I am applying the patch in libgeotiff (r 2037) and GDAL trunk (r22424). Also pulled libgeotiff changes down into GDAL trunk (r22424).