[R-sig-Geo] Help when moving from PROJ4 to PROJ6
Roger Bivand
Roger@B|v@nd @end|ng |rom nhh@no
Tue Jul 28 15:34:31 CEST 2020
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 using r-project.org
> 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: Roger.Bivand using nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
More information about the R-sig-Geo
mailing list