Merge pull request #197 from inbo/update_crs

Updating CRS specifications in existing tutorials

Author: florisvdh

content/tutorials/r_inla/spatial.pdf

@@ -16,6 +16,7 @@ output:
 library(knitr)
 opts_chunk$set(
   message = FALSE,
+  warning = FALSE, # suppress rgdal CRS warnings here and there
   dev = "cairo_pdf",
   dpi = 300
 )
@@ -329,8 +330,8 @@ PRborder %>%
   list() %>%
   Polygons(ID = 1) %>%
   list() %>%
-  SpatialPolygons(proj4string = CRS("+init=epsg:4326")) %>%
-  spTransform(CRS("+init=epsg:5880")) -> boundary
+  SpatialPolygons(proj4string = CRS(SRS_string = "EPSG:4326")) %>%
+  spTransform(CRS(SRS_string = "EPSG:5880")) -> boundary
 ```
 
 ```{r mesh-rainfall-1}
@@ -577,7 +578,7 @@ si <- inla.stack.index(stack_all, "grid")$data
 grid_data %>%
   bind_cols(model_grid$summary.fitted.values[si, ]) %>%
   `coordinates<-`(~X + Y) %>%
-  `proj4string<-`(CRS("+init=epsg:5880")) -> gd
+  `proj4string<-`(CRS(SRS_string = "EPSG:5880")) -> gd
 gd[!is.na(over(gd, boundary)), ] %>%
   as.data.frame() %>%
   ggplot() + geom_tile(aes(x = X, y = Y, fill = mean)) + coord_fixed()

content/tutorials/r_inla/spatial/index.Rmd

@@ -8,9 +8,12 @@ tags: ["gis", "r", "maps"]
 output: 
     md_document:
         preserve_yaml: true
-        variant: markdown_github
+        variant: gfm
 ---
 
+_General note: migration to the more actively developed `sf` package is currently advised by the `sp` maintainer._
+_The `sp` package, used in this tutorial, is still maintained in order to support the newest versions of the GDAL and PROJ backends._
+
 ```{r setup, include=FALSE}
 knitr::opts_chunk$set(echo = TRUE)
 ```
@@ -46,7 +49,7 @@ plot(data$lon, data$lat)
 We need to convert the `data.frame` to a `SpatialPointsDataFrame`:
 
 ```{r}
-crs_wgs84 <- CRS("+init=epsg:4326")
+crs_wgs84 <- CRS(SRS_string = "EPSG:4326")
 pts <- SpatialPointsDataFrame(data[c("lon","lat")], 
                               data[!(names(data) %in% c("lon","lat"))],
                               proj4string = crs_wgs84)
@@ -69,4 +72,4 @@ leaf_map
 
 **Nice, no?!**
 
-More information is provided [at the leaflet information website](https://rstudio.github.io/leaflet/)!
+More information is provided [at the leaflet information website](https://rstudio.github.io/leaflet/)!

content/tutorials/spatial_create_leaflet_map/index.Rmd

@@ -8,11 +8,15 @@ tags: ["gis", "r", "maps"]
 output: 
     md_document:
         preserve_yaml: true
-        variant: markdown_github
+        variant: gfm
 ---
 
-Introduction
-------------
+*General note: migration to the more actively developed `sf` package is
+currently advised by the `sp` maintainer.* *The `sp` package, used in
+this tutorial, is still maintained in order to support the newest
+versions of the GDAL and PROJ backends.*
+
+## Introduction
 
 The required packages are **leaflet** and **sp**.
 
@@ -21,10 +25,10 @@ library(leaflet)
 library(sp)
 ```
 
-Dummy data
-----------
+## Dummy data
 
-Let's create a dumy `data.frame` to play around, i.e. the three locations of INBO:
+Let’s create a dumy `data.frame` to play around, i.e. the three
+locations of INBO:
 
 ``` r
 names <- c("VAC HT","Geraardsbergen","Linkebeek")
@@ -39,15 +43,14 @@ We created three points:
 plot(data$lon, data$lat)
 ```
 
-![](index_files/figure-markdown_github/unnamed-chunk-3-1.png)
+![](index_files/figure-gfm/unnamed-chunk-3-1.png)<!-- -->
 
-Creating a map
---------------
+## Creating a map
 
 We need to convert the `data.frame` to a `SpatialPointsDataFrame`:
 
 ``` r
-crs_wgs84 <- CRS("+init=epsg:4326")
+crs_wgs84 <- CRS(SRS_string = "EPSG:4326")
 pts <- SpatialPointsDataFrame(data[c("lon","lat")], 
                               data[!(names(data) %in% c("lon","lat"))],
                               proj4string = crs_wgs84)
@@ -66,8 +69,9 @@ leaf_map <- leaflet(pts) %>%
 leaf_map
 ```
 
-![](index_files/figure-markdown_github/unnamed-chunk-6-1.png)
+![](index_files/figure-gfm/unnamed-chunk-6-1.png)<!-- -->
 
-**Nice, no?!**
+**Nice, no?\!**
 
-More information is provided [at the leaflet information website](https://rstudio.github.io/leaflet/)!
+More information is provided [at the leaflet information
+website](https://rstudio.github.io/leaflet/)\!

content/tutorials/spatial_create_leaflet_map/index.md

@@ -8,7 +8,7 @@ tags: ["gis", "r", "maps"]
 output: 
     md_document:
         preserve_yaml: true
-        variant: markdown_github
+        variant: gfm
 ---
 
 ```{r include = FALSE}
@@ -128,12 +128,15 @@ invasive_be10grid_sf %>%
 
 ## Point in polygon with the `sp` package
 
+_General note: migration to the more actively developed `sf` package is currently advised by the `sp` maintainer._
+_The `sp` package is still maintained in order to support the newest versions of the GDAL and PROJ backends._
+
 Instead of `sf` objects (= `data.frames` or `tibbles` with a geometry list-column), the `sp` package works with `Spatial` spatial data classes (which has many derived spatial data classes for points, polygons, ...).
 
 First, we need to convert the `data.frame` with point locations to a `SpatialPointsDataFrame`. We also need to ensure that the coordinate reference system (CRS) for both the point locations and the grid is the same. The data from GBIF are in WGS84 format.
 
 ```{r}
-crs_wgs84 <- CRS("+init=epsg:4326")
+crs_wgs84 <- CRS(SRS_string = "EPSG:4326")
 coord <- invasive_species %>%
     select(decimalLongitude, decimalLatitude)
 invasive_spatial <- SpatialPointsDataFrame(coord,
@@ -148,6 +151,19 @@ be10grid <- readOGR(dsn = file.path(tempdir(), "Belgium.sqlite"),
                   layer = "be_10km")
 ```
 
+Note the warning: it is because some PROJ.4 information, i.e. the string to represent the coordinate reference system, is not supported anymore in the current geospatial GDAL and PROJ libraries (the background workhorses for spatial R packages).
+The spatialite database from the EEA website (with the 10 km x 10 km grid) still uses the older PROJ.4 string .
+Because the `rgdal` package is still backwards compatible, we should not (yet) worry about this: `rgdal` does the translation for the newer GDAL 3 and PROJ >= 6.
+Do know that, instead of _PROJ.4_ strings, the _WKT2_ string is now used in R to better represent coordinate reference systems (so it would best be incorporated in the EEA data source).
+Just compare these:
+
+```{r}
+# PROJ.4 string = old; used by PROJ 4
+proj4string(be10grid)  # or: be10grid@proj4string
+# WKT2 string = new and much better; used by PROJ >= 6
+wkt(be10grid) %>% cat  # or: be10grid@proj4string %>% comment %>% cat
+```
+
 
 We transform the 10 km x 10 km grid to the same CRS system: