[R-sig-Geo] +towgs84 in st_write

[Roger Bivand]

On Tue, 12 Dec 2017, Edzer Pebesma wrote:

>
>
> On 12/12/2017 03:16 PM, manuel.schneider at agroscope.admin.ch wrote:
>> Dear list
>> 
>> I have a feature projected in EPSG:2056 (http://spatialreference.org/ref/epsg/2056/). When st_write this to a shapefile I get a prj-File without the +towgs84 parameters and this results in an offset if the shapefile is projected on the fly in a QGIS project with OpenLayers (EPSG:3857). Placement is correct if I open in a project with EPSG:2056 or if I manually assign EPSG:2056 to the shapefile.
>> This looks a bit similar to this older thread https://www.mail-archive.com/r-sig-geo@r-project.org/msg06462.html but in this case the +towgs84 parameters are well defined.
>> The feature has
>>> st_crs(p)
>> Coordinate Reference System:
>>   EPSG: 2056
>>   proj4string: "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"
>> st_write generates a prj reading PROJCS["Hotine_Oblique_Mercator_Azimuth_Center",GEOGCS["GCS_Bessel 1841",DATUM["D_unknown",SPHEROID["bessel",6377397.155,299.1528128]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],PARAMETER["latitude_of_center",46.95240555555556],PARAMETER["longitude_of_center",7.439583333333333],PARAMETER["azimuth",90],PARAMETER["scale_factor",1],PARAMETER["false_easting",2600000],PARAMETER["false_northing",1200000],UNIT["Meter",1]]
>> The result is that QGIS interprets this as a user defined CRS with +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +units=m +no_defs . The +towgs84 are assumed to be 0, I guess, and this results in the offset.
>> 
>> Does anybody know if there is a way to specify the information written by st_write to the .prj-File?
>
> I can reproduce that, and see the same when writing to ESRI Shapefile
> using rgdal::writeOGR. The funny thing is that from gdal one also gets
> the following back, when asking for proj.4 -> wkt conversion:

I think this is the writeOGR() morphToESRI argument (but commented out, 
maybe in the "ESRI Shapefile" driver?). Does sf use morphToESRI()? This 
changes the WKT SRS representation:

library(rgdal)
CRSargs(CRS("+init=epsg:2056"))
showWKT(CRSargs(CRS("+init=epsg:2056")), morphToESRI=TRUE)
showWKT(CRSargs(CRS("+init=epsg:2056")), morphToESRI=FALSE)

Roger

>
> cat(st_as_text(st_crs(2056), pretty = TRUE))
> PROJCS["unnamed",
>    GEOGCS["Bessel 1841",
>        DATUM["unknown",
>            SPHEROID["bessel",6377397.155,299.1528128],
>            TOWGS84[674.374,15.056,405.346,0,0,0,0]],
>        PRIMEM["Greenwich",0],
>        UNIT["degree",0.0174532925199433]],
>    PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
>    PARAMETER["latitude_of_center",46.95240555555556],
>    PARAMETER["longitude_of_center",7.439583333333333],
>    PARAMETER["azimuth",90],
>    PARAMETER["rectified_grid_angle",90],
>    PARAMETER["scale_factor",1],
>    PARAMETER["false_easting",2600000],
>    PARAMETER["false_northing",1200000],
>    UNIT["Meter",1]]
>
> which has the towgs parameters.
>
> If I write this file to GeoPackage, i.e. by
>
> st_write(nc, "xx.gpkg")
>
> I do see the towgs parameters when looking at it with ogrinfo (or
> reading it back with st_read).
>
> One more reason to ditch shapefiles?
>
>> 
>> Many thanks for suggestions
>> Manuel
>> 
>> 
>> 
>> -----
>> Manuel Schneider (Ing.-Agr. ETH, Dr. sc. nat.)
>> Team leader, Mountain grassland ecology & management
>> 
>> Federal Department of Economic Affairs, Education and Research EAER
>> Agroscope
>> 
>> Reckenholzstrasse 191, CH-8046 Z???rich
>> Ph. +41 58 468 75 98
>> Fax  +41 58 468 72 01
>> manuel.schneider at agroscope.admin.ch<mailto:manuel.schneider at art.admin.ch>
>> https://www.agroscope.admin.ch/agroscope/en/home/topics/plant-production/forage-grassland-grazing-systems.html
>> http://scholar.google.com/citations?user=a0CE8xoAAAAJ&hl=de
>> 
>> www.agroscope.ch<http://www.agroscope.ch/>  I good food, healthy environment
>> 
>> 
>>
>> 	[[alternative HTML version deleted]]
>> 
>> 
>> 
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> 
>
> -- 
> Edzer Pebesma
> Institute for Geoinformatics
> Heisenbergstrasse 2, 48151 Muenster, Germany
> Phone: +49 251 8333081
>
> _______________________________________________
> 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, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand at nhh.no
Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en