updated interpolation materials for 2023

Author: David O'Sullivan

labs/interpolation/01-overview-of-the-approach.Rmd

@@ -1,4 +1,4 @@
-#### GISC 422 T1 2021
+#### GISC 422 T2 2023
 
 # The overall approach to interpolation
 The *R* ecosystem's approach to spatial interpolation seems pretty complicated at first. It's not exactly simple, although the overall concept is straightforward enough and some of the apparent complexity has more to do with making different data types interact with one another successfully.

labs/interpolation/01-overview-of-the-approach.md

@@ -1,4 +1,4 @@
-#### GISC 422 T1 2021
+#### GISC 422 T2 2023
 
 # The overall approach to interpolation
 The *R* ecosystem's approach to spatial interpolation seems pretty complicated at first. It's not exactly simple, although the overall concept is straightforward enough and some of the apparent complexity has more to do with making different data types interact with one another successfully.

labs/interpolation/02-example-dataset.Rmd

@@ -1,19 +1,19 @@
-#### GISC 422 T1 2021
+#### GISC 422 T2 2023
 # Basics of working with raster data
 First load some libraries
 
 ```{r}
 library(sf)
 library(tmap)
-library(raster)
+library(terra)
 ```
 The new(ish) to us kid on the block here is `raster` for handling gridded raster datasets. One thing to be very aware of is that `raster` masks the `select` function from `dplyr` so you have to specify `dplyr::select` when using the `select` to tidy up datasets during data preparation.
 
 ## Read a raster dataset
 We are using a simple example of the elevation of Maungawhau (Mt Eden) in Auckland to demonstrate the interpolation methods. There is a version of this dataset available as standard in *R* but I made raster version of it for us to work with. So load this with the raster package:
 
 ```{r}
-volcano <- raster("data/maungawhau.tif")
+volcano <- rast("data/maungawhau.tif")
 ```
 
 Confusingly, when you read in a raster dataset it names the associated numerical data using the filename, so we rename that to `height` which is more appropriate for our purposes.

labs/interpolation/02-example-dataset.md

@@ -1,19 +1,19 @@
-#### GISC 422 T1 2021
+#### GISC 422 T2 2023
 # Basics of working with raster data
 First load some libraries
 
 ```{r}
 library(sf)
 library(tmap)
-library(raster)
+library(terra)
 ```
 The new(ish) to us kid on the block here is `raster` for handling gridded raster datasets. One thing to be very aware of is that `raster` masks the `select` function from `dplyr` so you have to specify `dplyr::select` when using the `select` to tidy up datasets during data preparation.
 
 ## Read a raster dataset
 We are using a simple example of the elevation of Maungawhau (Mt Eden) in Auckland to demonstrate the interpolation methods. There is a version of this dataset available as standard in *R* but I made raster version of it for us to work with. So load this with the raster package:
 
 ```{r}
-volcano <- raster("data/maungawhau.tif")
+volcano <- rast("data/maungawhau.tif")
 ```
 
 Confusingly, when you read in a raster dataset it names the associated numerical data using the filename, so we rename that to `height` which is more appropriate for our purposes.

labs/interpolation/03-preparing-for-interpolation.Rmd

@@ -1,14 +1,14 @@
-#### GISC 422 T1 2021
+#### GISC 422 T2 2023
 
 # Preparing for interpolation
 Run this first to make sure all the data and packages you need are loaded:
 ```{r}
 library(sf)
 library(tmap)
-library(raster)
+library(terra)
 library(dplyr)
 
-volcano <- raster("data/maungawhau.tif")
+volcano <- rast("data/maungawhau.tif")
 names(volcano) <- "height"
 ```
 
@@ -34,14 +34,19 @@ interp_ext <- controls %>%
 ## Control points
 For the demonstration data, we already know the result.
 
-Normally, we would have a set of control points in some spatial format and would simply read them with `sf::st_read`. Here, we will make set of random control points to work with in the other steps of these instructions when we are using the Maungawhau data. We start from the interpolation extent, and use `st_sample` to get a specified random number of points in the extent, then convert it to an `sf` dataset. Finally we use the `raster::extract` function to extract values from the raster dataset and assign their values to a height attribute of the `sf` dataset.
+Normally, we would have a set of control points in some spatial format and would simply read them with `sf::st_read`. Here, we will make set of random control points to work with in the other steps of these instructions when we are using the Maungawhau data. We start from the interpolation extent, and use `st_sample` to get a specified random number of points in the extent, then convert it to an `sf` dataset. Finally we use the `terra::extract` function to extract values from the raster dataset and assign their values to a height attribute of the `sf` dataset.
 
 ```{r}
 controls <- interp_ext %>%
   st_sample(size = 250) %>%
   st_sf() %>%
-  st_set_crs(st_crs(interp_ext)) %>%
-  mutate(height = raster::extract(volcano, .))
+  st_set_crs(st_crs(interp_ext))
+
+heights <- controls %>%
+  extract(x = volcano)
+
+controls <- controls %>%
+  mutate(height = heights$height)
 ```
 
 Every time you run the above you will get a different random set of the specified number of control point locations. It is useful to map them on top of the underlying data and think about how many you might need to get a reasonable representation of the height map of Maungawhau.
@@ -66,7 +71,7 @@ st_read("data/controls.gpkg") %>%
 Remember that if you change the control points, you should also change this file to keep them matching.
 
 ## Make a set of locations to interpolate
-Unlike the previous step which may not be necessary when you are provided with control points to interpolate directly, this step is always required. Basically, *R* wants a raster layer to *interpolate into*. We'll call this `sites` and make `sf`, `raster` and simple xyz versions of these.
+Unlike the previous step which may not be necessary when you are provided with control points to interpolate directly, this step is always required. Basically, *R* wants a raster layer to *interpolate into*. We'll call this `sites` and make `sf`, `terra` and simple xyz versions of these.
 
 ```{r}
 sites_sf <- interp_ext %>% # start with the extent
@@ -75,12 +80,12 @@ sites_sf <- interp_ext %>% # start with the extent
   st_set_crs(st_crs(interp_ext))
 
 sites_xyz <- sites_sf %>%
-  cbind(st_coordinates(.)) %>%
+  bind_cols(st_coordinates(.)) %>%
   st_drop_geometry() %>%
   mutate(Z = 0)
 
 sites_raster <- sites_xyz %>%
-  rasterFromXYZ()
+  rast(type = "xyz")
 crs(sites_raster) <- st_crs(controls)$wkt
 ```
 

labs/interpolation/03-preparing-for-interpolation.md

@@ -1,14 +1,14 @@
-#### GISC 422 T1 2021
+#### GISC 422 T2 2023
 
 # Preparing for interpolation
 Run this first to make sure all the data and packages you need are loaded:
 ```{r}
 library(sf)
 library(tmap)
-library(raster)
+library(terra)
 library(dplyr)
 
-volcano <- raster("data/maungawhau.tif")
+volcano <- rast("data/maungawhau.tif")
 names(volcano) <- "height"
 ```
 
@@ -34,14 +34,19 @@ interp_ext <- controls %>%
 ## Control points
 For the demonstration data, we already know the result.
 
-Normally, we would have a set of control points in some spatial format and would simply read them with `sf::st_read`. Here, we will make set of random control points to work with in the other steps of these instructions when we are using the Maungawhau data. We start from the interpolation extent, and use `st_sample` to get a specified random number of points in the extent, then convert it to an `sf` dataset. Finally we use the `raster::extract` function to extract values from the raster dataset and assign their values to a height attribute of the `sf` dataset.
+Normally, we would have a set of control points in some spatial format and would simply read them with `sf::st_read`. Here, we will make set of random control points to work with in the other steps of these instructions when we are using the Maungawhau data. We start from the interpolation extent, and use `st_sample` to get a specified random number of points in the extent, then convert it to an `sf` dataset. Finally we use the `terra::extract` function to extract values from the raster dataset and assign their values to a height attribute of the `sf` dataset.
 
 ```{r}
 controls <- interp_ext %>%
   st_sample(size = 250) %>%
   st_sf() %>%
-  st_set_crs(st_crs(interp_ext)) %>%
-  mutate(height = raster::extract(volcano, .))
+  st_set_crs(st_crs(interp_ext))
+
+heights <- controls %>%
+  extract(x = volcano)
+
+controls <- controls %>%
+  mutate(height = heights$height)
 ```
 
 Every time you run the above you will get a different random set of the specified number of control point locations. It is useful to map them on top of the underlying data and think about how many you might need to get a reasonable representation of the height map of Maungawhau.
@@ -66,7 +71,7 @@ st_read("data/controls.gpkg") %>%
 Remember that if you change the control points, you should also change this file to keep them matching.
 
 ## Make a set of locations to interpolate
-Unlike the previous step which may not be necessary when you are provided with control points to interpolate directly, this step is always required. Basically, *R* wants a raster layer to *interpolate into*. We'll call this `sites` and make `sf`, `raster` and simple xyz versions of these.
+Unlike the previous step which may not be necessary when you are provided with control points to interpolate directly, this step is always required. Basically, *R* wants a raster layer to *interpolate into*. We'll call this `sites` and make `sf`, `terra` and simple xyz versions of these.
 
 ```{r}
 sites_sf <- interp_ext %>% # start with the extent
@@ -75,12 +80,12 @@ sites_sf <- interp_ext %>% # start with the extent
   st_set_crs(st_crs(interp_ext))
 
 sites_xyz <- sites_sf %>%
-  cbind(st_coordinates(.)) %>%
+  bind_cols(st_coordinates(.)) %>%
   st_drop_geometry() %>%
   mutate(Z = 0)
 
 sites_raster <- sites_xyz %>%
-  rasterFromXYZ()
+  rast(type = "xyz")
 crs(sites_raster) <- st_crs(controls)$wkt
 ```
 

labs/interpolation/04-nn-and-idw.Rmd

@@ -1,20 +1,20 @@
-#### GISC 422 T1 2021
+#### GISC 422 T2 2023
 # Near neighbour and inverse-distance weighted interpolation
 Run this first to make sure all the data and packages you need are loaded. If any data are missing you probably didn't make them in one of the previous instruction pages.
 
 ```{r}
 library(sf)
 library(tmap)
-library(raster)
+library(terra)
 library(dplyr)
 library(gstat)
 
-volcano <- raster("data/maungawhau.tif")
+volcano <- rast("data/maungawhau.tif")
 names(volcano) <- "data/height"
 
 controls <- st_read("data/controls.gpkg")
 sites_sf <- st_read("data/sites-sf.gpkg")
-sites_raster <- raster("data/sites-raster.tif")
+sites_raster <- rast("data/sites-raster.tif")
 ```
 
 ## Inverse-distance weighted interpolation
@@ -23,9 +23,9 @@ These two methods are very similar, and IDW is actually *more general* so we'll
 As with all the `gstat` methods we use the `gstat::gstat` function to make a statistical model, and then apply it using the `raster::interpolate` function.
 
 ```{r}
-fit_IDW <- gstat(                 # makes a model 
-  formula = height ~ 1,           # The column `height` is what we are interested in
-  data = as(controls, "Spatial"), # using sf but converting to sp, which is required
+fit_IDW <- gstat(                   # makes a model 
+  formula = height ~ 1,             # The column `height` is what we are interested in
+  data = as(controls, "Spatial"),   # using sf but converting to sp, which is required
   set = list(idp = 2),
   # nmax = 12, maxdist = 100        # you can experiment with these options later...
 )
@@ -37,7 +37,7 @@ Having made the model (called `fit_IDW`) we pass it to the `predict` function to
 
 ```{r}
 interp_pts_IDW <- predict(fit_IDW, sites_sf)
-interp_IDW <- rasterize(as(interp_pts_IDW, "Spatial"), sites_raster, "var1.pred")
+interp_IDW <- rasterize(as(interp_pts_IDW, "SpatVector"), sites_raster, "var1.pred")
 names(interp_IDW) <- "height" # rename the variable to something more friendly
 ```
 
@@ -70,7 +70,7 @@ fit_NN <- gstat( # makes a model
 
 # and interpolate like before
 interp_pts_NN <- predict(fit_NN, sites_sf)
-interp_NN <- rasterize(as(interp_pts_NN, "Spatial"), sites_raster, "var1.pred")
+interp_NN <- rasterize(as(interp_pts_NN, "SpatVector"), sites_raster, "var1.pred")
 names(interp_NN) <- "height"
 
 # and display it

labs/interpolation/04-nn-and-idw.md

@@ -1,20 +1,20 @@
-#### GISC 422 T1 2021
+#### GISC 422 T2 2023
 # Near neighbour and inverse-distance weighted interpolation
 Run this first to make sure all the data and packages you need are loaded. If any data are missing you probably didn't make them in one of the previous instruction pages.
 
 ```{r}
 library(sf)
 library(tmap)
-library(raster)
+library(terra)
 library(dplyr)
 library(gstat)
 
-volcano <- raster("data/maungawhau.tif")
+volcano <- rast("data/maungawhau.tif")
 names(volcano) <- "data/height"
 
 controls <- st_read("data/controls.gpkg")
 sites_sf <- st_read("data/sites-sf.gpkg")
-sites_raster <- raster("data/sites-raster.tif")
+sites_raster <- rast("data/sites-raster.tif")
 ```
 
 ## Inverse-distance weighted interpolation
@@ -23,9 +23,9 @@ These two methods are very similar, and IDW is actually *more general* so we'll
 As with all the `gstat` methods we use the `gstat::gstat` function to make a statistical model, and then apply it using the `raster::interpolate` function.
 
 ```{r}
-fit_IDW <- gstat(                 # makes a model 
-  formula = height ~ 1,           # The column `height` is what we are interested in
-  data = as(controls, "Spatial"), # using sf but converting to sp, which is required
+fit_IDW <- gstat(                   # makes a model 
+  formula = height ~ 1,             # The column `height` is what we are interested in
+  data = as(controls, "Spatial"),   # using sf but converting to sp, which is required
   set = list(idp = 2),
   # nmax = 12, maxdist = 100        # you can experiment with these options later...
 )
@@ -37,7 +37,7 @@ Having made the model (called `fit_IDW`) we pass it to the `predict` function to
 
 ```{r}
 interp_pts_IDW <- predict(fit_IDW, sites_sf)
-interp_IDW <- rasterize(as(interp_pts_IDW, "Spatial"), sites_raster, "var1.pred")
+interp_IDW <- rasterize(as(interp_pts_IDW, "SpatVector"), sites_raster, "var1.pred")
 names(interp_IDW) <- "height" # rename the variable to something more friendly
 ```
 
@@ -70,7 +70,7 @@ fit_NN <- gstat( # makes a model
 
 # and interpolate like before
 interp_pts_NN <- predict(fit_NN, sites_sf)
-interp_NN <- rasterize(as(interp_pts_NN, "Spatial"), sites_raster, "var1.pred")
+interp_NN <- rasterize(as(interp_pts_NN, "SpatVector"), sites_raster, "var1.pred")
 names(interp_NN) <- "height"
 
 # and display it

labs/interpolation/05-trend-surfaces-and-kriging.Rmd

@@ -1,4 +1,4 @@
-#### GISC 422 T1 2021
+#### GISC 422 T2 2023
 
 # Trend surfaces and kriging
 
@@ -7,16 +7,16 @@ Run this first to make sure all the data and packages you need are loaded. If an
 ```{r}
 library(sf)
 library(tmap)
-library(raster)
+library(terra)
 library(dplyr)
 library(gstat)
 
-volcano <- raster("data/maungawhau.tif")
+volcano <- rast("data/maungawhau.tif")
 names(volcano) <- "data/height"
 
 controls <- st_read("data/controls.gpkg")
 sites_sf <- st_read("data/sites-sf.gpkg")
-sites_raster <- raster("data/sites-raster.tif")
+sites_raster <- rast("data/sites-raster.tif")
 ```
 
 There are many different styles of kriging. We'll work here with universal kriging which models variation in the data with two components a *trend surface* and a *variogram* which models how control points vary with distance from one another. So... to perform kriging we have to consider each of these elements in turn.
@@ -38,7 +38,7 @@ The form of the trend surface function is specified by the `degree` parameter an
 
 ```{r}
 interp_pts_TS <- predict(fit_TS, sites_sf)
-interp_TS <- rasterize(as(interp_pts_TS, "Spatial"), sites_raster)
+interp_TS <- rasterize(as(interp_pts_TS, "SpatVector"), sites_raster, field = c("var1.pred", "var1.var"))
 
 persp(interp_TS$var1.pred, scale = FALSE, expand = 2, theta = 35, phi = 30, lwd = 0.5)
 ```
@@ -100,7 +100,7 @@ I have found that it is important to set a fairly low `nmax` with these data. I
 
 ```{r}
 interp_pts_K <- predict(fit_K, sites_sf)
-interp_K <- rasterize(as(interp_pts_K, "Spatial"), sites_raster)
+interp_K <- rasterize(as(interp_pts_K, "SpatVector"), sites_raster, field = c("var1.pred", "var1.var"))
 
 persp(interp_K$var1.pred, scale = FALSE, expand = 2, theta = 35, phi = 30, lwd = 0.5)
 ```

labs/interpolation/05-trend-surfaces-and-kriging.md

@@ -1,4 +1,4 @@
-#### GISC 422 T1 2021
+#### GISC 422 T2 2023
 
 # Trend surfaces and kriging
 
@@ -7,16 +7,16 @@ Run this first to make sure all the data and packages you need are loaded. If an
 ```{r}
 library(sf)
 library(tmap)
-library(raster)
+library(terra)
 library(dplyr)
 library(gstat)
 
-volcano <- raster("data/maungawhau.tif")
+volcano <- rast("data/maungawhau.tif")
 names(volcano) <- "data/height"
 
 controls <- st_read("data/controls.gpkg")
 sites_sf <- st_read("data/sites-sf.gpkg")
-sites_raster <- raster("data/sites-raster.tif")
+sites_raster <- rast("data/sites-raster.tif")
 ```
 
 There are many different styles of kriging. We'll work here with universal kriging which models variation in the data with two components a *trend surface* and a *variogram* which models how control points vary with distance from one another. So... to perform kriging we have to consider each of these elements in turn.
@@ -38,7 +38,7 @@ The form of the trend surface function is specified by the `degree` parameter an
 
 ```{r}
 interp_pts_TS <- predict(fit_TS, sites_sf)
-interp_TS <- rasterize(as(interp_pts_TS, "Spatial"), sites_raster)
+interp_TS <- rasterize(as(interp_pts_TS, "SpatVector"), sites_raster, field = c("var1.pred", "var1.var"))
 
 persp(interp_TS$var1.pred, scale = FALSE, expand = 2, theta = 35, phi = 30, lwd = 0.5)
 ```
@@ -100,7 +100,7 @@ I have found that it is important to set a fairly low `nmax` with these data. I
 
 ```{r}
 interp_pts_K <- predict(fit_K, sites_sf)
-interp_K <- rasterize(as(interp_pts_K, "Spatial"), sites_raster)
+interp_K <- rasterize(as(interp_pts_K, "SpatVector"), sites_raster, field = c("var1.pred", "var1.var"))
 
 persp(interp_K$var1.pred, scale = FALSE, expand = 2, theta = 35, phi = 30, lwd = 0.5)
 ```