Question on adding wkt<- to set CRS of a Spatial* object (get more consistent R code)
Opened this issue · 5 comments
This is about avoiding the use of proj4strings in R code, in the light of changes in PROJ/GDAL and spatial R packages. I understand that WKT2 strings are preferred and proj4strings to be avoided.
We 'had' (and of course still have) the proj4string() / proj4string<- approach of getting and setting a CRS by using the (by PROJ>=6 less/un-) supported PROJ.4 string. They work with WKT2 strings as well [not really, see next message].
I think it would still be nicer if the current wkt() function to get an object's CRS as WKT2, could also be extended to set the CRS of a Spatial object, by using wkt<- (currently not possible). I believe that would result in more consistent code. Now one still has to 'mix' proj4string<- and wkt() in R code (for sp). Or maybe, I am missing something 😉 .
library(sp)
data(meuse)
meuse.xy = meuse[c("x", "y")]
coordinates(meuse.xy) <- ~x+y
meuse.xy2 <- meuse.xy
crs_rdh <- CRS(SRS_string = "EPSG:28992")
#> Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
#> datum Amersfoort in CRS definition
# setting CRS of a Spatial object
proj4string(meuse.xy) <- crs_rdh
#
# returning CRS of an object
proj4string(meuse.xy) # this works, ...
#> Warning in proj4string(meuse.xy): CRS object has comment, which is lost in
#> output
#> [1] "+proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"
#
# ... although we now prefer:
wkt(meuse.xy)
#> [1] "PROJCRS[\"Amersfoort / RD New\",\n BASEGEOGCRS[\"Amersfoort\",\n DATUM[\"Amersfoort\",\n ELLIPSOID[\"Bessel 1841\",6377397.155,299.1528128,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4289]],\n CONVERSION[\"RD New\",\n METHOD[\"Oblique Stereographic\",\n ID[\"EPSG\",9809]],\n PARAMETER[\"Latitude of natural origin\",52.1561605555556,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",5.38763888888889,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9999079,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",155000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",463000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]],\n ID[\"EPSG\",19914]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]],\n USAGE[\n SCOPE[\"unknown\"],\n AREA[\"Netherlands - onshore\"],\n BBOX[50.75,3.2,53.7,7.22]]]"
#
# so the following would be more consistent code, but it errors:
wkt(meuse.xy2) <- crs_rdh
#> Error in wkt(meuse.xy2) <- crs_rdh: could not find function "wkt<-"
wkt(meuse.xy2)
#> Warning in wkt(meuse.xy2): CRS object has no comment
#> NULLCreated on 2020-09-16 by the reprex package (v0.3.0)
Session info
devtools::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#> setting value
#> version R version 4.0.2 (2020-06-22)
#> os Linux Mint 20
#> system x86_64, linux-gnu
#> ui X11
#> language nl_BE:nl
#> collate nl_BE.UTF-8
#> ctype nl_BE.UTF-8
#> tz Europe/Brussels
#> date 2020-09-16
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#> package * version date lib source
#> assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.0.2)
#> backports 1.1.10 2020-09-15 [1] CRAN (R 4.0.2)
#> callr 3.4.4 2020-09-07 [1] CRAN (R 4.0.2)
#> cli 2.0.2 2020-02-28 [1] CRAN (R 4.0.2)
#> crayon 1.3.4 2017-09-16 [1] CRAN (R 4.0.2)
#> desc 1.2.0 2018-05-01 [1] CRAN (R 4.0.2)
#> devtools 2.3.1 2020-07-21 [1] CRAN (R 4.0.2)
#> digest 0.6.25 2020-02-23 [1] CRAN (R 4.0.2)
#> ellipsis 0.3.1 2020-05-15 [1] CRAN (R 4.0.2)
#> evaluate 0.14 2019-05-28 [1] CRAN (R 4.0.2)
#> fansi 0.4.1 2020-01-08 [1] CRAN (R 4.0.2)
#> fs 1.5.0 2020-07-31 [1] CRAN (R 4.0.2)
#> glue 1.4.2 2020-08-27 [1] CRAN (R 4.0.2)
#> highr 0.8 2019-03-20 [1] CRAN (R 4.0.2)
#> htmltools 0.5.0 2020-06-16 [1] CRAN (R 4.0.2)
#> knitr 1.29 2020-06-23 [1] CRAN (R 4.0.2)
#> lattice 0.20-41 2020-04-02 [4] CRAN (R 4.0.0)
#> magrittr 1.5 2014-11-22 [1] CRAN (R 4.0.2)
#> memoise 1.1.0 2017-04-21 [1] CRAN (R 4.0.2)
#> pkgbuild 1.1.0 2020-07-13 [1] CRAN (R 4.0.2)
#> pkgload 1.1.0 2020-05-29 [1] CRAN (R 4.0.2)
#> prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.0.2)
#> processx 3.4.4 2020-09-03 [1] CRAN (R 4.0.2)
#> ps 1.3.4 2020-08-11 [1] CRAN (R 4.0.2)
#> R6 2.4.1 2019-11-12 [1] CRAN (R 4.0.2)
#> remotes 2.2.0 2020-07-21 [1] CRAN (R 4.0.2)
#> reprex 0.3.0 2019-05-16 [1] CRAN (R 4.0.2)
#> rgdal 1.5-16 2020-08-07 [1] CRAN (R 4.0.2)
#> rlang 0.4.7 2020-07-09 [1] CRAN (R 4.0.2)
#> rmarkdown 2.3 2020-06-18 [1] CRAN (R 4.0.2)
#> rprojroot 1.3-2 2018-01-03 [1] CRAN (R 4.0.2)
#> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 4.0.2)
#> sp * 1.4-2 2020-05-20 [1] CRAN (R 4.0.2)
#> stringi 1.5.3 2020-09-09 [1] CRAN (R 4.0.2)
#> stringr 1.4.0 2019-02-10 [1] CRAN (R 4.0.2)
#> testthat 2.3.2 2020-03-02 [1] CRAN (R 4.0.2)
#> usethis 1.6.1 2020-04-29 [1] CRAN (R 4.0.2)
#> withr 2.2.0 2020-04-20 [1] CRAN (R 4.0.2)
#> xfun 0.17 2020-09-09 [1] CRAN (R 4.0.2)
#> yaml 2.2.1 2020-02-01 [1] CRAN (R 4.0.2)
#>
#> [1] /home/floris/lib/R/library
#> [2] /usr/local/lib/R/site-library
#> [3] /usr/lib/R/site-library
#> [4] /usr/lib/R/libraryUpdate: here I am wrong in my words:
proj4string() / proj4string<- [...]
They work with WKT2 strings as well
> proj4string(meuse.xy) <- wkt(crs_rdh)
Error in CRS(value) :
PROJ4 argument-value pairs must begin with +: PROJCRS["Amersfoort / RD New",
BASEGEOGCRS["Amersfoort",
DATUM["Amersfoort",
ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4289]],
CONVERSION["RD New",
METHOD["Oblique Stereographic",
ID["EPSG",9809]],
PARAMETER["Latitude of natural origin",52.1561605555556,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",5.38763888888889,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9999079,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",155000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["Fa I think the error demonstrates that proj4string<- still requires a proj4string input (at least for the sp CRAN release), wrapped in a CRS object, even though the WKT2 is provided as a comment of the CRS object.
> comment(crs_rdh)
[1] "PROJCRS[\"Amersfoort / RD New\",\n BASEGEOGCRS[\"Amersfoort\",\n DATUM[\"Amersfoort\",\n ELLIPSOID[\"Bessel 1841\",6377397.155,299.1528128,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4289]],\n CONVERSION[\"RD New\",\n METHOD[\"Oblique Stereographic\",\n ID[\"EPSG\",9809]],\n PARAMETER[\"Latitude of natural origin\",52.1561605555556,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",5.38763888888889,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9999079,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",155000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",463000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]],\n ID[\"EPSG\",19914]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]],\n USAGE[\n SCOPE[\"unknown\"],\n AREA[\"Netherlands - onshore\"],\n BBOX[50.75,3.2,53.7,7.22]]]"And I understand that wkt() will extract the comment on the proj4string slot. I.e. these are equal:
wkt(meuse.xy)
comment(meuse.xy@proj4string)Then, is a wkt<- method (not going by proj4strings) feasible or perhaps planned?
I think the following could be a more general solution than adding wkt<-, while also making code easier to understand. It could work with the CRS objects the way these are now implemented:
- add
set_CRS(): an alias forproj4string<-() - adding a
CRSargument in Spatial* functions, to be used as a synonym of theproj4stringargument (both supported, only one allowed)
Essentially, 'proj4string' is masked here. The above could be promoted
After all, my main point was to avoid referencing the term 'proj4string' in R code, not necessarily to have 'wkt' equivalents in function and argument names. And the input for proj4string<- and for proj4string arguments is the CRS object, not the string, which this suggestion would also make less confusing.
Only and always use slot(x, "proj4string") <- CRS(...). Assigning through proj4string<- was always fragile. Create the object first and assign that object to the slot. Should we formally deprecate the proj4string<- assign to prevent misunderstandings? Unfortunately, the slot has to be called "proj4string", because S4 class object definitions should not be changed.
Thank you for explaining. I'm trying to look at it from a user's perspective.
Create the object first and assign that object to the slot.
Still I think slot creation (or checking for its existence) and subsequent CRS assignment would best be taken care of by a simple replacement function, which does not require the user to know about the proj4string slot. I believe that slot(x, "proj4string") <- is less intuitive than, for example, doing something like crs() <- (as in raster). It's just a thought - of course your advice works.
I'll continue to look at this, but in general users now should instantiate "CRS" by reading data, not by themselves.