1. Introduction

This notebook illustrates spatial interpolation using data obtained in a field campaign over a potato crop farm. A lot of code is borrowed from Introduction to Spatial Data Programming with R by Michael Dorman (2021). Its main purpose is to guide data processing tasks for the HERMES project: Agro-Geoinformática: Detección temprana in situ de enfermedades a nivel foliar empleando imágenes espectrales.

Let’s start by cleaning memory and loading the required libraries:

rm(list=ls())
library(sf)
library(stars)
library(leaflet)
library(gstat)
library(automap)
library(raster)
library(RColorBrewer)

2. Reading input data

Let’s read a shapefile representing the study area:

(aoi <- st_read("C:\\datos\\crops\\tenjo\\muestreo\\area_estudio.shp"))
## Reading layer `area_estudio' from data source 
##   `C:\datos\crops\tenjo\muestreo\area_estudio.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 990223.4 ymin: 1027816 xmax: 990363.4 ymax: 1027945
## CRS:           NA

Note that the CRS has not been defined yet.

Now, let’s read the sampling points:

(ptos <- st_read("../crops/tenjo/muestreo/puntos_13_03.shp"))
## Reading layer `puntos_13_03' from data source 
##   `C:\datos\crops\tenjo\muestreo\puntos_13_03.shp' using driver `ESRI Shapefile'
## Simple feature collection with 60 features and 19 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 990235 ymin: 1027825 xmax: 990357.3 ymax: 1027938
## Projected CRS: MAGNA-SIRGAS / Colombia Bogota zone

Note that the points are georeferenced in a plane coordinate reference system:

st_crs(ptos)
## Coordinate Reference System:
##   User input: MAGNA-SIRGAS / Colombia Bogota zone 
##   wkt:
## PROJCRS["MAGNA-SIRGAS / Colombia Bogota zone",
##     BASEGEOGCRS["MAGNA-SIRGAS",
##         DATUM["Marco Geocentrico Nacional de Referencia",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4686]],
##     CONVERSION["Colombia MAGNA Bogota zone",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",4.59620041666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-74.0775079166667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",1000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",1000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping (large scale)."],
##         AREA["Colombia - onshore between 1°30'W and 1°30'E of Bogota (75°35'W and 72°35'W of Greenwich)."],
##         BBOX[-2.51,-75.59,11.82,-72.58]],
##     USAGE[
##         SCOPE["Topographic mapping (small scale)."],
##         AREA["Colombia - mainland onshore."],
##         BBOX[-4.23,-79.1,12.52,-66.87]],
##     ID["EPSG",3116]]

Now, let’s check the values of soil moisture in the top 5 cm as measured by a time-domain reflectometer (TDR) device. This property is stored in the attribute T of ptos:

var = "T"
ptos[var]

Then, visualization time:

getColor <- function(ptos) {
  sapply(ptos$T, function(T) {
  if(T >= 50) {
    "blue"
  } else if(T < 10) {
    "orange"
  } else {
    "green"
  } })
}

icons <- awesomeIcons(
  icon = 'ios-close',
  iconColor = 'black',
  library = 'ion',
  markerColor = getColor(ptos)
)

leaflet(data = ptos) %>% addTiles() %>%
  addAwesomeMarkers(~W, ~N, icon=icons, popup = ~as.character(T), label = ~as.character(T))

Note that, at this level of spatial detail, there is not OpenStreetMap data. We could try another type of backdrop tiles.

We also need a raster for the study area. Let´s load a subset of a multispectral imagery collected over the field:

x = read_stars("../crops/tenjo/marzo13/output/V1_B6_lote.tif")

Let’s check what is x:

plot(x)
## downsample set to c(2,2)

This is a thermal band collected using an Altum camera from Micasense. Note the range of temperature values in Celsius.

3. Spatial interpolation

3.1 Inverse weighted distance (IDW)

To interpolate using the IDW method we create the following gstat object, specifying just the formula and data:

g = gstat(formula = T ~ 1, data = ptos)

Now that our model is defined, we can use the predict function to actually interpolate. The predict function accepts:

  • A raster—stars object, such as x
  • A model—gstat object, such as g

The raster serves for two purposes:

  • Specifying the locations where we want to make predictions (in all methods)
  • Specifying covariate values (in Universal Kriging only)

The following expression interpolates soil moisture values according to the model defined in g and the raster template defined in x:

z = predict(g, x)
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]

The resulting stars object has two attributes:

  • var1.pred—the predicted values
  • var1.var—the variance (for Kriging only)

Look at z:

z
## stars object with 2 dimensions and 2 attributes
## attribute(s), summary of first 1e+05 cells:
##               Min.  1st Qu.   Median     Mean  3rd Qu.     Max.   NA's
## var1.pred  23.3339 23.60441 23.82117 23.86856 24.11519 24.73005  99060
## var1.var        NA       NA       NA      NaN       NA       NA 100000
## dimension(s):
##   from   to  offset delta                       refsys point values x/y
## x    1 3497  990223  0.04 MAGNA-SIRGAS / Colombia B...    NA   NULL [x]
## y    1 3227 1027945 -0.04 MAGNA-SIRGAS / Colombia B...    NA   NULL [y]

We can subset just the first attribute and rename it to “soil_moisture”:

z = z["var1.pred",,]
names(z) = "soil_moisture"

The interpolated soil moisture looks like this:

b = seq(0, 70, 2)
plot(z, breaks = b, col = hcl.colors(length(b)-1, "Spectral"), reset = FALSE)
## downsample set to c(2,2)
plot(st_geometry(ptos), pch = 3, add = TRUE)
contour(z, breaks = b, add = TRUE)

We can write the estimated property to disk:

write_stars(z, "SM_IDW.tif")

3.2 Ordinary Kriging (OK)

Kriging methods require a variogram model. The variogram model is an objective way to quantify the autocorrelation pattern in the data, and assign weights accordingly when making predictions.

As a first step, we can calculate and examine the empirical variogram using the variogram function. The function requires two arguments:

  • formula: Specifies the dependent variable and the covariates, just like in gstat
  • data: The point layer with the dependent variable and covariates as point attributes

For example, the following expression calculates the empirical variogram of annual, with no covariates:

v_emp_ok = variogram(T ~ 1, ptos)

Using plot to examine it we can examine the variogram:

plot(v_emp_ok)

There are several ways to fit a variogram model to an empirical variogram. We will use the function autofitVariogram from package automap:

v_mod_ok = autofitVariogram(T ~ 1, as(ptos, "Spatial"))

The function chooses the best fitting type of model, and also fine tunes its parameters.

Note that the autofitVariogram function does not work on sf objects, which is why we convert the object to a SpatialPointsDataFrame (package sp).

The fitted model can be plotted with plot:

plot(v_mod_ok)

Note that the fitted model is a Ste model, that is, a Matern, M. Stein’s parameterization model.

We can use show.vgms() to display variogram model types:

show.vgms()

The resulting object is actually a list with several components, including the empirical variogram and the fitted variogram model. The $var_model component of the resulting object contains the actual model:

v_mod_ok$var_model

The variogram model can then be passed to the gstat function in order to carry on with the Ordinary Kriging interpolation:

g = gstat(formula = T ~ 1, model = v_mod_ok$var_model, data = ptos)
z = predict(g, x)
## [using ordinary kriging]

Again, we will subset the predicted values attribute and rename it:

z = z["var1.pred",,]
names(z) = "soil_moisture"

The Ordinary Kriging predictions are shown here:

b = seq(0, 70, 1)
plot(z, breaks = b, col = hcl.colors(length(b)-1, "Spectral"), reset = FALSE)
## downsample set to c(2,2)
plot(st_geometry(ptos), pch = 3, add = TRUE)
contour(z, breaks = b, add = TRUE)

What is z?

z
## stars object with 2 dimensions and 1 attribute
## attribute(s), summary of first 1e+05 cells:
##                    Min.  1st Qu.  Median     Mean  3rd Qu.     Max.  NA's
## soil_moisture  18.97136 19.35219 19.7157 19.79669 20.16822 21.24577 99060
## dimension(s):
##   from   to  offset delta                       refsys point values x/y
## x    1 3497  990223  0.04 MAGNA-SIRGAS / Colombia B...    NA   NULL [x]
## y    1 3227 1027945 -0.04 MAGNA-SIRGAS / Colombia B...    NA   NULL [y]

Again, write the interpolated raster:

write_stars(z, "SM_OK.tif")

That´s all for now.

4. Reproducibility

If you reuse this code for your research please cite as: Lizarazo, I., 2021. Looking again at spatial interpolation. Available at https://rpubs.com/ials2un/laaspint.

sessionInfo()
## R version 4.0.5 (2021-03-31)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.1252  LC_CTYPE=Spanish_Colombia.1252   
## [3] LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] RColorBrewer_1.1-2 raster_3.4-10      automap_1.0-14     sp_1.4-5          
## [5] gstat_2.0-7        leaflet_2.0.4.1    stars_0.5-3        abind_1.4-5       
## [9] sf_1.0-1          
## 
## loaded via a namespace (and not attached):
##  [1] zoo_1.8-9          tidyselect_1.1.1   xfun_0.22          bslib_0.2.5.1     
##  [5] purrr_0.3.4        lattice_0.20-41    vctrs_0.3.8        generics_0.1.0    
##  [9] htmltools_0.5.1.1  spacetime_1.2-5    yaml_2.2.1         utf8_1.2.1        
## [13] rlang_0.4.10       e1071_1.7-7        jquerylib_0.1.4    pillar_1.6.1      
## [17] glue_1.4.2         DBI_1.1.1          plyr_1.8.6         lifecycle_1.0.0   
## [21] stringr_1.4.0      htmlwidgets_1.5.3  codetools_0.2-18   evaluate_0.14     
## [25] knitr_1.32         crosstalk_1.1.1    parallel_4.0.5     class_7.3-18      
## [29] fansi_0.5.0        highr_0.9          xts_0.12.1         Rcpp_1.0.6        
## [33] KernSmooth_2.23-18 classInt_0.4-3     lwgeom_0.2-6       jsonlite_1.7.2    
## [37] FNN_1.1.3          digest_0.6.27      stringi_1.6.2      dplyr_1.0.7       
## [41] grid_4.0.5         rgdal_1.5-23       tools_4.0.5        magrittr_2.0.1    
## [45] sass_0.4.0         proxy_0.4-26       tibble_3.1.2       crayon_1.4.1      
## [49] pkgconfig_2.0.3    ellipsis_0.3.2     reshape_0.8.8      rmarkdown_2.7     
## [53] R6_2.5.0           units_0.7-2        intervals_0.15.2   compiler_4.0.5