Changes to CRS in R

Written on January 21, 2021

Keywords: R, spatial data, coordinates

There is nothing more deceptive than an obvious fact.
- Arthur Conan Doyle, The Boscombe Valley Mystery

If you haven’t heard already, big changes are afoot in the R-spatial community.

…if you were/are like me you experienced a mix of emotions. But not to worry there are loads of resources and a lot of really smart people working the issues right now.

…so expect lots of blog posts and resources.

The cliff notes version (short, short version) is that changes in how the representation of coordinate reference systems (CRS) have finally caught up with how spatial data is handled in R packages (or maybe its the otherway around). In a vignette titled title “Why have CRS, projections and transformations”, Roger Bivand explains the nitty gritty.

Here are some more resources:

YouTube lecture by Roger Bivand (link)
Associated material (link)
Bivand, R.S. Progress in the R ecosystem for representing and handling spatial data. J Geogr Syst (2020). https://doi.org/10.1007/s10109-020-00336-0

Roger also penned this post explaining the migration specific for the rgdal, sp and raster packages specific to read, write, project, and transform objects using PROJ strings (“Migration to PROJ6/GDAL3”). It gets rather complex but a good resource.

In another resource I came across in my sleuthing and troubleshooting by Edzer Pebesma and Roger Bivand discussing how GDAL and PROJ (formerly proj.4) relates to geospatial tools including several R packages in a post titled “R spatial follows GDAL and PROJ developement”. As an example, they outline the dependency for the sf package, pictured here:

Also something worth reiterating here, briefly:

PROJ provides methods for coordinate representation, conversion (projection) and transformation, and
GDAL allows reading and writing of spatial raster and vector data in a standardized form, and provides a high-level interface to PROJ for these data structures, including the representation of coordinate reference systems (CRS)

We are ultimately dealing with coordinate reference systems (or CRS) but it also goes by another name…spatial reference system (SRS). This might make more sense soon. As summarized by INBO, CRS are defined by several elements:

a coordinate system,
a ‘datum’; it localizes the geodetic coordinate system relative to the Earth and needs a geometric definition of the ellipsoid, and
only for projected CRSes coordinate conversion parameters that determine the conversion from the geodetic to the projected coordinates.

INBO did a fantastic tutorial (https://inbo.github.io/tutorials/tutorials/spatial_crs_coding/) briefly discussing on the changes and walking through the how-to for sp, sf and raster packages. The rgdal package leans heavily on the sp package…incase you were worried.

Here are some examples and things that I have learned dealing with this issue. Nothing special and I suggest visiting the resources identified above (especially https://inbo.github.io/tutorials/tutorials/spatial_crs_coding/). I am partial to the sp and rgdal packages, this is what I initially learned and got comfortable using. So lets load rgdal.

library(rgdal)

In the “good’ol days” you could define a CRS with this

utm17 <- CRS("+proj=utm +zone=17 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")

Do this now and you get…

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition

Fast-forward to now. There might be several ways to do this but the easiest I found is

utm17 <- CRS(SRS_string="EPSG:4326")

Notice the the argument SRS_string … as in spatial reference system! (I just picked that up writing this post).

Another thing in the update is the use of WKT (well-known text) over that of PROJ strings. WKT strings are interesting and provides lots of good information on the CRS (or SRS) if your into that kind of thing. To make a WKT you use the wkt() function.

utm17 <- CRS(SRS_string="EPSG:4326")

utm17.wkt=wkt(utm17)
utm17.wkt

## [1] "GEOGCRS[\"WGS 84 (with axis order normalized for visualization)\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ID[\"EPSG\",6326]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433],\n        ID[\"EPSG\",8901]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]]"

or you can print the WKT to be more readable/organized with:

cat(utm17.wkt)

## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
##     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["geodetic longitude (Lon)",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]

Further down the road when you are doing analyses or even plotting in some packages (i.e. tmap) you might get a bunch of warnings like:

Warning message:
In sp::proj4string(obj) : CRS object has comment, which is lost in output

This shouldn’t stop any of the operations but you can “mute” the warnings by running options("rgdal_show_exportToProj4_warnings"="none") in your console. I keep my “un-muted” to make sure I don’t inadvertently miss something.

If your wanting to transform a dataset from one datum to another you will need to use the WKT string. For instance I use several different state agency spatial datasets, one of which uses NAD83 HARN (which is a discarded datum…still learning about what this means) and I usually work in UTM. I find UTM CRSes easier to work with in general. Going back to the example dataset…if I read the file into R I get:

dat<-readOGR(shapefile) #just as an example

Warning message:
In OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS = dumpSRS,  :
  Discarded datum NAD83_High_Accuracy_Reference_Network in CRS definition: +proj=tmerc +lat_0=24.3333333333333 +lon_0=-81 +k=0.999941177 +x_0=200000.0001016 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs

That was enough to make my head spin…but if you notice its just a warning message and it will still read the file into the R environment. Now to transform the CRS:

dat.tran<-spTransform(dat,utm17.wkt)

But lets say you are making a SpatialPointsDataFrame, one of the arguments is proj4string (which we are moving away from and the motivation for this whole post!).

Here is some data…

dat2<-data.frame(SITE=c(1,2,3),
                 UTMX=c(590382,583910,585419),
                 UTMY=c(2830587,2821685,2819900))
dat2

##   SITE   UTMX    UTMY
## 1    1 590382 2830587
## 2    2 583910 2821685
## 3    3 585419 2819900

dat2.shp<-SpatialPointsDataFrame(dat2[,c("UTMX","UTMY")],
                                data=dat2,
                                proj4string=utm17)

This is as much as I have been able to work through these changes. It’s not huge scale changes to existing work-flows but enough to cause some heartburn.

Hope this was helpful (sorry for all the Sherlock gifs)…keep coding friends.

Share this: