Fixed rank kriging
Ivan Lizarazo
24.05.2020
1. Introduction
This notebook is a quick exploration of fixed rank kriging (FRK), a technique for spatial/spatio-temporal modelling and prediction with large datasets. This notebook is based on the vignette of the R package FRK developed by Andrew Zammit-Mangion and Noel Cressie (2020). This notebook aims to help Geomatica Basica students at Universidad Nacional de Colombia to understand R geospatial functionalities.
Its authors state that “FRK differs from existing packages for spatial modelling and prediction by avoiding stationary and isotropic covariance and variogram models, instead constructing a spatial random effects (SRE) model on a discretised spatial domain. The discrete element is known as a basic areal unit (BAU), whose introduction in the software leads to several practical advantages”.
A detailed technical description of FRK can be found here
Use of FRK is illustrated here using point precipitation data obtained from a CHIRPS dataset. More information here.
Preliminaries
To install the most recent (development) version, first please install INLA:
#install.packages("INLA", repos=c(getOption("repos"), INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)Then please load devtools and type:
#library(devtools)
#install_github("andrewzm/FRK", dependencies = TRUE, build_vignettes = TRUE)Now, load the packages
## Loading required package: sp## rgdal: version: 1.4-8, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
## Linking to sp version: 1.3-2## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.0 ✓ purrr 0.3.4
## ✓ tibble 3.0.1 ✓ dplyr 0.8.5
## ✓ tidyr 1.0.3 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks raster::select()##
## Attaching package: 'FRK'## The following object is masked from 'package:raster':
##
## distance2. Reading point data
Let’s read the precipitation point dataset. It has been previously saved as a local file in geojson format.
precip.points <- geojsonio::geojson_read("./chirps/ppoints.geojson", what="sp")Let’s check what is the precip.points object:
precip.points## class : SpatialPointsDataFrame
## features : 213
## extent : -75.675, -74.825, 9.274999, 10.075 (xmin, xmax, ymin, ymax)## Warning in proj4string(x): CRS object has comment, which is lost in output## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 1
## names : rain
## min values : 1.06151795387268
## max values : 5.76676845550537# conversion
p.sf <- st_as_sf(precip.points)Let’s read the shapefile with the area of interest’s boundaries:
(aoi <- shapefile("/Users/ivanlizarazo/Documents/ivan/pr/montes/MontesdeMaria_WGS84.shp"))## class : SpatialPolygonsDataFrame
## features : 15
## extent : -75.70418, -74.77214, 9.22535, 10.14548 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 11
## names : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR, MPIO_CRSLC, MPIO_NAREA, MPIO_NANO, DPTO_CNMBR, Shape_Leng, Shape_Area, layer, path
## min values : 13, 13212, CHALÁN, 1780, 83.71001417, 2017, BOLÍVAR, 0.5720244142, 0.00689090159, bol_montes, /Users/ivanlizarazo/Documents/ivan/PR/montes/bol_montes.shp|layername=bol_montes
## max values : 70, 70823, ZAMBRANO, Ordenanza 6 de Octubre 30 de 1968, 1035.07793303, 2017, SUCRE, 2.08052417308, 0.08525488217, suc_montes, /Users/ivanlizarazo/Documents/ivan/PR/montes/suc_montes.shp|layername=suc_montesWe need to convert it to a spatial feature:
montes_sf <- sf::st_as_sf(aoi)Now, dissolve the internal boundaries:
(border_sf <-
montes_sf %>%
summarise(area = sum(MPIO_NAREA)))## Simple feature collection with 1 feature and 1 field
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -75.70418 ymin: 9.22535 xmax: -74.77214 ymax: 10.14548
## CRS: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## area geometry
## 1 6513.845 POLYGON ((-75.11082 9.44062...Two reprojection tasks:
p.sf.magna <- st_transform(p.sf, crs=3116)montes.sf.magna <- st_transform(montes_sf, crs=3116)A conversion:
(precip2 <- as(p.sf.magna, 'Spatial'))## class : SpatialPointsDataFrame
## features : 213
## extent : 824685.3, 918049.7, 1517694, 1606178 (xmin, xmax, ymin, ymax)## Warning in proj4string(x): CRS object has comment, which is lost in output## crs : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## variables : 1
## names : rain
## min values : 1.06151795387268
## max values : 5.76676845550537A new name for an attribute:
names(precip2) <- "rain"Check the output:
precip2## class : SpatialPointsDataFrame
## features : 213
## extent : 824685.3, 918049.7, 1517694, 1606178 (xmin, xmax, ymin, ymax)## Warning in proj4string(x): CRS object has comment, which is lost in output## crs : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## variables : 1
## names : rain
## min values : 1.06151795387268
## max values : 5.76676845550537Let’s write the precipitation dataset (just in case anything goes wrong):
#shapefile(precip2, filename='./chirps/precip2.shp', overwrite=TRUE)If needed, read it. Otherwise, keep going
Let’s round the precipitation values:
precip2$rainfall <- round(precip2$rain, 1)
##precip2$rainfall <- round(precip2$chirps.v2.0.2020.03.6, 1)Check what we got:
precip2## class : SpatialPointsDataFrame
## features : 213
## extent : 824685.3, 918049.7, 1517694, 1606178 (xmin, xmax, ymin, ymax)## Warning in proj4string(x): CRS object has comment, which is lost in output## crs : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## variables : 2
## names : rain, rainfall
## min values : 1.06151795387268, 1.1
## max values : 5.76676845550537, 5.8Again, object conversion:
(montes2 <- as(montes.sf.magna, 'Spatial'))## class : SpatialPolygonsDataFrame
## features : 15
## extent : 821468.2, 923749.4, 1512192, 1614054 (xmin, xmax, ymin, ymax)## Warning in proj4string(x): CRS object has comment, which is lost in output## crs : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## variables : 11
## names : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR, MPIO_CRSLC, MPIO_NAREA, MPIO_NANO, DPTO_CNMBR, Shape_Leng, Shape_Area, layer, path
## min values : 13, 13212, CHALÁN, 1780, 83.71001417, 2017, BOLÍVAR, 0.5720244142, 0.00689090159, bol_montes, /Users/ivanlizarazo/Documents/ivan/PR/montes/bol_montes.shp|layername=bol_montes
## max values : 70, 70823, ZAMBRANO, Ordenanza 6 de Octubre 30 de 1968, 1035.07793303, 2017, SUCRE, 2.08052417308, 0.08525488217, suc_montes, /Users/ivanlizarazo/Documents/ivan/PR/montes/suc_montes.shp|layername=suc_montesWrite the new dataset if you want to have a backup:
#shapefile(montes2, filename='./chirps/montes2.shp', overwrite=TRUE)Make sure the two objects’ extents match:
# Replace point boundary extent with that of Montes
precip2@bbox <- montes2@bboxPlotting the input data with tmap
tm_shape(montes2) + tm_polygons() +
tm_shape(precip2) +
tm_dots(col="rainfall", palette = "RdBu", midpoint = 3.0,
title="Sampled precipitation \n(in mm)", size=0.2) +
tm_text("rainfall", just="center", xmod=.6, size = 0.5) +
tm_legend(legend.outside=TRUE)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
3. Interpolating precipitation data
Let’s run the FRK interpolation with no manual intervention:
## Run FRK
S <- FRK(f = rainfall~1, # Formula to FRK
list(precip2), # All datasets are supplied in list
n_EM = 30) # Max number of EM iterations## Warning in proj4string(data): CRS object has comment, which is lost in output## Loading required namespace: INLA## Warning in proj4string(d): CRS object has comment, which is lost in output## Warning in proj4string(x): CRS object has comment, which is lost in output## Warning in proj4string(sp_pts): CRS object has comment, which is lost in output## Warning in proj4string(data_sp): CRS object has comment, which is lost in outputPred <- SRE.predict(SRE_model = S) # Prediction stage## Warning in SRE.predict(SRE_model = S): SRE.predict is deprecated. Please use
## predict.The object Pred contains the prediction vector, the square of the prediction standard error, and the prediction standard error at the BAU level in the fields mu, var, and sd respectively.
summary(Pred)## Object of class SpatialPixelsDataFrame
## Coordinates:
## min max
## coords.x1 805545.6 937189.4
## coords.x2 1499554.2 1624317.8
## Is projected: TRUE
## proj4string :
## [+proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1
## +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
## +no_defs]
## Number of points: 8861
## Grid attributes:
## cellcentre.offset cellsize cells.dim
## coords.x1 806012.5 933.6436 141
## coords.x2 1499996.6 884.8478 141
## Data attributes:
## coords.x1 coords.x2 fs var
## Min. :820017 Min. :1513269 Min. :1 Min. :0.02329
## 1st Qu.:854562 1st Qu.:1547778 1st Qu.:1 1st Qu.:0.02777
## Median :875102 Median :1568130 Median :1 Median :0.03003
## Mean :875784 Mean :1567485 Mean :1 Mean :0.03599
## 3rd Qu.:897510 3rd Qu.:1588481 3rd Qu.:1 3rd Qu.:0.03435
## Max. :922718 Max. :1610603 Max. :1 Max. :0.19526
## mu sd
## Min. :1.199 Min. :0.1526
## 1st Qu.:2.530 1st Qu.:0.1666
## Median :3.414 Median :0.1733
## Mean :3.352 Mean :0.1862
## 3rd Qu.:4.146 3rd Qu.:0.1853
## Max. :5.596 Max. :0.4419Note the position of the layers containing the prediction (5) and the standard error
xy <- data.frame(coordinates(Pred)) # Extract info from predictionshead(xy)## coords.x1 coords.x2
## 2d871998be5df427fd4062b5342eb9b7 846159.1 1513269
## f262f98cc558e3925ace1941bf943028 847092.8 1513269
## de4842b70120fbc594e58a59eea74234 848026.4 1513269
## b4cf734ee7b115313e9a1668ee111cdd 848960.1 1513269
## 8974789bbd522f2afbdd8c60b20549a1 849893.7 1513269
## 077f56bc33f32a3e8867f12053c1025b 850827.4 1513269class(Pred)## [1] "SpatialPixelsDataFrame"
## attr(,"package")
## [1] "sp"Convert the SpatialPixelsDataFrame into a RasterLayer:
(rPred <- raster(Pred, layer=5, values=TRUE))## class : RasterLayer
## dimensions : 141, 141, 19881 (nrow, ncol, ncell)
## resolution : 933.6436, 884.8478 (x, y)
## extent : 805545.6, 937189.4, 1499554, 1624318 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## source : memory
## names : mu
## values : 1.198822, 5.595727 (min, max)names(rPred) <- "rain"xy$precip <- Pred$mu
xy$se <- Pred$sdConvert the points spatial data frame to a data frame for plotting purposes:
zdf <- as.data.frame(precip2)head(zdf)## rain rainfall coords.x1 coords.x2
## 1 3.638230 3.6 863247.2 1606178
## 2 3.631608 3.6 868729.1 1606158
## 3 3.505968 3.5 874210.9 1606138
## 4 4.947118 4.9 879692.6 1606119
## 5 5.308706 5.3 885174.2 1606102
## 6 5.208715 5.2 890655.7 1606084These are the points used as input:
## Plotting
ggplot(zdf) + geom_point(aes(coords.x1,coords.x2,colour=rainfall)) +
scale_colour_distiller(palette="Spectral") + theme_bw() + coord_fixed()
This the output of the RFK algoritm (i.e. the interpolated precipitation surface):
ggplot(xy) + geom_raster(aes(coords.x1,coords.x2,fill=precip)) +
scale_fill_distiller(palette="Spectral") + theme_bw() + coord_fixed()
ggplot(xy) + geom_tile(aes(coords.x1,coords.x2,fill=se)) +
geom_point(data=zdf,aes(coords.x1,coords.x2),pch=46) +
scale_fill_distiller(palette="Spectral") + theme_bw() + coord_fixed()
head(xy$precip)## [1] 2.391815 2.415911 2.446010 2.480465 2.516934 2.554164head(xy$mu)## NULL# Plot the map
tm_shape(rPred) +
tm_raster(n=10, palette="RdBu", auto.palette.mapping=FALSE,
title="Fixed Rank Kriging\nPredicted precipitation \n(in mm)") +
tm_shape(precip2) + tm_dots(size=0.2) +
tm_legend(legend.outside=TRUE)
A better visualization of the interpolated surface:
library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("red", "orange", "yellow", "cyan", "blue"), values(rPred),
na.color = "transparent")
leaflet() %>% addTiles() %>%
addRasterImage(rPred, colors = pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(rPred),
title = "Fixed Rank Kriging Interpolation [mm]")Let’s convert the standard error values of the Pred object into a raster layer:
(rError <- raster(Pred, layer=6, values=TRUE))## class : RasterLayer
## dimensions : 141, 141, 19881 (nrow, ncol, ncell)
## resolution : 933.6436, 884.8478 (x, y)
## extent : 805545.6, 937189.4, 1499554, 1624318 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## source : memory
## names : sd
## values : 0.1526234, 0.4418811 (min, max)A better visualization of the uncertainty on the estimation:
library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("green", "yellow", "red"), values(rError),
na.color = "transparent")
leaflet() %>% addTiles() %>%
addRasterImage(rError, colors = pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(rError),
title = "Uncertainty in fixed rank kriging [mm]")You may compare your own FRK results with the ones obtained using OK and IDW in your first interpolation notebook.
sessionInfo()## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/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.3 FRK_0.2.3 gstat_2.0-6
## [5] tmap_3.0 forcats_0.5.0 stringr_1.4.0 dplyr_0.8.5
## [9] purrr_0.3.4 readr_1.3.1 tidyr_1.0.3 tibble_3.0.1
## [13] ggplot2_3.3.0 tidyverse_1.3.0 sf_0.9-3 raster_3.1-5
## [17] rgdal_1.4-8 sp_1.4-2
##
## loaded via a namespace (and not attached):
## [1] leafem_0.1.1 colorspace_1.4-1 ellipsis_0.3.1
## [4] class_7.3-17 htmlTable_1.13.3 base64enc_0.1-3
## [7] fs_1.4.1 dichromat_2.0-0 httpcode_0.3.0
## [10] rstudioapi_0.11 farver_2.0.3 fansi_0.4.1
## [13] lubridate_1.7.8 xml2_1.3.2 codetools_0.2-16
## [16] splines_3.6.3 knitr_1.28 spam_2.5-1
## [19] Formula_1.2-3 jsonlite_1.6.1 tmaptools_3.0
## [22] broom_0.5.6 cluster_2.1.0 dbplyr_1.4.3
## [25] rgeos_0.5-3 png_0.1-7 compiler_3.6.3
## [28] httr_1.4.1 backports_1.1.7 lazyeval_0.2.2
## [31] assertthat_0.2.1 Matrix_1.2-18 cli_2.0.2
## [34] acepack_1.4.1 htmltools_0.4.0 tools_3.6.3
## [37] dotCall64_1.0-0 gtable_0.3.0 glue_1.4.1
## [40] geojson_0.3.2 V8_3.0.2 Rcpp_1.0.4.6
## [43] cellranger_1.1.0 vctrs_0.3.0 crul_0.9.0
## [46] nlme_3.1-147 leafsync_0.1.0 crosstalk_1.1.0.1
## [49] lwgeom_0.2-3 xfun_0.14 rvest_0.3.5
## [52] lifecycle_0.2.0 XML_3.99-0.3 jqr_1.1.0
## [55] zoo_1.8-8 scales_1.1.1 sparseinv_0.1.3
## [58] hms_0.5.3 parallel_3.6.3 curl_4.3
## [61] yaml_2.2.1 gridExtra_2.3 rpart_4.1-15
## [64] latticeExtra_0.6-29 stringi_1.4.6 maptools_0.9-9
## [67] e1071_1.7-3 checkmate_2.0.0 intervals_0.15.2
## [70] rlang_0.4.6 pkgconfig_2.0.3 evaluate_0.14
## [73] lattice_0.20-41 labeling_0.3 htmlwidgets_1.5.1
## [76] tidyselect_1.1.0 plyr_1.8.6 magrittr_1.5
## [79] R6_2.4.1 geojsonio_0.9.2 generics_0.0.2
## [82] Hmisc_4.4-0 DBI_1.1.0 pillar_1.4.4
## [85] haven_2.2.0 foreign_0.8-75 withr_2.2.0
## [88] units_0.6-6 stars_0.4-2 xts_0.12-0
## [91] survival_3.1-12 abind_1.4-5 nnet_7.3-14
## [94] INLA_20.03.17 spacetime_1.2-3 modelr_0.1.6
## [97] crayon_1.3.4 KernSmooth_2.23-17 rmarkdown_2.1
## [100] jpeg_0.1-8.1 grid_3.6.3 readxl_1.3.1
## [103] data.table_1.12.8 FNN_1.1.3 splancs_2.01-40
## [106] reprex_0.3.0 digest_0.6.25 classInt_0.4-3
## [109] munsell_0.5.0 viridisLite_0.3.0