inbo/tutorials

Check robustness of CRS specifications in tutorials

Closed this issue · 2 comments

Given current migrations of CRS representation in sp, sf>=0.9-0, rgdal and dependent packages, our tutorials will need to be checked for existing non-compliant CRS specifications. This may be needed repeatedly, as things are still evolving (both in PROJ & GDAL and in the major R packages).

Also some packages, e.g. raster (see https://github.com/rspatial/raster/issues), are currently not yet compliant, hence they must be used with care at least with PROJ >=6 (but e.g. with sf as well).

In general, changes in geospatial packages are ongoing and more changes are to be expected, at least in the function's output.

For now: avoid hard-coded proj4strings, including things like "+init=epsg:31370". More info and coding examples:

A grep on current master (9c04f82) reveals the following occurrences of obsolete CRS specifications (either in code or text):

./content/tutorials/r_inla/spatial/index.Rmd:332:  SpatialPolygons(proj4string = CRS("+init=epsg:4326")) %>%
./content/tutorials/r_inla/spatial/index.Rmd:333:  spTransform(CRS("+init=epsg:5880")) -> boundary
./content/tutorials/r_inla/spatial/index.Rmd:580:  `proj4string<-`(CRS("+init=epsg:5880")) -> gd
./content/tutorials/spatial_create_leaflet_map/index.Rmd:49:crs_wgs84 <- CRS("+init=epsg:4326")
./content/tutorials/spatial_point_in_polygon/index.Rmd:136:crs_wgs84 <- CRS("+init=epsg:4326")
./content/tutorials/spatial_standards_raster/index.Rmd:61:         crs = CRS("+init=epsg:31370")) %>% 
./content/tutorials/spatial_transform_crs/index.Rmd:39:proj4string(ngi) <- CRS("+init=epsg:3812")
./content/tutorials/spatial_transform_crs/index.Rmd:43:We could verify the correctness of this position by plotting in on a map. Here we use the `leaflet` package which requires the data to be in the "WGS84" coordinate reference system. Therefore we use `spTransform` to do this transformation. "WGS84" has EPSG number 4326. But here the coordinate reference system string itself is easier to memorise: "+proj=longlat".
./content/tutorials/spatial_transform_crs/index.Rmd:46:ngi_ll <- spTransform(ngi, CRS("+proj=longlat"))
./content/tutorials/spatial_transform_crs/index.Rmd:59:    shortened_PROJ.4_string = c("+init=epsg:4326", "+init=epsg:31370", "+init=epsg:3812", "+init=epsg:3857")

Planning to make a PR for this.

Also the long-standing #27 is affected. Planning to have a look at it so that it can be rendered + merged.