[gdal-dev] Comparing two SRS's, different Proj4

[Even Rouault]

Le vendredi 11 mars 2016 10:58:49, Rutger a écrit :
> Dear list,
> 
> I am trying to compare 2 different projections using OSR in Python. One SRS
> is initialized from a proj4 string given by PostGIS. It's a custom
> definition, so using EPSG is unfortunately not an option.
> 
> The other SRS is read from a Geotiff raster using GDAL, so initialized via
> Wkt from ds.GetProjection().
> 
> My method of comparing is maybe a bit naive, i tried:
> srs1.ExportToProj4() == srs2.ExportToProj4()
> 
> If there something better i'd be happy to hear.
> 
> The problem is however that both ExportToWkt() statements return a
> different string. The one initialized from Wkt defines the datum as
> "+datum=WGS84", which is identical to running "gdalinfo -proj4" and i
> assume this is correct.
> 
> The other srs changes this to "+ellps=WGS84 +towgs84=0,0,0,0,0,0,0".
> 
> Here is a notebook replicating the behavior, but without using a raster:
> https://gist.github.com/RutgerK/e4e40b1025f7279af430
> 
> Why is the ExportToProj4() different then the string from which it was
> initialized?

Rutger,

This is a bit tricky. There are 2 things to consider:
- one of the "culprit" is ImportFromProj4(). It uses a proj.4 function that 
"normalizes" the proj.4 parameters, and that may have some influence on +datum 
parameters being appended a +towgs84 clause  or not depending on the proj4. 
version.
I made a fix for the WGS84 case, a few months ago in trunk (so not in 2.0.X), 
so as not to happen the +towgs84=0,0,0 clause :
https://trac.osgeo.org/gdal/changeset/29775
- the other culprit is  ExportToProj4(). Normally the default behaviour of 
ExportToProj4() when it sees a +towgs84 node is to remove the +datum and let 
only +ellps +towgs84. I don't recall precisely the reason for this, but this 
behaviour can be altered by setting the OVERRIDE_PROJ_DATUM_WITH_TOWGS84 config 
option/environmenet variable to NO.

I've tested on GDAL 2.0.X :
$ OVERRIDE_PROJ_DATUM_WITH_TOWGS84=NO gdalsrsinfo '+proj=aea +lat_1=20 
+lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs '

PROJ.4 : '+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 
+datum=WGS84 +units=m +no_defs '

In trunk, no need to add OVERRIDE_PROJ_DATUM_WITH_TOWGS84=NO since due to the 
first fix I mentionned ImportFromProj4() of +datum=WGS84 no longer generates 
+towgs84=0,0,0 node.

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com