[R-sig-Geo] towgs snd prj files

[Roger Bivand]

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