1. Why this notebook?

This is an update for a previous R Notebook written to illustrate how to conduct spatial interpolation in R which may be found here.

This update aims to help Geomatica Basica students at Universidad Nacional to fix a problem originated from the shift from PROJ4 to PROJ6 in the recent versions of the spatial packages of R (and gdal for that matter). In the previous notebook I used sp::proj4string() to get the projection code / CRS argument from an sp object and use that code further to create an empty grid for conducting IDW interpolation. Now I get a warning stating that the grid and the input data have different CRSs.

In order to solve such a problem, I use here the WKT format (instead of the PROJ4 format) to pass the CRS to the empty grid. In addition, I use the UTM 18N CRS to project the input data from geographic coordinates into map coordinates (instead of the MAGNA SIRGAS Colombia Bogota Zone). Otherwise, the problem persisted.

Note that this notebook does not cover all items contained in the previous one. You may need to go back to it to complete your technical report.

Again, I focus here on explaining the technical procedures rather than on interpreting or analyzing the outputs.

2. Input data

I will use a data set of monthly precipitation over Departamento de Boyacá during July 2019. It has been obtained using climateR as explained in a previous notebook.

Let’s start by cleaning the memory and reading the data:

rm(list=ls())
library(sf)
library(rgdal)
library(rasterVis)
library(dplyr)
library(raster)
library(gstat)
library(geoR)
library(berryFunctions)
(boyaca <- st_read("../../datos/shapefiles/15_BOYACA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))
## Reading layer `MGN_MPIO_POLITICO' from data source 
##   `/Users/ivanlizarazo/Documents/ivan/GB/datos/shapefiles/15_BOYACA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 123 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## Geodetic CRS:  WGS 84
## Read the precipitation dataset
(precip <- raster("../../datos/precip_boyaca_X2019.07.tif"))
## class      : RasterLayer 
## dimensions : 59, 66, 3894  (nrow, ncol, ncell)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -74.66667, -71.91667, 4.625, 7.083333  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : precip_boyaca_X2019.07.tif 
## names      : precip_boyaca_X2019.07 
## values     : 18, 399  (min, max)

Note a few things about this object:

## change attribute name
names(precip) <- "X2019.07"
library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("red", "orange", "#fcc000","yellow", "cyan", "blue", "#3240cd"), values(precip$X2019.07), 
                    na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(precip$X2019.07 , colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(precip$X2019.07),
    title = "Rainfall-Jul.2019 [mm]")

3. Pre-processing raster data

We need to convert the raster data into point data. But, first, let’s mask the data to cover only the Boyacá departament:

(precip.mask <- mask(x = precip, mask = boyaca))
## class      : RasterLayer 
## dimensions : 59, 66, 3894  (nrow, ncol, ncell)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -74.66667, -71.91667, 4.625, 7.083333  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : X2019.07 
## values     : 23, 335  (min, max)

Let’s display the new raster layer:

leaflet() %>% addTiles() %>%
  addRasterImage(precip.mask$X2019.07 , colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(precip.mask$X2019.07),
    title = "Rainfall in Boyacá-Jul.2019 [mm]")

Now, the conversion from raster into points:

(precip.points <- rasterToPoints(precip.mask, spatial = TRUE))
## class       : SpatialPointsDataFrame 
## features    : 1083 
## extent      : -74.64583, -71.97917, 4.6875, 7.020833  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 1
## names       : X2019.07 
## min values  :       23 
## max values  :      335
names(precip.points) <- "rain"

Let’s remember that, for the interpolation task, we need project the data into map coordinates:

# crs=32618
(p.sf.magna <- st_transform(st_as_sf(precip.points), crs=32618))
(boyaca.sf.magna <- st_transform(boyaca, crs=32618))

Now, let’s convert the simple features into spatial data frames:

precip2 <- as(p.sf.magna, 'Spatial')
boyaca2 <- as(boyaca.sf.magna, 'Spatial')

Make sure the two extents match:

# Replace point boundary extent with that of Montes
precip2@bbox <- boyaca2@bbox

Let’s split the data into training and validation sets

train_index <- sample(1:nrow(precip2), 0.7 * nrow(precip2))
test_index <- setdiff(1:nrow(precip2), train_index)
ptos_train <- precip2[train_index, ]
ptos_test  <- precip2[test_index,]
ptos_train
## class       : SpatialPointsDataFrame 
## features    : 758 
## extent      : 539208.9, 833820.1, 518381.1, 777128.2  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## variables   : 1
## names       : rain 
## min values  :   23 
## max values  :  335

4. Inverse distance weighted interpolation

library(gstat)
# Create an empty grid where n is the total number of cells
grd   <- as.data.frame(spsample(precip2, "regular", n=319000))
# You need to figure out what is the expected size of the output grd
names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE  # Create SpatialPixel object
fullgrid(grd)    <- TRUE  # Create SpatialGrid object
# Define the CRS using an EPSG Code
x <- CRS(SRS_string='EPSG:32618')

# Display the stored CRS using comment()
#cat(comment(x), "\n")

# Store the wkt in a variable
wkt <- comment(x)

# Use this to assign the CRS of another sp-object
(y <- CRS(SRS_string = wkt))
## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
# Add P's projection information to the empty grid
sp::proj4string(grd) <- y
# Interpolate the grid cells using a power value of 2 (idp=2.0)
P.idw <- gstat::idw(rain ~ 1, ptos_train, newdata=grd, idp=2.0)
## [inverse distance weighted interpolation]
# Convert to raster object then clip to AOI
r <- raster(P.idw)
(r.m   <- raster::mask(r, boyaca2))
## class      : RasterLayer 
## dimensions : 532, 599, 318668  (nrow, ncol, ncell)
## resolution : 500.217, 500.217  (x, y)
## extent     : 537338.8, 836968.8, 514625.6, 780741  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : var1.pred 
## values     : 23.44648, 333.833  (min, max)
leaflet() %>% addTiles() %>%
  addRasterImage(r.m$var1.pred, colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(r.m$var1.pred),
    title = "IDW rainfall in July 2019 [mm]")

5. Kriging interpolation

This section is a shortened version of this notebook.

Let’s start by fitting the variogram model:

p.sf.magna
(the_crs <- crs(p.sf.magna))
## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
p.sf.magna$x <- st_coordinates(p.sf.magna)[,1]
p.sf.magna$y <- st_coordinates(p.sf.magna)[,2]
# Semivariance:
geoprec <- as.geodata(cbind(p.sf.magna$x,p.sf.magna$y,p.sf.magna$rain)) 
vario <- variog(geoprec, max.dist=300000)
## variog: computing omnidirectional variogram
fit <- variofit(vario)
## variofit: covariance model used is matern 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim
## Warning in variofit(vario): initial values not provided - running the default
## search
## variofit: searching for best initial value ... selected values:
##               sigmasq   phi     tausq     kappa
## initial.value "5505.06" "46106" "1101.01" "0.5"
## status        "est"     "est"   "est"     "fix"
## loss value: 161180927038.702
plot(vario) ; lines(fit)

# distance to closest other point:
d <- sapply(1:nrow(p.sf.magna), function(i) min(berryFunctions::distance(
                          p.sf.magna$x[i], p.sf.magna$y[i], p.sf.magna$x[-i], p.sf.magna$y[-i])))
hist(d/1000, breaks=20, main="distance to closest gauge [km]")

Now, kriging interpolation using Ordinary Kriging:

# Kriging:
res <- 1000 # 0.5 km, since stations are 4.6 km apart on average 
grid <- expand.grid(seq(min(p.sf.magna$x),max(p.sf.magna$x),res),
                    seq(min(p.sf.magna$y),max(p.sf.magna$y),res))
krico <- krige.control(type.krige="OK", obj.model=fit)
krobj <- krige.conv(geoprec, locations=grid, krige=krico)
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
grid_sf <- sf::st_as_sf(grid, coords=1:2, crs=sf::st_crs(boyaca.sf.magna))
isinp <- sapply(sf::st_within(grid_sf, boyaca.sf.magna), length) > 0
krobj2 <- krobj
krobj2$predict[!isinp] <- NA
# conversion from OK output to raster
(k_raster <- raster::rasterFromXYZ(cbind(grid, temp=krobj2$predict),
crs=y))
## class      : RasterLayer 
## dimensions : 259, 295, 76405  (nrow, ncol, ncell)
## resolution : 1000, 1000  (x, y)
## extent     : 538708.9, 833708.9, 517881.1, 776881.1  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : temp 
## values     : 23.01427, 320.7807  (min, max)
# interactive map
leaflet() %>% addTiles() %>%
  addRasterImage(k_raster$temp, colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(r.m$var1.pred),
    title = "OK rainfall in July 2019 [mm]")

6. Reproducibility

Cite as: Lizarazo, I., 2021. A tutorial on spatial interpolation. Available at https://rpubs.com/ials2un/interpol.

sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] RColorBrewer_1.1-2    leaflet_2.0.4.1       berryFunctions_1.20.1
##  [4] geoR_1.8-1            gstat_2.0-7           dplyr_1.0.7          
##  [7] rasterVis_0.50.3      latticeExtra_0.6-29   lattice_0.20-44      
## [10] terra_1.3-4           raster_3.4-13         rgdal_1.5-23         
## [13] sp_1.4-5              sf_1.0-2             
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.7              FNN_1.1.3               png_0.1-7              
##  [4] class_7.3-19            zoo_1.8-9               digest_0.6.27          
##  [7] utf8_1.2.2              R6_2.5.0                RandomFields_3.3.8     
## [10] evaluate_0.14           e1071_1.7-8             highr_0.9              
## [13] pillar_1.6.2            rlang_0.4.11            jquerylib_0.1.4        
## [16] hexbin_1.28.2           rmarkdown_2.9           stringr_1.4.0          
## [19] htmlwidgets_1.5.3       munsell_0.5.0           proxy_0.4-26           
## [22] compiler_4.1.0          xfun_0.24               base64enc_0.1-3        
## [25] pkgconfig_2.0.3         htmltools_0.5.1.1       tcltk_4.1.0            
## [28] tidyselect_1.1.1        tibble_3.1.3            intervals_0.15.2       
## [31] codetools_0.2-18        fansi_0.5.0             spacetime_1.2-5        
## [34] viridisLite_0.4.0       crayon_1.4.1            MASS_7.3-54            
## [37] grid_4.1.0              jsonlite_1.7.2          lifecycle_1.0.0        
## [40] DBI_1.1.1               magrittr_2.0.1          scales_1.1.1           
## [43] units_0.7-2             KernSmooth_2.23-20      stringi_1.7.3          
## [46] farver_2.1.0            bslib_0.2.5.1           ellipsis_0.3.2         
## [49] xts_0.12.1              generics_0.1.0          vctrs_0.3.8            
## [52] tools_4.1.0             glue_1.4.2              purrr_0.3.4            
## [55] splancs_2.01-42         crosstalk_1.1.1         jpeg_0.1-9             
## [58] abind_1.4-5             parallel_4.1.0          yaml_2.2.1             
## [61] colorspace_2.0-2        classInt_0.4-3          knitr_1.33             
## [64] RandomFieldsUtils_0.5.3 sass_0.4.0