+towgs84 in st_write

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

+towgs84 in st_write

Manuel Schneider
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@.../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?

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
[hidden email]<mailto:[hidden email]>
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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: +towgs84 in st_write

edzer


On 12/12/2017 03:16 PM, [hidden email] 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@.../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
> [hidden email]<mailto:[hidden email]>
> 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
> [hidden email]
> 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: +towgs84 in st_write

Roger Bivand
Administrator
On Tue, 12 Dec 2017, Edzer Pebesma wrote:

>
>
> On 12/12/2017 03:16 PM, [hidden email] 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@.../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
>> [hidden email]<mailto:[hidden email]>
>> 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
>> [hidden email]
>> 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
> [hidden email]
> 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: [hidden email]
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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: +towgs84 in st_write

edzer


On 12/12/2017 11:27 PM, Roger Bivand wrote:
> Does sf use morphToESRI()?

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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: +towgs84 in st_write

Roger Bivand
Administrator
On Tue, 12 Dec 2017, Edzer Pebesma wrote:

>
>
> On 12/12/2017 11:27 PM, Roger Bivand wrote:
>> Does sf use morphToESRI()?
>
> No.
>

I think the "ESRI Shapefile" driver changed - at one stage it was needed
during writing some time ago, I think. Now the driver simply does it
itself (line 763 in ogrsf_frmts/ogrshapedatasource.cpp:

* -------------------------------------------------------------------- */
/*      Create the .prj file, if required.                              */
/* -------------------------------------------------------------------- */
     if( poSRS != NULL )
     {
         CPLString osPrjFile =
             CPLFormFilename( NULL, pszFilenameWithoutExt, "prj");

         // The shape layer needs its own copy.
         poSRS = poSRS->Clone();
         poSRS->morphToESRI();

         char *pszWKT = NULL;
         VSILFILE *fp = NULL;
         if( poSRS->exportToWkt( &pszWKT ) == OGRERR_NONE
             && (fp = VSIFOpenL( osPrjFile, "wt" )) != NULL )
         {
             VSIFWriteL( pszWKT, strlen(pszWKT), 1, fp );
             VSIFCloseL( fp );
         }

         CPLFree( pszWKT );

         poSRS->morphFromESRI();
     }

So the driver is doing what it believes ArcGIS would like to read - the
*.prj file isn't well specified. In https://issues.qgis.org/issues/2154 
Frank Warmerdam wrote eight years ago: "Yes, OGR does strip the towgs84
parameter when writing the .prj file. TOWGS84 is not a legal construct in
an ESRI Projection Engine string (for .prj files)."

Sounds like another reason to abandon shapefiles as not fit for purpose.

Roger

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: +towgs84 in st_write

Manuel Schneider
In reply to this post by Manuel Schneider
Dear Roger and Edzer

Many thanks for looking into this. Projection is perfectly transferred from R to QGIS when using a geopackage, so this is definitively a shapefile issue. It's absolutely time to change ... and since st_write to GPKG works perfectly now (thanks to your efforts) there is no obstacle to do so.
All best wishes
Manuel

-----Ursprüngliche Nachricht-----
Von: R-sig-Geo [mailto:[hidden email]] Im Auftrag von Roger Bivand
Gesendet: Mittwoch, 13. Dezember 2017 07:09
An: Edzer Pebesma <[hidden email]>
Cc: [hidden email]
Betreff: Re: [R-sig-Geo] +towgs84 in st_write

On Tue, 12 Dec 2017, Edzer Pebesma wrote:

>
>
> On 12/12/2017 11:27 PM, Roger Bivand wrote:
>> Does sf use morphToESRI()?
>
> No.
>

I think the "ESRI Shapefile" driver changed - at one stage it was needed
during writing some time ago, I think. Now the driver simply does it
itself (line 763 in ogrsf_frmts/ogrshapedatasource.cpp:

* -------------------------------------------------------------------- */
/*      Create the .prj file, if required.                              */
/* -------------------------------------------------------------------- */
     if( poSRS != NULL )
     {
         CPLString osPrjFile =
             CPLFormFilename( NULL, pszFilenameWithoutExt, "prj");

         // The shape layer needs its own copy.
         poSRS = poSRS->Clone();
         poSRS->morphToESRI();

         char *pszWKT = NULL;
         VSILFILE *fp = NULL;
         if( poSRS->exportToWkt( &pszWKT ) == OGRERR_NONE
             && (fp = VSIFOpenL( osPrjFile, "wt" )) != NULL )
         {
             VSIFWriteL( pszWKT, strlen(pszWKT), 1, fp );
             VSIFCloseL( fp );
         }

         CPLFree( pszWKT );

         poSRS->morphFromESRI();
     }

So the driver is doing what it believes ArcGIS would like to read - the
*.prj file isn't well specified. In https://issues.qgis.org/issues/2154 
Frank Warmerdam wrote eight years ago: "Yes, OGR does strip the towgs84
parameter when writing the .prj file. TOWGS84 is not a legal construct in
an ESRI Projection Engine string (for .prj files)."

Sounds like another reason to abandon shapefiles as not fit for purpose.

Roger

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo