[R-sig-Geo] towgs snd prj files

[Agustin Lobo]

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