Help when moving from PROJ4 to PROJ6
This seems to have been resolved now by reinstalling sf. Re: axis order: sf will not swap coordinates by itself, it is only the way x/y coordinates are interpreted that changes (for EPSG:4326 as lng/lat by default, but as lat/lng after giving st_axis_order(TRUE)). You need to do something to get them interpreted as lat/lng.
On 7/28/20 3:34 PM, Roger Bivand wrote:
On Tue, 28 Jul 2020, Gilberto Camara wrote:
Dear R-SIG-GEO (esp. Roger and Edzer)
I am having problems with ?raster?, ?rgdal?, and ?sf? when moving from
PROJ4 to PROJ6. Roger?s explanation on the r-spatial blog
("https://www.r-spatial.org/r/2020/03/17/wkt.html?) and the rgdal blog
(http://rgdal.r-forge.r-project.org/articles/PROJ6_GDAL3.html
<http://rgdal.r-forge.r-project.org/articles/PROJ6_GDAL3.html>) are
clear on the reasons why there are problems with CRS. As Roger points
out, ?package maintainers will need to review their use of crs
objects?. So very true!
I appreciate the tremendous effort of Edzer and Roger for moving from
PROJ4 to PROJ6.
However, I am failing to understand what is happening in the following
example:
================================================
library(sf)
Could you please report: sf_extSoftVersion() Mine are (on Linux, own PROJ & GDAL installations): sf_extSoftVersion() ????????? GEOS?????????? GDAL???????? proj.4 GDAL_with_GEOS???? USE_PROJ_H ?????? "3.8.1"??????? "3.1.2"??????? "7.1.0"???????? "true"???????? "true" where there is no change in the coordinates (there should not be a change as sf::st_crs() only sets the CRS (but may swap the axes as per specification). I see:
st_crs(ll_sfc2)
Coordinate Reference System: ? User input: +proj=longlat +datum=WGS84 +no_defs ? wkt: GEOGCRS["unknown", ??? DATUM["World Geodetic System 1984", ??????? ELLIPSOID["WGS 84",6378137,298.257223563, ??????????? LENGTHUNIT["metre",1]], ??????? ID["EPSG",6326]], ??? PRIMEM["Greenwich",0, ??????? ANGLEUNIT["degree",0.0174532925199433], ??????? ID["EPSG",8901]], ??? CS[ellipsoidal,2], ??????? AXIS["longitude",east, ??????????? ORDER[1], ??????????? ANGLEUNIT["degree",0.0174532925199433, ??????????????? ID["EPSG",9122]]], ??????? AXIS["latitude",north, ??????????? ORDER[2], ??????????? ANGLEUNIT["degree",0.0174532925199433, ??????????????? ID["EPSG",9122]]]] but
st_crs(ll_sfc)
Coordinate Reference System: ? User input: EPSG:4326 ? wkt: GEOGCRS["WGS 84", ??? DATUM["World Geodetic System 1984", ??????? ELLIPSOID["WGS 84",6378137,298.257223563, ??????????? LENGTHUNIT["metre",1]]], ??? PRIMEM["Greenwich",0, ??????? ANGLEUNIT["degree",0.0174532925199433]], ??? CS[ellipsoidal,2], ??????? AXIS["geodetic latitude (Lat)",north, ??????????? ORDER[1], ??????????? ANGLEUNIT["degree",0.0174532925199433]], ??????? AXIS["geodetic longitude (Lon)",east, ??????????? ORDER[2], ??????????? ANGLEUNIT["degree",0.0174532925199433]], ??? USAGE[ ??????? SCOPE["unknown"], ??????? AREA["World"], ??????? BBOX[-90,-180,90,180]], ??? ID["EPSG",4326]] I can't see from https://www.r-spatial.org/r/2020/03/17/wkt.html#axis-order how to instantiate a crs object with non-authorized order. In sp/rgdal:
sp <- SpatialPoints(cbind(longitude, latitude),
proj4string=CRS("+proj=longlat +datum=WGS84"))
sp
SpatialPoints: ???? longitude? latitude [1,] -55.62277 -11.75932 Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84 +no_defs
cat(wkt(sp), "\n")
GEOGCRS["unknown", ??? DATUM["World Geodetic System 1984", ??????? ELLIPSOID["WGS 84",6378137,298.257223563, ??????????? LENGTHUNIT["metre",1]], ??????? ID["EPSG",6326]], ??? PRIMEM["Greenwich",0, ??????? ANGLEUNIT["degree",0.0174532925199433], ??????? ID["EPSG",8901]], ??? CS[ellipsoidal,2], ??????? AXIS["longitude",east, ??????????? ORDER[1], ??????????? ANGLEUNIT["degree",0.0174532925199433, ??????????????? ID["EPSG",9122]]], ??????? AXIS["latitude",north, ??????????? ORDER[2], ??????????? ANGLEUNIT["degree",0.0174532925199433, ??????????????? ID["EPSG",9122]]]]
sp2 <- SpatialPoints(cbind(longitude, latitude),
proj4string=CRS(SRS_string="EPSG:4326"))
sp2
SpatialPoints: ???? longitude? latitude [1,] -55.62277 -11.75932 Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84 +no_defs
cat(wkt(sp2), "\n")
GEOGCRS["WGS 84", ??? DATUM["World Geodetic System 1984", ??????? ELLIPSOID["WGS 84",6378137,298.257223563, ??????????? LENGTHUNIT["metre",1]], ??????? ID["EPSG",6326]], ??? PRIMEM["Greenwich",0, ??????? ANGLEUNIT["degree",0.0174532925199433], ??????? ID["EPSG",8901]], ??? CS[ellipsoidal,2], ??????? AXIS["longitude",east, ??????????? ORDER[1], ??????????? ANGLEUNIT["degree",0.0174532925199433, ??????????????? ID["EPSG",9122]]], ??????? AXIS["latitude",north, ??????????? ORDER[2], ??????????? ANGLEUNIT["degree",0.0174532925199433, ??????????????? ID["EPSG",9122]]], ??? USAGE[ ??????? SCOPE["unknown"], ??????? AREA["World"], ??????? BBOX[-90,-180,90,180]]] where sp/rgdal enforces legacy/GIS/visualization axis order. Coercing from the EPSG:4326 sp/rgdal object to sfc gives the intended order:
ll_sfc3 <- st_as_sfc(sp2) ll_sfc3
Geometry set for 1 feature geometry type:? POINT dimension:????? XY bbox:?????????? xmin: -55.62277 ymin: -11.75932 xmax: -55.62277 ymax: -11.75932 geographic CRS: WGS 84 POINT (-55.62277 -11.75932)
st_crs(ll_sfc3)
Coordinate Reference System: ? User input: WGS 84 ? wkt: GEOGCRS["WGS 84", ??? DATUM["World Geodetic System 1984", ??????? ELLIPSOID["WGS 84",6378137,298.257223563, ??????????? LENGTHUNIT["metre",1]], ??????? ID["EPSG",6326]], ??? PRIMEM["Greenwich",0, ??????? ANGLEUNIT["degree",0.0174532925199433], ??????? ID["EPSG",8901]], ??? CS[ellipsoidal,2], ??????? AXIS["longitude",east, ??????????? ORDER[1], ??????????? ANGLEUNIT["degree",0.0174532925199433, ??????????????? ID["EPSG",9122]]], ??????? AXIS["latitude",north, ??????????? ORDER[2], ??????????? ANGLEUNIT["degree",0.0174532925199433, ??????????????? ID["EPSG",9122]]], ??? USAGE[ ??????? SCOPE["unknown"], ??????? AREA["World"], ??????? BBOX[-90,-180,90,180]]] @Edzer - should we raise an sf issue, as st_crs<- appears not to respect st_axis_order()? I'm not sure that this clarifies enough, but it's a start ... Roger
longitude <- -55.62277 latitude? <- -11.75932
st_point <- sf::st_point(c(longitude, latitude)) ll_sfc?? <- sf::st_sfc(st_point, crs = "EPSG:4326")
ll_sfc[[1]]
POINT (-55.55527 -11.51782)
ll_sfc2 <- sf::st_sfc(st_point, crs = "+proj=longlat +datum=WGS84 +no_defs?)
ll_sfc2[[1]]
POINT (-55.62277 -11.75932) ==================================================
sessionInfo()
R version 4.0.1 (2020-06-06) Platform: x86_64-apple-darwin17.0 (64-bit) Running under: macOS Catalina 10.15.5 attached base packages: [1] stats???? graphics? grDevices utils???? datasets? methods?? base other attached packages: [1] sf_0.9-5????? sits_0.9.5.1? raster_3.3-13 rgdal_1.5-12? sp_1.4-2 ===================================================================== Why is using the EPSG code resulting in a different behaviour that using the old PROJ4 text? Many thanks Gilberto =========================== Prof Dr Gilberto Camara Secretariat Director GEO - Group on Earth Observations 7 bis, Avenue de La Paix CH-1211 Geneva - Switzerland Tel: +41227308480 www.earthobservations.org ????[[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
_______________________________________________ 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, 48149 Muenster, Germany Phone: +49 251 8333081 -------------- next part -------------- A non-text attachment was scrubbed... Name: pEpkey.asc Type: application/pgp-keys Size: 3110 bytes Desc: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20200729/5bb3203e/attachment.bin>