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 _______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
