On Sat, 3 Nov 2012, Even Rouault wrote:
Selon Roger Bivand <roger.biv...@nhh.no>:
On Sat, 3 Nov 2012, Even Rouault wrote:
Selon Roger Bivand <roger.biv...@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.biv...@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.biv...@nhh.no
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
