How should I handle CRS properly after the major change in PROJ library?

Question

As explained here, the PROJ library has changed a lot and one should stop using proj.4 strings. To me it is not clear, what I should do instead. This answer explains, how I could deal with the problem without any changes in one concrete use case, but to me it seems, this is more a hack than a good idea on a long term. (At the moment, this is the only thing, I can do. The only alternative is not to update important packages.)
Is there any advice for the following process:
For example, I download data (e.g. boundaries of cantons) as geojson-files from map.geo.admin.ch.
As we don't want to do the process each time, we import the data and save it as RData-file:
boundaries <- "my_boundaries.geojson"
data <- geojsonio::geojson_read(boundaries, what = "sp")
# Warnmeldung:
# In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
# Discarded datum Unknown based on Bessel 1841 ellipsoid in CRS definition,
# but +towgs84= values preserved
data2 <- sp::spTransform(
  data, 
  sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
save(data2, "boundaries.RData")

When I analyse data according to polygons, I do the following:
boundaries <- load("boundaries.RData")
data_coords <- data.frame( # from some source mydata
  lat = mydata$lat,
  lng = mydata$lng, stringsAsFactors = FALSE)

spdf <- sp::SpatialPointsDataFrame(
  coords = data_coords,
  data = data_coords,
  proj4string = raster::crs(boundaries))
df_in <- sp::over(spdf , boundaries , fn = NULL)

So my question is:

What is the proper process for GIS-data after the major change in PROJ library?
What does this mean in our concrete case above?
I also found a nice example with plots, but it also uses proj4string. What should I do with this example?

(Unfortunately, I am not a GIS-expert)

Using CRS throws an error, if used as suggested in the answer by Roger Bivand below:
data2 <- sp::spTransform(
  data,
  sp::CRS(SRS_string = "OGC:CRS84"))
# Fehler in h(simpleError(msg, call)) : 
#   Fehler bei der Auswertung des Argumentes 'CRSobj' bei der Methodenauswahl
#   für Funktion 'spTransform': Fehlender Wert, wo TRUE/FALSE nötig ist

This is working well:
spdf <- sp::SpatialPointsDataFrame(
  coords = data_coords,
  data = data_coords)
slot(spdf, "proj4string") <- slot(boundaries, "proj4string")

Roger Bivand · Answer

Please see: https://cran.r-project.org/web/packages/rgdal/vignettes/CRS_projections_transformations.html rather than the exploratory PROJ6_GDAL3 vignette, which described processes of seeing what seemed to work.

One recent note, based on https://lists.osgeo.org/pipermail/proj/2020-September/009802.html, is that we should be extremely careful with EPSG:4326 when exchanging data with others. PROJ >= 6 and GDAL >= 3 now respect the choice of the authorities that EPSG:4326 has axis order Latitude/Longitude, unless we try to swap axes. sp/rgdal try to convert to GIS/visualization order, but may get things wrong - using "OGC:CRS84" instead of "EPSG:4326" gives us an authorised code that is Longitude first followed by Latitude by definition.

The underlying shift is from Proj4 strings including instructions for transformation from the implicit datum to WGS84 for each CRS (PROJ < 6/GDAL < 3), to a bare specification of the CRS, including its explicit datum if an authority is known.

Transformation and projection then takes place not by going from source CRS to WGS84 GEOGCRS, and on to target CRS, but searching the PROJ database for the best candidate operation for getting from source CRS to target CRS directly, reducing cumulated transformation error.

The initial question - worry about warnings - well, the warnings are there to grab your attention, and have clearly worked. We (sf/sp/rgdal) simply cannot know what is in scripts that users run. We do run extensive reverse dependency checks on CRAN packages using sf and/or sp classes, and have dramatically reduced the number of breakages since late last year (from about 100 of about 900 packages to a handful), so we hope that user scripts should not be affected by replacing Proj4 strings with WKT2_2019 strings in sf and sp (and raster) objects, but we can't know.

In particular, sp::proj4string() only reports the now deprecated Proj4 string, and many of the warnings come from using it. The quoted warning comes from the absence of a +datum= key-value pair in the output Proj4 string, although it is present in the WKT2_2019 string, and was present in the CRS object read by GDAL from the geojson file. Please trust +towgs84= much less as well, it should not be relied on either, and should be treated as deprecated.

So legacy usages like proj4string(x) <- proj4string(y) now do not do what they did before PROJ >= 6/GDAL >= 3. Use rather slot(x, "proj4string") <- slot(y, "proj4string"). If possible, avoid instantiating "CRS" objects with Proj4 strings, rather use CRS(SRS_string="OGC:CRS84") for example, instead of CRS("+proj=longlat +datum=WGS84").

How should I handle CRS properly after the major change in PROJ library?

One Answer

Add your own answers!

Ask a Question