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

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

[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


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: OGRSpatialReference::IsSame returning TRUE unexpectedly

Even Rouault-2

> I assume my understanding of IsSame is flawed so any clarification would be

> a big help.

>

 

Matt,

 

Your understanding of IsSame() is correct. You just hit a bug. I've fixed it per https://trac.osgeo.org/gdal/ticket/7029. You can see the changeset code as a workaround to add in your own code if you need.

 

IsSame() is a "smart" comparison, ie it tries to abstract away unsignificant differences like the PROJCS text value that don't participate really to the SRS and can vary over formats, data producers, etc... The issue here is that web-mercator (EPSG:3857) and WGS84-mercator (EPSG:3395) are represented the same way at the GEOGCS and projection levels. The only difference is that GDAL creates a EXTENSION["PROJ4", ...] node in the web-mercator case to really use the spherical definition. And IsSame() before my fix neglected to compare that extension node. So you were a bit unlucky since that must be pretty much the only pair of EPSG entries in that situation.

 

My initial tought was to compare the authority codes/nodes when they exist on both sides, but I was wondering if there are not cases where you'd have the same SRS registered in 2 different catalogs. Presumably within the same catalog, a different code should mean a different CRS, but I'm not even completely sure you couldn't find oddities. I was thinking to EPSG:3857 vs the deprecated EPSG:3785, but this isn't a good example, since they differ on the DATUM. So currently IsSame() returns FALSE on that pair, although I guess one might consider them as equivalent (the deprecated version is just a buggy expression on the current EPSG:3857). Or not...

You also have the situation where really different CRS (IsSame() == FALSE) ends up being translated the same way as proj.4 strings...

 

So it is actually not obvious to define the semantics of IsSame(). It should at least return FALSE for CRS that are clearly different (coordinate transformation results in differences).

 

Even

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: OGRSpatialReference::IsSame returning TRUE unexpectedly

Matthew Amato
Thanks for the detailed write up and quick fix!

Just to clarify your last sentence, are you saying there are still potential cases where IsSame returns TRUE when it shouldn't?  Or were you just emphasizing that this is a very tricky check to get right in all cases?  I'm just wondering if I should ditch my optimization of checking the spatial reference and always do a reproject (via VRT) even if it results in a no-op in some cases.  I assume the overhead of a no-op VRT is fairly small and I'm concerned with correctness more than anything else.

Thanks again.

On Tue, Sep 5, 2017 at 5:02 PM, Even Rouault <[hidden email]> wrote:

> I assume my understanding of IsSame is flawed so any clarification would be

> a big help.

>

 

Matt,

 

Your understanding of IsSame() is correct. You just hit a bug. I've fixed it per https://trac.osgeo.org/gdal/ticket/7029. You can see the changeset code as a workaround to add in your own code if you need.

 

IsSame() is a "smart" comparison, ie it tries to abstract away unsignificant differences like the PROJCS text value that don't participate really to the SRS and can vary over formats, data producers, etc... The issue here is that web-mercator (EPSG:3857) and WGS84-mercator (EPSG:3395) are represented the same way at the GEOGCS and projection levels. The only difference is that GDAL creates a EXTENSION["PROJ4", ...] node in the web-mercator case to really use the spherical definition. And IsSame() before my fix neglected to compare that extension node. So you were a bit unlucky since that must be pretty much the only pair of EPSG entries in that situation.

 

My initial tought was to compare the authority codes/nodes when they exist on both sides, but I was wondering if there are not cases where you'd have the same SRS registered in 2 different catalogs. Presumably within the same catalog, a different code should mean a different CRS, but I'm not even completely sure you couldn't find oddities. I was thinking to EPSG:3857 vs the deprecated EPSG:3785, but this isn't a good example, since they differ on the DATUM. So currently IsSame() returns FALSE on that pair, although I guess one might consider them as equivalent (the deprecated version is just a buggy expression on the current EPSG:3857). Or not...

You also have the situation where really different CRS (IsSame() == FALSE) ends up being translated the same way as proj.4 strings...

 

So it is actually not obvious to define the semantics of IsSame(). It should at least return FALSE for CRS that are clearly different (coordinate transformation results in differences).

 

Even

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com



_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: OGRSpatialReference::IsSame returning TRUE unexpectedly

Even Rouault-2

On mardi 5 septembre 2017 20:21:17 CEST Matthew Amato wrote:

> Thanks for the detailed write up and quick fix!

>

> Just to clarify your last sentence, are you saying there are still

> potential cases where IsSame returns TRUE when it shouldn't?

 

Hopefully not.

 

> I'm just wondering if I should ditch my optimization of checking

> the spatial reference and always do a reproject (via VRT) even if it

> results in a no-op in some cases. I assume the overhead of a no-op VRT is

> fairly small and I'm concerned with correctness more than anything else.

 

If "same for reprojection" is what matters for you, a safe way of checking is calling exportToProj4() on both SRS and checking for equality on the returned strings.

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev