[R-sig-Geo] +towgs84 in st_write

[Edzer Pebesma]

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:

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