towgs snd prj files

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.

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

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

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