+towgs84 in st_write

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

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

-- 
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081