[R-sig-Geo] towgs snd prj files

Wed Nov 7 16:38:26 CET 2012

On Wed, 7 Nov 2012, Agustin Lobo wrote:

> ok, thanks, just got your message while I was typing my second one. It 
> seems we have to live with this issue. The best is probably make a 
> version of the object in both datums via spTransform() and then write 2 
> shapefiles as well.

The issue is that where the +towgs84 coefficients are not known with 
complete certainty (the case for ED50, where they differ a good deal 
depending on where you are), they are dropped when converting from PROJ.4 
to the internal representation used by OGR:

library(rgdal)
pt <- matrix(c(10000, 3000000), ncol=2)
spdf <- SpatialPointsDataFrame(SpatialPoints(pt), data=data.frame(ID=1))
proj4string(spdf) <- CRS("+init=epsg:23031")
proj4string(spdf)
td <- tempdir()
writeOGR(spdf, td, "tmp", driver="ESRI Shapefile") 
readLines(paste(td, "tmp.prj", sep="/"))
proj4string(readOGR(td, "tmp"))

If we use OSGB instead, with D_OSGB_1936 rather than D_unknown for its 
datum, the +towgs84 values are preserved:

spdf <- SpatialPointsDataFrame(SpatialPoints(pt), data=data.frame(ID=1))
proj4string(spdf) <- CRS("+init=epsg:27700")
proj4string(spdf)
writeOGR(spdf, td, "tmposgb", driver="ESRI Shapefile")
readLines(paste(td, "tmposgb.prj", sep="/"))
proj4string(readOGR(td, "tmposgb"))

So the logic in OGR is trying to ensure that only information that is 
authoritative is encoded in the output file. QGIS has a "custom" way of 
handling spatial reference, trying to make things "easy" for the user, 
whereas OGR prefers not to guess when it doesn't have authority. ED50 is 
not a standard at all, it is a collection of national (sometimes 
sub-national) standards that have also changed over time. It is only since 
GRS80/WGS84 that we have satellite-based measurements on which to base 
spheroid and datum definitions, hence the importance of having authority 
for transformations to and from WGS84. It doesn't make matters easier that 
the Earth also changes, but these changes are small in comparison with the 
differences induced by poor datum transformation. The datum values 
recognised by PROJ.4 are:

projInfo("datum")

which includes OSGB36, but not ED50. This:

http://www.ogp.org.uk/pubs/373-10.pdf

gives a flavour of the difficulties involved in oil exploration in the 
North Sea when also transforming between ED50 and WGS84 - note that these 
coefficients are very different from those shown in another source:

http://earth-info.nga.mil/GandG/coordsys/onlinedatum/CountryEuropeTable.html

referred to in this thread:

http://lists.osgeo.org/pipermail/gdal-dev/2009-May/020656.html

This is a never-ending story ...

Roger

>
> Or perhaps using spatialite objects instead of shapfiles would be a 
> solution? 
> http://www.gaia-gis.it/spatialite-2.0/SpatiaLite2-tutorial.html#t5

No, nothing to do with this.

>
> Agus
>
> On Wed, Nov 7, 2012 at 1:21 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
>> On Wed, 7 Nov 2012, Agustin Lobo wrote:
>>
>>> This might be a bit off-topic.
>>
>>
>> There is no guarantee of any kind that GDAL/OGR will reproduce the same
>> projection parameters on re-reading, and there is nothing in sp/rgdal that
>> can be done about this. The thread at:
>>
>> https://stat.ethz.ch/pipermail/r-sig-geo/2012-November/016581.html
>>
>> did permit an extra argument to be passed to readGDAL() to preserve datum in
>> addition to equivalent ellps and towgs84 tags.
>>
>> The underlying difficulties are that none of the various spatial reference
>> systems are unequivocal, and changing between PROJ.4, WKT, and others (EPSG
>> codes in different versions) may drop tags, especially where they are not
>> uniquely defined.
>>
>> On R-forge (revision 383) I have added a morphToESRI= logical argument to
>> writeOGR(), NULL default, set automatically TRUE if the driver is "ESRI
>> Shapefile" and FALSE for all other drivers (this is the same as internal
>> code has always been). The user may try to see whether this makes a
>> difference, I didn't see any difference in the exported WKT written in the
>> *.prj file.
>>
>> The main issue is that for any given specification, the towgs84 tag may not
>> be unique (or other tags), so that the PROJ.4/GDAL policy is not to do
>> anything that might not be accurate in conversion. There are no easy
>> solutions here.
>>
>> In your case, for the copy og EPSG distributed with rgdal in the binary
>> Windows and OSX packages, and for my PROJ.4/GDAL combination:
>>
>>> CRS("+init=epsg:23031")
>>
>> CRS arguments:
>>  +init=epsg:23031 +proj=utm +zone=31 +ellps=intl
>> +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs
>>
>> but http://spatialreference.org/ref/epsg/23031/ does not include a +towgs84
>> term. Which is credible? This is an ED50 projection/datum, so I'd actually
>> trust dropping towgs84 more than keeping it. The web page suggests that the
>> ED50 datum should be defined by country with respect to towgs84; I often
>> look at Grids & Datums to cross-check:
>>
>> http://www.asprs.org/Grids-Datums.html
>>
>> Hope this helps,
>>
>> Roger
>>
>>> I do:
>>>>
>>>> wpBert700 at proj4string
>>>
>>> CRS arguments:
>>> +init=epsg:23031 +proj=utm +zone=31 +ellps=intl
>>> +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs
>>>>
>>>> writeOGR(wpBert700,dsn="wpBert700v2",layer="wpBert700v2",driver="ESRI
>>>> Shapefile",overwrite=T)
>>>
>>>
>>> But the created prj file:
>>> PROJCS["UTM_Zone_31_Northern_Hemisphere",GEOGCS["GCS_International
>>> 1909
>>> (Hayford)",DATUM["D_unknown",SPHEROID["intl",6378388,297]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]
>>>
>>> lacks the +towgs84=-87,-98,-121,0,0,0,0 information, which is a
>>> problem if I do not remember I have to add it for a GIS display
>>> with reprojection on the fly.
>>> Is there any vector format that can  be written by rgdal in which the
>>> CRS info would include the +towgs84=-87,-98,-121,0,0,0,0 ?
>>>
>>> Thanks!
>>>
>>> Agus
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>> --
>> 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
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
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