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)
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)
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:

plot(precip.mask, main= "CHIRPS precipitaciones en Norte de Santander desde 26.04 hasta 30.04 de 2020 [mm]")
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:

plot(precip.mask, main= "Datos CHIRPS de precipitaciones desde el 26 al 30.04.2020 [mm]")
plot(aoi, add=TRUE)

points(precip.points$x, precip.points$y, col = "darkblue", cex = 0.4)

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")

Registered S3 method overwritten by 'geojson':
  method        from     
  print.geojson geojsonsf
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
1 21856.44
                        geometry
1 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
1     CÚCUTA  54001
2      ABREGO  54003
3   ARBOLEDAS  54051
4   BOCHALEMA  54099
5  BUCARASICA  54109
6  CÃ\u0081COTA  54125
7  CÃ\u0081CHIRA  54128
8  CHINÃ\u0081COTA  54172
9   CUCUTILLA  54223
10    DURANIA  54239
                         geometry
1  MULTIPOLYGON (((-72.4778 8....
2  MULTIPOLYGON (((-73.01687 8...
3  MULTIPOLYGON (((-72.73134 7...
4  MULTIPOLYGON (((-72.60265 7...
5  MULTIPOLYGON (((-72.95019 8...
6  MULTIPOLYGON (((-72.62101 7...
7  MULTIPOLYGON (((-73.04222 7...
8  MULTIPOLYGON (((-72.58771 7...
9  MULTIPOLYGON (((-72.79776 7...
10 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
36     EL CARMEN  54245
36.1   EL CARMEN  54245
36.2   EL CARMEN  54245
35   CONVENCIÓN  54206
36.3   EL CARMEN  54245
36.4   EL CARMEN  54245
36.5   EL CARMEN  54245
36.6   EL CARMEN  54245
35.1 CONVENCIÓN  54206
35.2 CONVENCIÓN  54206
        Lluvia
36   21.525385
36.1 13.994514
36.2 16.871891
35   16.749319
36.3  8.823729
36.4  9.523139
36.5 10.936153
36.6 12.797283
35.1 10.394583
35.2 10.106214
                     geometry
36   POINT (-73.025 9.274999)
36.1 POINT (-73.125 9.224999)
36.2 POINT (-73.075 9.224999)
35   POINT (-73.025 9.224999)
36.3 POINT (-73.225 9.174999)
36.4 POINT (-73.175 9.174999)
36.5 POINT (-73.125 9.174999)
36.6 POINT (-73.075 9.174999)
35.1 POINT (-73.025 9.174999)
35.2 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.5) + tm_legend(legend.outside = TRUE) + tm_shape(precip2)

Error: tm_shape element missing