towgs snd prj files
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