[gdal-dev] Fwd: OSRExportToProj4 losing datum

[Even Rouault]

Selon Roger Bivand <Roger.Bivand at nhh.no>:

> Even,
>
> I've opened #4880 on this - it is problematic in the R-spatial setting
> because the CRS (coordinate reference system) object recorded in each
> Spatial object uses the PROJ.4 string to represent spatial reference. A user
> can create a CRS, write out a raster (which expands the description to
> include , read it back, and the CRS are not the same. The same does not
> appear to happen with OGR write/read (ESRI Shapefile driver). Is the
> conclusion that calling programs must first find "NAD83" in the raster input
> WKT, if it is there set the environment variable if the default GDAL
> behaviour is not desired, export to Proj4, then remove the environment
> variable if set? Does this only affect NAD83 - it does seem to be limited to
> this datum, as WGS84 does not show the same behaviour.

Roger,

Hum, I somehow regret of having tried to fix the issue that was originally
raised in trac.osgeo.org/gdal/ticket/3450 since it seems that all GRASS and now
RGdal developers hate me now ;-) But hopefully they will realize that they
should consider using more reliable ways of representing SRS.

In the case of NAD83, +datum=nad83 is expanded as
GEOGCS[\"NAD83\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS
1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6269\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9108\"]],AUTHORITY[\"EPSG\",\"4269\"]]
. This expansion is done at import time (importFromProj4).

This explains why at export time the +datum isn't exported since there's a
TOWGS84 node.

You would also have a similar behaviour with +datum=WGS72 since +datum=WGS72
expansion implicitely adds a TOWGS84 node.

In the case of NAD83, I'm not sure why SetWellKnownGeogCS("NAD83") adds the
TOWGS84[0,0,0,0,0,0,0]. If this was removed, then I believe that exportToProj4()
would again generate +datum=nad83. But I'm now more than hesitant in touching
anything in that area, since this could potentially break something else.

Sorry to repeat myself, but you should NOT use proj4 strings as a way of
representing a SRS. It is just not reliable. It is not, and has never been. You
must realize that proj4 only knows 10 datum names (see pj_datums.c). All other
datums must be translated into +ellps= +towgs84= anyway. You shouldn't assume
that round-tripping importFromProj4() / exportToProj4() will be perfect. proj4
strings are only usefull for proj4, to do SRS transformations. They are not
meant at being a comprehensive way of representing SRS. I'm not sure how hard it
is for RGdal, but you should consider migrating to using WKT strings / model.
WKT has been designed for that.

Best regards,

Even