[gdal-dev] Fwd: OSRExportToProj4 losing datum

[Roger Bivand]

On Sat, 3 Nov 2012, Even Rouault wrote:

> Selon Roger Bivand <Roger.Bivand at nhh.no>:
>
>> On Sat, 3 Nov 2012, Even Rouault wrote:
>>
>>> 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.
>>
>> Thanks, Even.
>>
>> It isn't easy to make design decisions, and in the case of sp, the R
>> spatial class definition package, our concern was that spatial data
>> analysts were most often unaware of spatial reference issues. PROJ.4
>> strings are a low-threshold text string entry point for users who have no
>> grasp of other representations, and when we began (2003), were among the
>> only games in town. Inside R, the objects cannot reliably be pointers to
>> C/C++ objects cross-platform; indeed, it is only a few month since we
>> achieved PROJ.4/GDAL CRAN (archive network) support for OSX on two intel
>> architectures. I'll try out adding a user-settable argument in calling
>> functions to set and unset the environment variable, but this isn't very
>> satisfactory.
>
> If you were happy with GDAL < 1.8 behaviour, you could also always set
> OVERRIDE_PROJ_DATUM_WITH_TOWGS84=NO before calling ExportToProj4(). Only the
> cases where a TOWGS84 node is defined in the OGRSpatialReference object that
> contradicts the default +towgs84 hard-coded in the datum definition in PROJ.4
> pj_datums.c would be affected, but that's probably a minor use case.

I've committed a user argument to relevant functions in rgdal for setting 
the behaviour to datum-preserving, either on a case-by-case or global 
level. If the environment variable is present, it will have precedence and 
will not be overwritten.

>
>>
>> I agree that in a better world, in sp 2.0, we should replace the
>> representation in the CRS class with (maybe) WKT. But I don't think that
>> this will happen any time soon, and in any case, we would still need to be
>> able to let users to enter the simpler PROJ.4 representations, as writing
>> WKT is much more error-prone.
>
> Using EPSG codes can also be a convenient way for users to specify a SRS, but it
> does not cover all possible SRS.

We use EPSG to look for PROJ.4, but which then retain the +init= tag 
within R.

Best wishes,

Roger

>
>>
>> Best wishes,
>>
>> Roger
>>
>>>
>>> Best regards,
>>>
>>> Even
>>>
>>
>> --
>> Roger Bivand
>> Department of Economics, NHH Norwegian School of Economics,
>> Helleveien 30, N-5045 Bergen, Norway.
>> voice: +47 55 95 93 55; fax +47 55 95 95 43
>> e-mail: Roger.Bivand at nhh.no
>>
>>
>
>
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no