|
| 1 | +#### GISC 422 T1 2021 |
| 2 | + |
| 3 | +# Two digressions on trend surfaces and kriging |
| 4 | + |
| 5 | +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. |
| 6 | + |
| 7 | +```{r message = FALSE} |
| 8 | +library(sf) |
| 9 | +library(tmap) |
| 10 | +library(raster) |
| 11 | +library(dplyr) |
| 12 | +library(gstat) |
| 13 | +
|
| 14 | +## The warnings are a bit out of hand on this page so |
| 15 | +options("rgdal_show_exportToProj4_warnings"="none") |
| 16 | +
|
| 17 | +volcano <- raster("data/maungawhau.tif") |
| 18 | +names(volcano) <- "data/height" |
| 19 | +
|
| 20 | +controls <- st_read("data/controls.gpkg") |
| 21 | +sites_sf <- st_read("data/sites-sf.gpkg") |
| 22 | +sites_raster <- raster("data/sites-raster.tif") |
| 23 | +``` |
| 24 | + |
| 25 | +## Some thoughts on kriging in `gstat` |
| 26 | + |
| 27 | +Although `gstat` is the 'go to' package for geostatistics in *R* and although it is very flexible, it has some limitations. Among these are: |
| 28 | + |
| 29 | +- Awkward interfaces to spatial data (not exclusive to `gstat`!) |
| 30 | +- Making the variogram by hand is tricky |
| 31 | +- Inclusion of a trend surface to run *universal kriging* is particularly challenging to get right |
| 32 | +- And, even if you can get *universal kriging* to work, you can't use a localised trend surface (admittedly this is not supported by many platforms) |
| 33 | + |
| 34 | +## Digression 1: evenly space control points |
| 35 | + |
| 36 | +Some of the challenges encountered in the main instructions are mitigated with better control points. You can use `spatstat` spatial processes to control the `st_sample` function so the code chunk below shows how this can improve kriging results. |
| 37 | + |
| 38 | +```{r} |
| 39 | +# new control set using the rSSI point process from spatstat |
| 40 | +controls_ssi <- st_read("data/interp-ext.gpkg") %>% |
| 41 | + st_sample(size = 250, type = "SSI", r = 30, n = 250) %>% |
| 42 | + st_sf() %>% |
| 43 | + st_set_crs(st_crs(controls)) %>% |
| 44 | + mutate(height = raster::extract(volcano, .)) |
| 45 | +``` |
| 46 | + |
| 47 | +Now put them on a web map and note how much more evenly spaced the `controls_ssi` points are. |
| 48 | +```{r} |
| 49 | +tmap_mode('view') |
| 50 | +tm_shape(controls) + tm_dots(col = "black") + tm_shape(controls_ssi) + tm_dots(col = "red") |
| 51 | +``` |
| 52 | + |
| 53 | +```{r} |
| 54 | +# make a new variogram |
| 55 | +v_ssi <- variogram( |
| 56 | + height ~ 1, |
| 57 | + data = as(controls, "Spatial"), |
| 58 | + cutoff = 500, |
| 59 | + width = 25, |
| 60 | +) |
| 61 | +# fit the variogram |
| 62 | +v.fit_ssi <- fit.variogram(v_ssi, vgm(model = "Gau")) |
| 63 | +# make the kriging model |
| 64 | +fit_K_ssi <- gstat( |
| 65 | + formula = height ~ 1, |
| 66 | + data = as(controls_ssi, "Spatial"), |
| 67 | + model = v.fit_ssi, |
| 68 | + nmax = 8 |
| 69 | +) |
| 70 | +# interpolate! |
| 71 | +interp_pts_K_ssi <- predict(fit_K_ssi, sites_sf) |
| 72 | +interp_K_ssi <- rasterize(as(interp_pts_K_ssi, "Spatial"), sites_raster) |
| 73 | +
|
| 74 | +persp(interp_K_ssi$var1.pred, scale = FALSE, expand = 2, theta = 35, phi = 30, lwd = 0.5) |
| 75 | +``` |
| 76 | + |
| 77 | +The lesson here is that we *always* could use better data! |
| 78 | + |
| 79 | +## Digression 2: Using a local trend surface in universal kriging |
| 80 | + |
| 81 | +In theory universal kriging models overall trends with a trend surface of some kind, then interpolates the residuals from that surface using a variogram and kriging. |
| 82 | + |
| 83 | +But we saw before that with a localised trend surface you can already get a pretty nice interpolation using that alone. However, `gstat` doesn't let you use localised trend surfaces in kriging---at least not in any simple way. In the code chunk below, I show how this limitation can potentially be sidestepped. |
| 84 | + |
| 85 | +```{r} |
| 86 | +# make a trend surface interpolation |
| 87 | +fit_TS <- gstat( |
| 88 | + formula = height ~ 1, |
| 89 | + data = as(controls, "Spatial"), |
| 90 | + nmax = 24,, |
| 91 | + degree = 3, |
| 92 | +) |
| 93 | +interp_pts_TS <- predict(fit_TS, sites_sf) |
| 94 | +interp_TS <- rasterize(as(interp_pts_TS, "Spatial"), sites_raster)$var1.pred |
| 95 | +
|
| 96 | +persp(interp_TS, scale = FALSE, expand = 2, theta = 35, phi = 30, lwd = 0.5) |
| 97 | +``` |
| 98 | + |
| 99 | +Now we proceed to krige on the residuals from this surface. |
| 100 | +```{r} |
| 101 | +# get the ts values and include in controls also making a residual |
| 102 | +controls_resid <- controls %>% |
| 103 | + mutate(ts = extract(interp_TS, .), |
| 104 | + resid = height - ts) |
| 105 | +
|
| 106 | +# now proceeed with kriging on the residual values |
| 107 | +v_resid <- variogram( |
| 108 | + resid ~ 1, |
| 109 | + data = as(controls_resid, "Spatial"), |
| 110 | + cutoff = 500, |
| 111 | + width = 25, |
| 112 | +) |
| 113 | +# fit the variogram |
| 114 | +v.fit_resid <- fit.variogram(v_resid, vgm(model = "Sph")) |
| 115 | +# make the kriging model |
| 116 | +fit_K_resid <- gstat( |
| 117 | + formula = resid ~ 1, |
| 118 | + data = as(controls_resid, "Spatial"), |
| 119 | + model = v.fit_resid, |
| 120 | + nmax = 8 |
| 121 | +) |
| 122 | +# interpolate! |
| 123 | +interp_pts_K_resid <- predict(fit_K_resid, sites_sf) |
| 124 | +interp_K_resid <- rasterize(as(interp_pts_K_resid, "Spatial"), sites_raster)$var1.pred |
| 125 | +interp_K_final <- interp_TS + interp_K_resid |
| 126 | +
|
| 127 | +persp(interp_K_final, scale = FALSE, expand = 2, theta = 35, phi = 30, lwd = 0.5) |
| 128 | +``` |
| 129 | +Back to [trend surfaces and kriging](05-trend-surfaces-and-kriging.md) | On to [splines](06-splines.md) |
0 commit comments