[gdal-dev] How to check if projections are the same?

[Even Rouault]

Selon Paul Meems <bontepaarden at gmail.com>:

> We use GDAL/OGR/Proj4 in our application: MapWindow GIS.
>
> We're testing the MapWinGIS.GeoProjection class. It's just a wrapper around
> OGRSpatialReference class, so it's GDAL behavior infact.
>
> The test code:
>
>         private void button8_Click(object sender, EventArgs e)
>         {
>             MapWinGIS.GeoProjection proj = new GeoProjection();
>             proj.ImportFromEPSG(4326);
>             System.Diagnostics.Debug.Print(proj.ExportToWKT());
>
>             string s = proj.ExportToProj4();
>             proj.ImportFromProj4(s);
>             System.Diagnostics.Debug.Print(proj.ExportToWKT());
>
>             s = proj.ExportToProj4();
>             proj.ImportFromProj4(s);
>             System.Diagnostics.Debug.Print(proj.ExportToWKT());
>         }
>
>
> Results:
>
> 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"]]
>
>
> 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"]],AUTHORITY["EPSG","4326"]]
>
>
> GEOGCS["WGS
>
84",DATUM["unknown",SPHEROID["WGS84",6378137,298.257223563],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]
>
>
> If we test the first and the third string with GeoProjection.IsSame() it
> will return false, but they are the same.

Looking at  see OGRSpatialReference::IsSameGeogCS(), I see :

/* -------------------------------------------------------------------- */
/*      Does the datum name match?  Note that we assume                 */
/*      compatibility if either is missing a datum.                     */
/* -------------------------------------------------------------------- */
    pszThisValue = this->GetAttrValue( "DATUM" );
    pszOtherValue = poOther->GetAttrValue( "DATUM" );

    if( pszThisValue != NULL && pszOtherValue != NULL
        && !EQUAL(pszThisValue,pszOtherValue) )
        return FALSE;

So clearly the DATUM name is a criteria. In your example WGS_1984 != unknown,
hence the result.

>
> How should we bypass this?

You can reimplement a variant of IsSame() for fit your needs (remove the above
test).

But I'm not too sure if you can just use the other informations (spheroid,
towgs84, prime meridian, etc...) to test if 2 GEOGCS are identical. For example
:

WKT[EPSG:3824] =
GEOGCS["TWD97",
    DATUM["Taiwan_Datum_1997",
        SPHEROID["GRS 1980",6378137,298.257222101,
            AUTHORITY["EPSG","7019"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","1026"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","3824"]]

WKT[EPSG:3889] =
GEOGCS["IGRS",
    DATUM["Iraqi_Geospatial_Reference_System",
        SPHEROID["GRS 1980",6378137,298.257222101,
            AUTHORITY["EPSG","7019"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","1029"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","3889"]]

As far as proj4 is concerned, they are identical ("+proj=longlat +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +no_defs "). But are they the "same"... ?

So it really depends on what you do with the result of the comparision. If it is
just for reprojection purposes, you could ExportToProj4() and compare the
strings. But I guess there would be also some cases where it wouldn't be
perfect. For example "+datum=WGS84" is equivalent with "+ellps=WGS84
+towgs84=0,0,0,0,0,0,0".

Best regards,

Even