R Notebook

4.5 Kriging Ajustar el modelo de variograma Primero, necesitamos crear un modelo de variograma. Tenga en cuenta que el modelo de variograma se calcula sobre los datos de tendencia. Esto se implementa en el siguiente fragmento de código pasando el modelo de tendencia de primer orden (definido en un fragmento de código anterior como objeto de fórmula f.1) a la función variograma.

4.3 Ajuste polinomial de primer orden Para ajustar un modelo polinomial de primer orden de la forma precip.

# Define the 1st order polynomial equation
f.1 <- as.formula(precipitaciones ~ X + Y) 
 
# Add X and Y to P
P$X <- coordinates(P)[,1]
P$Y <- coordinates(P)[,2]

# Run the regression model
lm.1 <- lm( f.1, data=P)

# Use the regression model output to interpolate the surface
dat.1st <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.1, newdata=grd))) 
# Clip the interpolated raster to Texas
r   <- raster(dat.1st)

CRS object has comment, which is lost in outputDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

r.m <- raster::mask(r, meta2)

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

# Plot the map
tm_shape(r.m) + 
  tm_raster(n=10, palette="RdBu", auto.palette.mapping=FALSE, 
            title="1st order polinomial fir \nPredicted precipitation \n(in mm)") +
  tm_shape(P) + tm_dots(size=0.2) +
  tm_legend(legend.outside=TRUE)

The argument auto.palette.mapping is deprecated. Please use midpoint for numeric data and stretch.palette for categorical data to control the palette mapping.CRS object has comment, which is lost in output

Ahora si, kriging

# Definir la ecuación polinómica de primer orden
f.1 <- as.formula (precipitaciones ~ X + Y)

# Calcular el variograma de muestra; tenga en cuenta que el modelo de tendencia f.1 es uno de los
# parámetros pasados a variogram (). Esto le dice a la función que cree el
# variograma en los datos de tendencia.
var.smpl <- variogram(f.1, P, cloud = FALSE, cutoff=200000, width=10000)

# Calcular el modelo de variograma pasando los valores de pepita, umbral y rango
# to fit.variogram () a través de la función vgm ().
dat.fit  <- fit.variogram(var.smpl, fit.ranges = TRUE, fit.sills = TRUE,
                          vgm(psill=5, model="Mat", range=120000, nugget=-0.2))

## modelo puede ser Exp, Gau, Sph, Mat
## Consulte el paquete gstat - función vgm - para obtener los nombres de estos modelos
## Revise la teoría de kriging para comprender su significado

dat.fit

# La siguiente plot nos permite evaluar el ajuste
plot(var.smpl, dat.fit, xlim=c(0,150000))

Generar superficie Kriged Luego, se usa el modelo de variograma dat.fit para generar una superficie interpolada kriged. La función krige nos permite incluir el modelo de tendencia, lo que nos evita tener que reducir la tendencia de los datos, corregir los residuos y luego combinar los dos rásteres. En cambio, todo lo que tenemos que hacer es pasar krige la fórmula de tendencia f.1.

# Definir el modelo de tendencia
f.1 <- as.formula (precipitaciones ~ X + Y)

# Realice la interpolación de krige (observe el uso del modelo de variograma
# creado en el paso anterior)
dat.krg <- krige (f.1, P, grd, dat.fit)

CRS object has comment, which is lost in outputCRS object has comment, which is lost in output

[using universal kriging]

CRS object has comment, which is lost in outputDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

# Convertir superficie kriged en un objeto ráster para el recorte
r <- raster(dat.krg)

CRS object has comment, which is lost in outputDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

r.m <- raster::mask(r, meta2)

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

r.m

class      : RasterLayer 
dimensions : 652, 767, 500084  (nrow, ncol, ncell)
resolution : 543.8037, 543.8037  (x, y)
extent     : 911335.7, 1328433, 671186.5, 1025747  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +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      : var1.pred 
values     : 0.8829389, 34.13062  (min, max)

# Plot the map
tm_shape(r.m) + 
  tm_raster(n=10, palette="RdBu", auto.palette.mapping=FALSE, 
            title="Kriging universal \nPredicted precipitation \n(in mm)") +
  tm_shape(P) + tm_dots(size=0.2) +
  tm_legend(legend.outside=TRUE)

The argument auto.palette.mapping is deprecated. Please use midpoint for numeric data and stretch.palette for categorical data to control the palette mapping.CRS object has comment, which is lost in output

Una mejor visualización

library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("red", "orange", "yellow", "blue", "darkblue"), values(precip.mask),
  na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(r.m, colors = pal, opacity = 0.6) %>%
  addLegend(pal = pal, values = values(r.m),
    title = "Precipitación Interpolación Kriging en Meta\n Del  26.03 al 30.03 en 2020 [mm]")

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedUsing PROJ not WKT2 stringsDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedUsing PROJ not WKT2 stringsDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedUsing PROJ not WKT2 stringsDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsSome values were outside the color scale and will be treated as NASome values were outside the color scale and will be treated as NA

Generar los mapas de varianza e intervalo de confianza. El objeto dat.krg almacena no solo los valores interpolados, sino también los valores de varianza. Estos se pueden pasar al objeto ráster para la asignación de la siguiente manera:

r   <- raster(dat.krg, layer="var1.var")

CRS object has comment, which is lost in outputDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

r.m <- raster::mask(r, meta2)

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

tm_shape(r.m) + 
  tm_raster(n=7, palette ="Reds",
            title="Interpolación kriging \nMapa de varianza \n(en mm)") +tm_shape(P) + tm_dots(size=0.02) +
  tm_legend(legend.outside=TRUE)

CRS object has comment, which is lost in output

Un mapa más fácil de interpretar es el mapa del intervalo de confianza del 95% que se puede generar a partir del objeto de varianza de la siguiente manera (los valores del mapa deben interpretarse como el número de mm por encima y por debajo de la cantidad de lluvia estimada).

r   <- sqrt(raster(dat.krg, layer="var1.var")) * 1.96

CRS object has comment, which is lost in outputDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

r.m <- raster::mask(r, meta2)

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

tm_shape(r.m) + 
  tm_raster(n=7, palette ="Reds",
            title="Interpolación Kriging\n mapa con 95% IC \n(en mm)") +tm_shape(P) + tm_dots(size=0.2) +
  tm_legend(legend.outside=TRUE)

CRS object has comment, which is lost in output