Datos CHIRPS

Ana María Montaño Hernández

1. Descargar datos de precipitación global.

El Grupo de Riesgos Climáticos de precipitación infrarroja con datos de estación (CHIRPS) es un conjunto de datos de precipitación cuasi global de más de 35 años. Con un alcance de 50 ° S-50 ° N (y todas las longitudes) y desde 1981 hasta el presente, CHIRPS incorpora nuestra climatología interna, CHPclim, imágenes satelitales con resolución de 0.05 ° y datos de estaciones in situ para crear series de tiempo de lluvia cuadriculadas para análisis de tendencias y monitoreo estacional de sequías.

library(rgdal)

rgdal: version: 1.4-8, (SVN revision 845)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
 Path to GDAL shared files: C:/Users/user/Documents/R/win-library/3.6/rgdal/gdal
 GDAL binary built with GEOS: TRUE 
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: C:/Users/user/Documents/R/win-library/3.6/rgdal/proj
 Linking to sp version: 1.4-1 

library(raster)
library(sf)

Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

library(tidyverse)

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages --------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.3.0     v purrr   0.3.3
v tibble  2.1.3     v dplyr   0.8.4
v tidyr   1.0.2     v stringr 1.4.0
v readr   1.3.1     v 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()

library(tmap)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(gstat)
library(sp)

2. Pre-procesamiento de los datos CHIRPS:

Primero hay que leer el archivo CHIRPS sin comprimir:

precip <- raster("C:/Users/user/Documents/Intro_to_R/chirps-v2.0.2020.04.6.tif")

Este archivo representa la lluvia acumulada de los últimos días de abril de 2020 (es decir, la precipitación del 26 al 30 de abril)

Observemos qué hay en el objeto precip:

precip

class      : RasterLayer 
dimensions : 2000, 7200, 14400000  (nrow, ncol, ncell)
resolution : 0.05, 0.05  (x, y)
extent     : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : C:/Users/user/Documents/Intro_to_R/chirps-v2.0.2020.04.6.tif 
names      : chirps.v2.0.2020.04.6 

Tenga en cuenta que este es un conjunto de datos con cobertura global. El CRS es el sistema de coordenadas geográficas. Ahora, se va a cargar un archivo shapefile que represente nuestra área de interés:

(aoi <- shapefile("C:/Users/user/Documents/Norte de Santander/54_NORTE DE SANTANDER/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 40 
extent      : -73.63379, -72.04761, 6.872201, 9.290847  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 9
names       : DPTO_CCDGO, MPIO_CCDGO,        MPIO_CNMBR,                          MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO,         DPTO_CNMBR,     Shape_Leng,       Shape_Area 
min values  :         54,      54001,            ABREGO,                                1555,   44.94540102,      2017, NORTE DE SANTANDER, 0.455016349208, 0.00368621836307 
max values  :         54,      54874, VILLA DEL ROSARIO, Ordenanza 9 de Noviembre 25 de 1948, 2680.07859017,      2017, NORTE DE SANTANDER,  4.08629755008,   0.220098910874 

Hay que tener en cuenta y asegurarse que los sistemas de referencia de coordenadas para ambos conjuntos de datos sean los mismos. Ahora, se van a Recortar los datos de precipitación:

#Datos de precipitación de cultivos por extensión, del área de interés.
precip.crop <- raster::crop(precip, extent(aoi))

Recortar la capa Raster

precip.mask <- mask(x = precip.crop, mask = aoi)

Chequear la salida:

precip.mask

class      : RasterLayer 
dimensions : 49, 32, 1568  (nrow, ncol, ncell)
resolution : 0.05, 0.05  (x, y)
extent     : -73.65, -72.05, 6.849999, 9.299999  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : chirps.v2.0.2020.04.6 
values     : 2.060457, 38.14587  (min, max)

Luego, se plotea la capa Raster recortada:

png("C:/Users/user/Documents/Intro_to_R/precipitaciones1.png")
plot(precip.mask, main= "CHIRPS precipitaciones en Norte de Santander desde 26.04 hasta 30.04 de 2020 [mm]", cex.main= 0.8 )
plot(aoi, add=TRUE)

También, se puede realizar un mapa mejor usando la librería Leaflet y la función addRasterImage

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

leaflet() %>% addTiles() %>%
  addRasterImage(precip.mask, colors = pal, opacity = 0.6) %>%
  addLegend(pal = pal, values = values(precip.mask),
            title = "CHIRPS precipitaciones en Norte de Santander desde 26.04 hasta 30.04 de 2020 [mm]")

ning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Inf

La conversión de Raster en puntos se puede realizar utilizando la función rasterToPoints de la librería raster:

precip.points <- rasterToPoints(precip.mask, spatial = TRUE)

precip.points

class       : SpatialPointsDataFrame 
features    : 721 
extent      : -73.575, -72.075, 6.924999, 9.274999  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       : chirps.v2.0.2020.04.6 
min values  :      2.06045699119568 
max values  :      38.1458702087402 

names(precip.points) <- "Lluvia"

precip.points

class       : SpatialPointsDataFrame 
features    : 721 
extent      : -73.575, -72.075, 6.924999, 9.274999  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       :           Lluvia 
min values  : 2.06045699119568 
max values  : 38.1458702087402 

str(precip.points)

Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 721 obs. of  1 variable:
  .. ..$ Lluvia: num [1:721] 21.53 13.99 16.87 16.75 8.82 ...
  ..@ coords.nrs : num(0) 
  ..@ coords     : num [1:721, 1:2] -73 -73.1 -73.1 -73 -73.2 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "x" "y"
  ..@ bbox       : num [1:2, 1:2] -73.57 6.92 -72.07 9.27
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

Ahora, se plotean los puntos:

png("C:/Users/user/Documents/Intro_to_R/puntos.png")
plot(precip.mask, main= "Datos CHIRPS de precipitaciones desde el 26 al 30.04.2020 [mm]", cex.main= 0.9)
plot(aoi, add=TRUE)
points(precip.points$x, precip.points$y, col = "darkblue", cex = 0.5)

Escribamos los objetos espaciales de puntos en un archivo de disco. Esto para evitar que algo salga mal. Se va a utilizar la librería rgdal

geojsonio::geojson_write(precip.points, file = "C:/Users/user/Documents/Intro_to_R/ppoints.geojson")

Success! File is at C:/Users/user/Documents/Intro_to_R/ppoints.geojson

<geojson-file>
  Path:       C:/Users/user/Documents/Intro_to_R/ppoints.geojson
  From class: SpatialPointsDataFrame

Se van a leer los puntos de precipitación. Hay que consultar la documentación de geojsonio para comprender qué parámetros deben pasarse a la función geojson_read.

precip.points <- geojsonio::geojson_read("C:/Users/user/Documents/Intro_to_R/ppoints.geojson", what= "sp")

Verifiquemos qué hay en el objeto precip.points:

precip.points

class       : SpatialPointsDataFrame 
features    : 721 
extent      : -73.575, -72.075, 6.924999, 9.274999  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       :           Lluvia 
min values  : 2.06045699119568 
max values  : 38.1458702087402 

Pero, ¿Cuál es el punto de convertir los datos de precipitación ráster en datos de precipitación de puntos? Bueno, es porque hay pocas estaciones meteorológicas de la Organización Meteorológica Mundial (OMM) en nuestro país. Por lo tanto, CHIRPS puede ser una opción para obtener los datos de precipitación puntuales necesarios para obtener una superficie de precipitación continua.

Preparación adicional

En caso de que se pierda la conexión, podemos leer nuevamente el archivo de forma con el área de interés:

(aoi <- shapefile("C:/Users/user/Documents/Norte de Santander/54_NORTE DE SANTANDER/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 40 
extent      : -73.63379, -72.04761, 6.872201, 9.290847  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 9
names       : DPTO_CCDGO, MPIO_CCDGO,        MPIO_CNMBR,                          MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO,         DPTO_CNMBR,     Shape_Leng,       Shape_Area 
min values  :         54,      54001,            ABREGO,                                1555,   44.94540102,      2017, NORTE DE SANTANDER, 0.455016349208, 0.00368621836307 
max values  :         54,      54874, VILLA DEL ROSARIO, Ordenanza 9 de Noviembre 25 de 1948, 2680.07859017,      2017, NORTE DE SANTANDER,  4.08629755008,   0.220098910874 

Ahora, es necesario convertirlo a una característica espacial.

NortedeSantander_sf <- sf::st_as_sf(aoi)

Disolver los límites internos:

(border_sf <- 
   NortedeSantander_sf %>%
   summarise(area = sum(MPIO_NAREA)))

Simple feature collection with 1 feature and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -73.63379 ymin: 6.872201 xmax: -72.04761 ymax: 9.290847
CRS:            +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
      area                       geometry
1 21856.44 POLYGON ((-72.4556 7.553288...

Convertir de característica espacial a dataframe espacial:

(border <- as(border_sf, "Spatial"))

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : -73.63379, -72.04761, 6.872201, 9.290847  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 
variables   : 1
names       :           area 
value       : 21856.43898003 

Otra conversión:

(NortedeSantander_sf <- st_as_sf(aoi) %>% mutate(MUNIC = MPIO_CNMBR, CODIGO = MPIO_CCDGO) %>% select(MUNIC, CODIGO))

Simple feature collection with 40 features and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -73.63379 ymin: 6.872201 xmax: -72.04761 ymax: 9.290847
CRS:            +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
First 10 features:
        MUNIC CODIGO                       geometry
1     CÚCUTA  54001 MULTIPOLYGON (((-72.4778 8....
2      ABREGO  54003 MULTIPOLYGON (((-73.01687 8...
3   ARBOLEDAS  54051 MULTIPOLYGON (((-72.73134 7...
4   BOCHALEMA  54099 MULTIPOLYGON (((-72.60265 7...
5  BUCARASICA  54109 MULTIPOLYGON (((-72.95019 8...
6  CÃ\u0081COTA  54125 MULTIPOLYGON (((-72.62101 7...
7  CÃ\u0081CHIRA  54128 MULTIPOLYGON (((-73.04222 7...
8  CHINÃ\u0081COTA  54172 MULTIPOLYGON (((-72.58771 7...
9   CUCUTILLA  54223 MULTIPOLYGON (((-72.79776 7...
10    DURANIA  54239 MULTIPOLYGON (((-72.63625 7...

Ahora se usará la función st_intersection de la libreria sf

p.sf <- st_as_sf(precip.points)

# Intersecta polígonos con puntos, manteniendo la información de ambos
(precip.sf = st_intersection(NortedeSantander_sf, p.sf))

although coordinates are longitude/latitude, st_intersection assumes that they are planar
attribute variables are assumed to be spatially constant throughout all geometries

Simple feature collection with 721 features and 3 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -73.575 ymin: 6.924999 xmax: -72.075 ymax: 9.274999
CRS:            +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
First 10 features:
           MUNIC CODIGO    Lluvia                 geometry
36     EL CARMEN  54245 21.525385 POINT (-73.025 9.274999)
36.1   EL CARMEN  54245 13.994514 POINT (-73.125 9.224999)
36.2   EL CARMEN  54245 16.871891 POINT (-73.075 9.224999)
35   CONVENCIÓN  54206 16.749319 POINT (-73.025 9.224999)
36.3   EL CARMEN  54245  8.823729 POINT (-73.225 9.174999)
36.4   EL CARMEN  54245  9.523139 POINT (-73.175 9.174999)
36.5   EL CARMEN  54245 10.936153 POINT (-73.125 9.174999)
36.6   EL CARMEN  54245 12.797283 POINT (-73.075 9.174999)
35.1 CONVENCIÓN  54206 10.394583 POINT (-73.025 9.174999)
35.2 CONVENCIÓN  54206 10.106214 POINT (-72.975 9.174999)

Dos reproyecciones:

p.sf.magna <- st_transform(precip.sf, crs = 3116)

NortedeSantander.sf.magna <- st_transform(NortedeSantander_sf, crs =3116)

(precip2 <- as(p.sf.magna, "Spatial"))

class       : SpatialPointsDataFrame 
features    : 721 
extent      : 1055435, 1221276, 1257812, 1517605  (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 
variables   : 3
names       :             MUNIC, CODIGO,           Lluvia 
min values  :            ABREGO,  54001, 2.06045699119568 
max values  : VILLA DEL ROSARIO,  54874, 38.1458702087402 

shapefile(precip2, filename = "C:/Users/user/Documents/Intro_to_R/precip2.shp", overwrite = TRUE)

precip2$precipitaciones <- round(precip2$Lluvia, 1)

precip2

class       : SpatialPointsDataFrame 
features    : 721 
extent      : 1055435, 1221276, 1257812, 1517605  (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 
variables   : 4
names       :             MUNIC, CODIGO,           Lluvia, precipitaciones 
min values  :            ABREGO,  54001, 2.06045699119568,             2.1 
max values  : VILLA DEL ROSARIO,  54874, 38.1458702087402,            38.1 

(NortedeSantander2 <- as(NortedeSantander.sf.magna, "Spatial"))

class       : SpatialPolygonsDataFrame 
features    : 40 
extent      : 1048947, 1224322, 1252020, 1519363  (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 
variables   : 2
names       :             MUNIC, CODIGO 
min values  :            ABREGO,  54001 
max values  : VILLA DEL ROSARIO,  54874 

precip2@bbox <- NortedeSantander2@bbox

Plotear los datos con tmpap

tm_shape(NortedeSantander2) + tm_polygons()+ tm_shape(precip2) + tm_dots(col = "precipitaciones", palette = "Paired", midpoint = 3.0,                               title = "Precipitación muestreada (en mm)", size = 0.2)

tm_text("precipitaciones", just = "center", xmod = 0.6, size = 0.4) + tm_legend(legend.outside = TRUE))

Error: inesperado ')' in "tm_text("precipitaciones", just = "center", xmod = 0.6, size = 0.4) + tm_legend(legend.outside = TRUE))"

3. Interpolación de datos de precipitación

3.1 Polígonos Thiessen

Los polígonos de Thiessen (o interpolación de proximidad) se pueden crear utilizando la función dirichlet de spatstat.

th <- as(dirichlet(as.ppp(precip2)), "SpatialPolygons")

crs(th) <- crs(precip2)
crs(NortedeSantander2) <- crs(precip2)

crs(th)

CRS arguments:
 +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 

crs(precip2)

CRS arguments:
 +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 

th.z <- over(th, precip2, fn=mean)
th.spdf <- SpatialPolygonsDataFrame(th, th.z)

th.clp <- raster::intersect(NortedeSantander2, th.spdf)

tm_shape(th.clp) + tm_polygons(col= "Lluvia", palette= "RdYlBu", midpoint= 3.0,
                               title = "Polígonos Thiessen\n Precipitación predicha\n [en mm]") +
tm_legend(legend.outside= TRUE)

3.2 Interpolación usando una ponderación basada en el inverso de la distancia (IDW)

grd  <- as.data.frame(spsample(precip2, "regular", n= 100000))
names(grd)  <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)  <- TRUE
fullgrid(grd) <- TRUE


proj4string(grd)  <- proj4string(precip2)

P.idw <- gstat::idw(Lluvia ~ 1, precip2, newdata= grd, idp = 2.0)

[inverse distance weighted interpolation]

r  <- raster(P.idw)

r

class      : RasterLayer 
dimensions : 390, 256, 99840  (nrow, ncol, ncell)
resolution : 684.7279, 684.7279  (x, y)
extent     : 1049175, 1224466, 1252357, 1519401  (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     : 2.085806, 37.93612  (min, max)

NortedeSantander2

class       : SpatialPolygonsDataFrame 
features    : 40 
extent      : 1048947, 1224322, 1252020, 1519363  (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 
variables   : 2
names       :             MUNIC, CODIGO 
min values  :            ABREGO,  54001 
max values  : VILLA DEL ROSARIO,  54874 

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

r.m

class      : RasterLayer 
dimensions : 390, 256, 99840  (nrow, ncol, ncell)
resolution : 684.7279, 684.7279  (x, y)
extent     : 1049175, 1224466, 1252357, 1519401  (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     : 2.085806, 37.93612  (min, max)

Ahora, se va a plotear:

tm_shape(r.m) + tm_raster(n=8, palette = "RdBu", auto.palette.mapping = FALSE,
                          title = "Distancia Inversa Ponderada\n Precipitación prevista [en mm]") + tm_shape(precip2) + tm_dots(size = 0.05) + 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.

Ahora, se obtendrá una mejor visualización de la superficie interpolada.

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 interpolada con IDW en Norte de Santander \ndel 26.04 al 30.04 del 2020 [mm] " )

ning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Inf

Afinando la interpolación.

Para ajustar la elección del parámetro de potencia, puede realizar una rutina de validación de omisión para medir el error en los valores interpolados.

precip2

class       : SpatialPointsDataFrame 
features    : 721 
extent      : 1048947, 1224322, 1252020, 1519363  (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 
variables   : 4
names       :             MUNIC, CODIGO,           Lluvia, precipitaciones 
min values  :            ABREGO,  54001, 2.06045699119568,             2.1 
max values  : VILLA DEL ROSARIO,  54874, 38.1458702087402,            38.1 

P <- precip2

IDW.out <- vector(length = length(P))

for (i in 1:length(P)) {
  IDW.out[i] <- gstat::idw(Lluvia ~ 1, P[-i,], P[i,], idp= 2.0) $var1.pred
}

Plotear las diferencias

OP <- par(pty = "s", mar=c(4,3,0,0))
plot(IDW.out ~ P$Lluvia, asp=1, xlab= "Observado", ylab ="Predicho", pch= 16, 
     col=rgb(0,0,0,0.5))
abline(lm(IDW.out ~ P$Lluvia), col="red", lw=2, lty=2)

abline(0,1)

par(OP)

El error cuadrático medio (RMSE) se puede calcular desde IDW.out de la siguiente manera:

sqrt( sum((IDW.out - P$Lluvia)^2) / length(P))

[1] 2.393775

Revisar este valor, y arreglarlo porque el error está muy alto. Tal vez revisar el parámetro de potencia de idp = 2.0

Validación cruzada (Cross-validation)

Además de generar una superficie interpolada, se puede crear un mapa de intervalo de confianza del 95% del modelo de interpolación. Aquí se creará un mapa de IC del 95% a partir de una interpolación IDW que utiliza un parámetro de potencia de 2 (idp = 2.0)

3.3 Kriging

3.3.1 Ajustar el modelo de variograma

Primero, es necesario crear un modelo de variograma. Hay que tener 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.

f.1 <- as.formula(precipitaciones ~ X + Y)
P$X <- coordinates(P)[,1]
P$Y <- coordinates(P)[,2]
var.smpl <- variogram(f.1, P, cloud = FALSE, cutoff = 100000, width = 8990)

dat.fit  <- fit.variogram(var.smpl, fit.ranges = TRUE, fit.sills = TRUE, 
                          vgm(psill = 3, model = "Mat", range = 150000, nugget = 0.2))

dat.fit

plot(var.smpl, dat.fit, xlim=c(0,110000))

Generar superficie Kriged Luego, use 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.

f.1 <- as.formula(Lluvia ~ X + Y)
dat.krg <- krige(f.1, P, grd, dat.fit)

[using universal kriging]

r <- raster(dat.krg)
r.m <- raster::mask(r, NortedeSantander2)

r.m

class      : RasterLayer 
dimensions : 390, 256, 99840  (nrow, ncol, ncell)
resolution : 684.7279, 684.7279  (x, y)
extent     : 1049175, 1224466, 1252357, 1519401  (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     : 1.982649, 38.13942  (min, max)

tm_shape(r.m) + 
  tm_raster(n=10, palette="RdBu", auto.palette.mapping=FALSE, 
            title="Universal Kriging\nPrecipitación predicha \n(en mm)") +
  tm_shape(P) + tm_dots(size=0.1) +
  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.

r   <- raster(dat.krg, layer="var1.var")
r.m <- raster::mask(r, NortedeSantander2)

tm_shape(r.m) + 
  tm_raster(n=7, palette ="Reds",
            title="Kriging Interpolation\nMapa de varianza \n(in squared mm)") +tm_shape(P) + tm_dots(size=0.2) +
  tm_legend(legend.outside=TRUE)

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
r.m <- raster::mask(r, NortedeSantander2)

tm_shape(r.m) + 
  tm_raster(n=7, palette ="Reds",
            title="Interpolación kriging\n Mapa de 95% IC\n(en mm)") +tm_shape(P) + tm_dots(size=0.1) +
  tm_legend(legend.outside=TRUE)

sessionInfo()

R version 3.6.3 (2020-02-29)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18363)

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] spatstat_1.64-1     rpart_4.1-15        nlme_3.1-148       
 [4] spatstat.data_1.4-3 RColorBrewer_1.1-2  leaflet_2.0.3      
 [7] gstat_2.0-6         tmap_3.0            forcats_0.5.0      
[10] stringr_1.4.0       dplyr_0.8.4         purrr_0.3.3        
[13] readr_1.3.1         tidyr_1.0.2         tibble_2.1.3       
[16] ggplot2_3.3.0       tidyverse_1.3.0     sf_0.9-3           
[19] rgdal_1.4-8         raster_3.0-12       sp_1.4-1           

loaded via a namespace (and not attached):
  [1] leafem_0.1.1          colorspace_1.4-1      deldir_0.1-25        
  [4] class_7.3-15          rsconnect_0.8.16      base64enc_0.1-3      
  [7] fs_1.3.2              dichromat_2.0-0       httpcode_0.3.0       
 [10] rstudioapi_0.11       farver_2.0.3          fansi_0.4.1          
 [13] lubridate_1.7.4       xml2_1.2.5            splines_3.6.3        
 [16] codetools_0.2-16      knitr_1.28            polyclip_1.10-0      
 [19] jsonlite_1.6.1        tmaptools_3.0         broom_0.5.5          
 [22] dbplyr_1.4.2          png_0.1-7             rgeos_0.5-2          
 [25] shiny_1.4.0           compiler_3.6.3        httr_1.4.1           
 [28] backports_1.1.5       Matrix_1.2-18         assertthat_0.2.1     
 [31] fastmap_1.0.1         lazyeval_0.2.2        cli_2.0.2            
 [34] later_1.0.0           htmltools_0.4.0       tools_3.6.3          
 [37] gtable_0.3.0          glue_1.3.1            geojson_0.3.2        
 [40] V8_3.1.0              Rcpp_1.0.3            cellranger_1.1.0     
 [43] vctrs_0.2.3           crul_0.9.0            leafsync_0.1.0       
 [46] crosstalk_1.0.0       lwgeom_0.2-1          xfun_0.12            
 [49] rvest_0.3.5           mime_0.9              lifecycle_0.2.0      
 [52] goftest_1.2-2         XML_3.99-0.3          jqr_1.1.0            
 [55] zoo_1.8-7             scales_1.1.0          spatstat.utils_1.17-0
 [58] hms_0.5.3             promises_1.1.0        parallel_3.6.3       
 [61] yaml_2.2.1            curl_4.3              stringi_1.4.6        
 [64] maptools_0.9-9        e1071_1.7-3           intervals_0.15.2     
 [67] rlang_0.4.5           pkgconfig_2.0.3       evaluate_0.14        
 [70] lattice_0.20-38       tensor_1.5            htmlwidgets_1.5.1    
 [73] tidyselect_1.0.0      magrittr_1.5          R6_2.4.1             
 [76] geojsonio_0.9.2       generics_0.0.2        DBI_1.1.0            
 [79] mgcv_1.8-31           pillar_1.4.3          haven_2.2.0          
 [82] foreign_0.8-75        withr_2.1.2           units_0.6-5          
 [85] stars_0.4-1           xts_0.12-0            abind_1.4-5          
 [88] spacetime_1.2-3       modelr_0.1.6          crayon_1.3.4         
 [91] KernSmooth_2.23-16    rmarkdown_2.1         grid_3.6.3           
 [94] readxl_1.3.1          FNN_1.1.3             reprex_0.3.0         
 [97] digest_0.6.25         classInt_0.4-3        xtable_1.8-4         
[100] httpuv_1.5.2          munsell_0.5.0         viridisLite_0.3.0