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@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@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
pEpkey.asc
Description: application/pgp-keys
_______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo