[gdal-dev] OGRSpatialReference::IsSame returning TRUE unexpectedly

[Matthew Amato]

Short version:

I don't understand why the below code snippet is returning TRUE.  Unless
I'm missing something, EPSG 3857 is not the same as 3395, why is GDAL
reporting that they are?

            auto epsg3857 = OGRSpatialReference();
            epsg3857.importFromEPSG(3857);

            auto epsg3395 = OGRSpatialReference();
            epsg3395.importFromEPSG(3395);

            if (epsg3857.IsSame(&epsg3395) == TRUE) {
                std::cout << "Bug?" << std::endl;
            }

Longer version:

I'm attempting to load in a GeoTIFF and detect whether it is already in
EPSG:3857 or needs to be reprojected.  My function looks like this:

        bool isSpatialReferenceSystem(GDALDataset &dataset, int epsg)
        {
            std::string in_srs_wkt = dataset.GetProjectionRef();

            if (in_srs_wkt.length() == 0)
            {
                return false;
            }

            auto dataset_srs = OGRSpatialReference(in_srs_wkt.c_str());
            auto target_srs = OGRSpatialReference();
            target_srs.importFromEPSG(epsg);
            return dataset_srs.IsSame(&target_srs) == TRUE;
        }

This function worked as expected until a few days ago when I was given a
GeoTIFF with the following information (via gdalinfo):

Driver: GTiff/GeoTIFF
Files: data.tif
Size is 11779, 8832
Coordinate System is:
PROJCS["WGS 84 / World Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","3395"]]
Origin = (723576.685867484660000,5474400.172804606100000)
Pixel Size = (33.077567560570841,-33.077567560572312)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  723576.686, 5474400.173) (  6d30' 0.00"E, 44d15' 0.34"N)
Lower Left  (  723576.686, 5182259.096) (  6d30' 0.00"E, 42d20' 0.00"N)
Upper Right ( 1113197.354, 5474400.173) ( 10d 0' 0.08"E, 44d15' 0.34"N)
Lower Right ( 1113197.354, 5182259.096) ( 10d 0' 0.08"E, 42d20' 0.00"N)
Center      (  918387.020, 5328329.634) (  8d15' 0.04"E, 43d17'57.55"N)
Band 1 Block=11779x1 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=11779x1 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=11779x1 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=11779x1 Type=Byte, ColorInterp=Alpha

Creating a GDALDataset from the image and calling my function returns TRUE
when I expect it to return false, after some debugging, I traced it back to
my simpler code snippet that shows 3857 and 3395 being treated as the same.

The odd thing is that if I tell GDAL to reproject the image to 3857 (via a
VRT), everything works.  So GDAL clearly knows they are different under the
hood.

I assume my understanding of IsSame is flawed so any clarification would be
a big help.

Thanks,

Matt
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20170905/894afe14/attachment.html>