Interpolación espacial de datos climáticos

Introducción

Las librerias stars, geostat y automap sirven como alternativa a la libreria geoR para realizar la interpolación espacial de datos puntuales.

Librerias

library(sf)

Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE

library(stars)

Loading required package: abind

library(gstat)
library(automap)
library(leaflet)
library(leafem)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:raster’:

    intersect, select, union

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

Lectura del shapefile del área de estudio

mun.tmp = st_read("C:\\Users\\ynata\\OneDrive\\Documentos\\GGB2022\\Datos\\DPTOCASANARE.shp")

Reading layer `DPTOCASANARE' from data source 
  `C:\Users\ynata\OneDrive\Documentos\GGB2022\Datos\DPTOCASANARE.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 19 features and 11 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -73.07777 ymin: 4.287476 xmax: -69.83591 ymax: 6.346111
Geodetic CRS:  WGS 84

mun.tmp %>% dplyr::select( MPIO_CCNCT, MPIO_CNMBR) -> casanare
rename(casanare, MPIO_CCDGO = MPIO_CCNCT)

Simple feature collection with 19 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -73.07777 ymin: 4.287476 xmax: -69.83591 ymax: 6.346111
Geodetic CRS:  WGS 84
First 10 features:
   MPIO_CCDGO MPIO_CNMBR                       geometry
1       85001      YOPAL POLYGON ((-72.39513 5.56853...
2       85010    AGUAZUL POLYGON ((-72.56545 5.36972...
3       85015    CHÁMEZA POLYGON ((-72.81017 5.36659...
4       85136  LA SALINA POLYGON ((-72.33885 6.34470...
5       85139       MANÍ POLYGON ((-72.34155 5.06495...
6       85162  MONTERREY POLYGON ((-72.89989 5.03465...
7       85225    NUNCHÍA POLYGON ((-72.19558 5.71924...
8       85230     OROCUÉ POLYGON ((-71.50965 5.21730...
9       85263       PORE POLYGON ((-72.04587 5.83819...
10      85279    RECETOR POLYGON ((-72.80501 5.38092...

Lectura del raster de precipitación

archivo <- ("C:\\Users\\ynata\\OneDrive\\Documentos\\GGB2022\\R\\precip.tif")
(precip <- read_stars(archivo))

stars object with 2 dimensions and 1 attribute
attribute(s):
            Min. 1st Qu. Median     Mean 3rd Qu.  Max.
precip.tif   0.2     2.4    7.2 14.30432    15.7 161.1
dimension(s):

precip.mask <- st_crop(precip, casanare)

precip.mask 

stars object with 2 dimensions and 1 attribute
attribute(s):
            Min. 1st Qu. Median     Mean 3rd Qu.  Max. NA's
precip.tif   0.2     1.7    3.1 6.656184     6.1 122.2 1892
dimension(s):

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      precip.mask,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("red", "orange", "#fcc000","yellow", "cyan", "blue", "#3240cd"), 
                                  domain = 4:83)
    ) 
m

Vectorización de la capa de precipitación

precip.pts <- st_as_sf(precip.mask, as_points = TRUE, merge = FALSE, long = TRUE)

precip.pts

Simple feature collection with 2086 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -73.0625 ymin: 4.270833 xmax: -69.85417 ymax: 6.354167
Geodetic CRS:  WGS 84
First 10 features:
   precip.tif                   geometry
1        26.8 POINT (-72.39583 6.270833)
2         0.5 POINT (-71.27083 6.270833)
3         0.5 POINT (-71.22917 6.270833)
4         0.5  POINT (-71.1875 6.270833)
5         1.1  POINT (-70.4375 6.270833)
6         1.0 POINT (-70.39583 6.270833)
7         0.9 POINT (-70.35417 6.270833)
8         1.3  POINT (-70.3125 6.270833)
9         1.1 POINT (-70.27083 6.270833)
10        1.0 POINT (-70.22917 6.270833)

Datos de entrenamiento y prueba

dt = sort(sample(2086, 2086*.7))
train<-precip.pts[dt,]
test<-precip.pts[-dt,]

train

Simple feature collection with 1460 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -73.0625 ymin: 4.3125 xmax: -69.85417 ymax: 6.270833
Geodetic CRS:  WGS 84
First 10 features:
   precip.tif                   geometry
1        26.8 POINT (-72.39583 6.270833)
2         0.5 POINT (-71.27083 6.270833)
3         0.5 POINT (-71.22917 6.270833)
6         1.0 POINT (-70.39583 6.270833)
7         0.9 POINT (-70.35417 6.270833)
10        1.0 POINT (-70.22917 6.270833)
11        1.4  POINT (-70.1875 6.270833)
13       23.1 POINT (-72.35417 6.229167)
14        0.3  POINT (-71.4375 6.229167)
15        0.3 POINT (-71.39583 6.229167)

Visualización de lugares de muestreo

longit <- st_coordinates(train)[,1]
latit <- st_coordinates(train)[,2]
rain <- train$precip.tif

id <- seq(1,1460,1) 

(sitios <- data.frame(id, longit, latit, rain))

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      precip.mask,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("red", "orange", "#fcc000","yellow", "cyan", "blue", "#3240cd"), 
                                  domain = 4:83)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$rain, clusterOptions = markerClusterOptions())
m

Lectura del modelo digital de elevación

## file.choose()

(dem = read_stars("C:\\Users\\ynata\\OneDrive\\Documentos\\GGB2022\\Datos\\ELV_CASANARE.tif"))

stars object with 2 dimensions and 1 attribute
attribute(s):
                  Min. 1st Qu. Median     Mean 3rd Qu. Max.  NA's
ELV_CASANARE.tif    84     128    157 259.2479     197 4208 43383
dimension(s):

pal.dem <- colorNumeric(palette = c("forestgreen", "green", "yellow", "brown", "lightcyan"), domain = 84:4208 , na.color = "transparent")

m <- leaflet() %>%
  addTiles() %>%  
  leafem:::addGeoRaster(
      dem,
      opacity = 0.7,                
      colorOptions = colorOptions(palette = c("forestgreen", "green", "yellow",  "brown",         "lightcyan"), 
                                  domain = 84:4208)
    ) %>%
  addMarkers(lng=sitios$longit,lat=sitios$latit, popup=sitios$rain, clusterOptions = markerClusterOptions()) %>%
  addLegend("bottomright", pal = pal.dem, values= dem$ELV_CASANARE.tif,
    title = "Elevation")

m

Extracción de los valores de elavación en los sitios de entrenamiento

st_crs(dem) <- st_crs(train)

(train = st_join(train, st_as_sf(dem)))

Simple feature collection with 1460 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -73.0625 ymin: 4.3125 xmax: -69.85417 ymax: 6.270833
Geodetic CRS:  WGS 84
First 10 features:
   precip.tif ELV_CASANARE.tif                   geometry
1        26.8             3965 POINT (-72.39583 6.270833)
2         0.5              152 POINT (-71.27083 6.270833)
3         0.5              153 POINT (-71.22917 6.270833)
6         1.0              105 POINT (-70.39583 6.270833)
7         0.9              104 POINT (-70.35417 6.270833)
10        1.0              101 POINT (-70.22917 6.270833)
11        1.4              102  POINT (-70.1875 6.270833)
13       23.1             2648 POINT (-72.35417 6.229167)
14        0.3              166  POINT (-71.4375 6.229167)
15        0.3              165 POINT (-71.39583 6.229167)

train = train[!is.na(train$ELV_CASANARE.tif), ]

names(train) <- c("precip", "elev", "geometry")

train

Simple feature collection with 1460 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -73.0625 ymin: 4.3125 xmax: -69.85417 ymax: 6.270833
Geodetic CRS:  WGS 84
First 10 features:
   precip elev                   geometry
1    26.8 3965 POINT (-72.39583 6.270833)
2     0.5  152 POINT (-71.27083 6.270833)
3     0.5  153 POINT (-71.22917 6.270833)
6     1.0  105 POINT (-70.39583 6.270833)
7     0.9  104 POINT (-70.35417 6.270833)
10    1.0  101 POINT (-70.22917 6.270833)
11    1.4  102  POINT (-70.1875 6.270833)
13   23.1 2648 POINT (-72.35417 6.229167)
14    0.3  166  POINT (-71.4375 6.229167)
15    0.3  165 POINT (-71.39583 6.229167)

Interpolación espacial

Creación del objeto gstat

Para interpolar, primero tenemos que crear un objeto de la clase gstat, utilizando una función del mismo nombre: gstat.

Un objeto gstat contiene toda la información necesaria para realizar la interpolación espacial, a saber.

• La definición del modelo
• Los datos de calibración

Basándose en sus argumentos, la función gstat “entiende” qué tipo de modelo de interpolación queremos utilizar:

• Sin modelo de variograma → IDW
• Modelo de variograma, sin covariables → Kriging ordinario
• Modelo de variograma, con covariables → Kriging universal

Vamos a utilizar tres parámetros de la función gstat:

• fórmula: la “fórmula” de predicción que especifica las variables dependientes e independientes (también llamadas covariables)
• datos: los datos de calibración (es decir, los datos del tren)
• modelo: el modelo del variograma

IDW

g1 = gstat(formula = precip ~ 1, data = train)

Ahora que nuestro modelo de interpolación g1 está definido, podemos utilizar la función predecir para interpolar realmente, es decir, para estimar los valores de precipitación.

La función predecir acepta:

• Un objeto raster-stars, como dem
• Un objeto model-gstat, como g1

El raster sirve para dos propósitos:

• Especificar las ubicaciones en las que queremos hacer predicciones (en todos los métodos)
• Especificar los valores de las covariables (sólo en Kriging Universal)

Por ejemplo, la siguiente expresión interpola los valores de las precipitaciones según el modelo definido en g1 y el modelo de trama definido en dem:

 ## [inverse distance weighted interpolation]
z1 = predict(g1, dem)

z1

stars object with 2 dimensions and 2 attributes
attribute(s):
           Min.  1st Qu.   Median     Mean  3rd Qu.  Max.  NA's
var1.pred   0.3 2.375452 3.839834 6.464584 6.091785 118.3 43383
var1.var     NA       NA       NA      NaN       NA    NA 95448
dimension(s):

Podemos subconjuntar sólo el primer atributo y renombrarlo como “precipitación”:

names(z1) = "precipitation"
z1 = z1["precipitation",,]

En la siguiente figura se muestra el ráster de precipitaciones interpolado mediante IDW:

b = seq(1, 100, 1)
plot(z1, breaks = b, col = hcl.colors(length(b)-1, "Spectral"), reset = FALSE)
plot(st_geometry(train), pch = 3, add = TRUE)

contour(z1, breaks = b, add = TRUE)

Kriging ordinario Los métodos de Kriging requieren un modelo de variograma. El modelo de variograma es una forma objetiva de cuantificar el patrón de autocorrelación en los datos, y asignar pesos en consecuencia al hacer predicciones.

Como primer paso, podemos calcular y examinar el variograma empírico utilizando la función variograma.

La función requiere dos argumentos

• Fórmula: especifica la variable dependiente y las covariables, al igual que en gstat
• Datos: la capa de puntos con la variable dependiente y las covariables como atributos de los puntos

Por ejemplo, la siguiente expresión calcula el variograma empírico, sin covariables:

v_emp_ok = variogram(precip ~ 1, data=train)
plot(v_emp_ok)

Hay varias formas de ajustar un modelo de variograma a un variograma empírico. Utilizaremos la más sencilla: el ajuste automático mediante la función autofitVariogram del paquete automap:

v_mod_ok = autofitVariogram(precip ~ 1, as(train, "Spatial"))

La función elige el tipo de modelo que mejor se ajusta, y también afina sus parámetros. Puede utilizar show.vgms() para mostrar los tipos de modelos de variograma.

Tenga en cuenta que la función autofitVariogram no funciona en objetos sf, por lo que convertimos el objeto en un SpatialPointsDataFrame (paquete sp).

El modelo ajustado puede representarse con plot:

plot(v_mod_ok)

El objeto resultante es en realidad una lista con varios componentes, incluyendo el variograma empírico y el modelo de variograma ajustado. El componente $var_model del objeto resultante contiene el modelo real:

v_mod_ok$var_model

NA

El modelo del variograma puede entonces pasarse a la función gstat, y podemos continuar con la interpolación Kriging ordinaria:

## [using ordinary kriging]
g2 = gstat(formula = precip ~ 1, model = v_mod_ok$var_model, data = train)
z2= predict(g2, dem)

Una vez más, vamos a subconjuntar el atributo de valores predichos y a renombrarlo:

names(z2) = "precip"

Las predicciones de Kriging ordinario se muestran en la siguiente figura:

b = seq(1, 100, 1)
plot(z2, breaks = b, col = hcl.colors(length(b)-1, "Spectral"), reset = FALSE)
plot(st_geometry(train), pch = 3, add = TRUE)

contour(z2, breaks = b, add = TRUE)

Universal Kriging

v_emp_uk = variogram(precip ~ elev, train)
train.sp = as(train, "Spatial")
v_mod_uk = autofitVariogram(precip ~ elev, train.sp)

plot(v_mod_uk)

plot(v_emp_ok, model = v_mod_ok$var_model, ylim = c(0, 150), main = "OK")

plot(v_emp_uk, model = v_mod_uk$var_model, ylim = c(0, 70), main = "UK")

g3 = gstat(formula = precip ~ elev, model = v_mod_uk$var_model, data = train.sp)

names(dem)

[1] "ELV_CASANARE.tif"

names(dem) <- "elev"

z3 = predict(g3, dem)

z3 = z3["var1.pred",,]
names(z3) = "precipitation"

b = seq(1, 100, 1)
plot(z3, breaks = b, col = hcl.colors(length(b)-1, "Spectral"), reset = FALSE)
plot(st_geometry(train), pch = 3, add = TRUE)

contour(z3, breaks = b, add = TRUE)

4. Evaluación cualitativa de los resultados

paleta <- colorNumeric(palette = c("coral2", "coral1","coral", "orange", "#fcc000","yellow", "cyan", "blue", "#3240cd"), domain = 10:100, na.color = "transparent")

colores <- colorOptions(palette = c("red", "orange", "#fcc000","yellow", "cyan", "blue", "#3240cd"), domain = 10:100, na.color = "transparent")

m <- leaflet() %>%
  addTiles() %>%  
  addGeoRaster(z1, opacity = 0.7, colorOptions = colores, group="IDW") %>%
  addGeoRaster(z2, colorOptions = colores, opacity = 0.7, group= "OK")  %>%
  addGeoRaster(z3, colorOptions = colores, opacity = 0.7, group= "UK")  %>%
  # Add layers controls
  addLayersControl(
    overlayGroups = c("UK", "OK", "IDW"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>% 
  
    addLegend("bottomright", pal = paleta, values= z1$precipitation,
    title = "Precipitation (2019)"
  )

Warning in pal(c(r[1], cuts, r[2])) :
  Some values were outside the color scale and will be treated as NA
Warning in write_stars.stars(x, dsn = fl) :
  all but first attribute are ignored

m

5. Validación cruzada

bubble(as(cv1[, "residual"], "Spatial"))

sqrt(sum((cv1$var1.pred - cv1$observed)^2) / nrow(cv1))

[1] 5.564086

bubble(as(cv2[, "residual"], "Spatial"))

sqrt(sum((cv2$var1.pred - cv2$observed)^2) / nrow(cv2))

[1] 2.287359

bubble(as(cv3[, "residual"], "Spatial"))

sqrt(sum((cv3$var1.pred - cv3$observed)^2) / nrow(cv3))

[1] 2.661829