Este cuaderno incluye las siguientes secciones:
- Descarga de datos de precipitación global de CHIRPS: La descarga se realizó siguiendo las indicaciones del profesor, siguiendo el enlace se descargaron los datos de precipitación correspondientes al global_pentad en el directorio tifs, específicamente los datos de los últimos cinco días del mes de Abril de 2020 del Grupo de Riesgos Climáticos de precipitación infrarroja con datos de estación (CHIRPS) y se descomprimió para obtener un archivo tif.
- Instalar los paquetes y cargar las librerías necesarías: Primero se limpirá la memoria y posterior se instalarán los paquetes necesarios como tidyverse, rgdal, sf, raster, tmap, gstat y sp.Así:
rm(list=ls())
list.of.packages <- c("tidyverse", "rgdal", "sf", "raster", "tmap", "gstat", "sp")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
Ahora, se cargan las librerías:
library(rgdal)
## 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.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/Brian/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/Brian/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)
## -- 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.5
## 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)
- Preprocesamiento de datos de precipitación para el departamento de Boyacá: Inicialmente se lee el archivo descomprimido como se muestra:
precipitacionfile <- raster("C:/Users/Brian/Desktop/Daniela/geomatica/Informe final/chirps-v2.0.2020.04.6.tif")
Este archivo representa la lluvia acumulada en los últimos días de abril de 2020, es decir, del 26-30 de abril. Ahora, es necesrio verificar el contenido del objeto precipitacion:
precipitacion
## class : RasterLayer
## dimensions : 667, 2400, 1600800 (nrow, ncol, ncell)
## resolution : 0.15, 0.15 (x, y)
## extent : -180, 180, -50.05, 50 (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 : -9999, 468.8154 (min, max)
Teniendo en cuenta que este es un conjunto de datos con cobertura global y su CRS es un sistema de coordenadas geográficas, debemos cargar un archivo shape que represente nuestra área de interés, en este caso Boyacá:
(aoi <- shapefile("C:/Users/Brian/Desktop/Daniela/Boyacá/Marco Geoestadístico Boyacá/15_BOYACA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))
## class : SpatialPolygonsDataFrame
## features : 123
## extent : -74.66496, -71.94885, 4.655196, 7.055557 (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 : 15, 15001, ALMEIDA, 1537, 25.35271459, 2017, BOYACÃ, 0.204497978002, 0.00206924119091
## max values : 15, 15897, ZETAQUIRÃ, Ordenanza 8 de Diciembre 4 de 1965, 1513.60104233, 2017, BOYACÃ, 2.40391520946, 0.123609862522
Se observa que los dos conjuntos de datos cuentan con sistema de coordenadas geográficas. Ahora, es necesario recortar los datos de precipitación de cultivos por extensión del área de interés “Boyacá”:
precipitacion.crop <- raster::crop(precipitacion, extent(aoi))
Posteriormente es necesario enmascarar la capa raster y verificamos la salida:
precipitacion.mask <- mask(x = precipitacion.crop, mask = aoi)
precipitacion.mask
## class : RasterLayer
## dimensions : 16, 18, 288 (nrow, ncol, ncell)
## resolution : 0.15, 0.15 (x, y)
## extent : -74.7, -72, 4.699999, 7.099999 (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 : 4.699088, 93.04938 (min, max)
Para obtener un mejor mapa, se instala el paquete Leaflet si es necesario y si ya lo tenemos se carga la librería, para usar la función addRasterImage que este paquete nos ofrece:
library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("red","orange", "yellow", "blue", "darkblue"), values(precipitacion.mask),
na.color = "transparent")
leaflet() %>% addTiles() %>%
addRasterImage(precipitacion.mask, colors = pal, opacity = 0.6) %>%
addLegend(pal = pal, values = values(precipitacion.mask),
title = "CHIRPS Precipitaciones en Boyacá del 26/04 al 30/04 en 2020 [mm]")
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
Teniendo en cuenta que los datos que tenemos están en datum WGS84, los cuales están en coordenadas geográficas, no se puede hacer interpolación o no se recomienda. Para hacer la interpolación debemos hacer la conversión a coordenadas planas (más recomendable). Adicional a la transformación de coordenadas o reproyección, se tienen que colocar los datos en el formato apropiado para que las librerías de R puedan trabajar, las librerías que hay trabajan bien con el tipo de estructura SpatialPointsDataFrame, en otras palabras, es la conversión de raster a puntos. Lo anterior se puede ralizar utilizando la función * rasterToPoints* de la librería raster:
precipitacion.points <- rasterToPoints(precipitacion.mask, spatial = TRUE)
precipitacion.points
## class : SpatialPointsDataFrame
## features : 83
## extent : -74.625, -72.075, 4.774999, 7.024999 (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 : 4.69908825556437
## max values : 93.0493799845378
names(precipitacion.points) <- "lluvia"
precipitacion.points
## class : SpatialPointsDataFrame
## features : 83
## extent : -74.625, -72.075, 4.774999, 7.024999 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 1
## names : lluvia
## min values : 4.69908825556437
## max values : 93.0493799845378
str(precipitacion.points)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 83 obs. of 1 variable:
## .. ..$ lluvia: num [1:83] 5.21 4.7 5.44 7.71 6.49 ...
## ..@ coords.nrs : num(0)
## ..@ coords : num [1:83, 1:2] -72.1 -72.4 -72.2 -72.1 -72.4 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "x" "y"
## ..@ bbox : num [1:2, 1:2] -74.62 4.77 -72.07 7.02
## .. ..- 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, tracemos los puntos:
plot(precipitacion.mask, main= "CHIRPS Precipitaciones en Boyacá del 26/04 al 30/04 en 2020 [mm]")
plot(aoi, add=TRUE)
points(precipitacion.points$x, precipitacion.points$y, col = "red", cex = .6)
Verifiquemos cuál es el objeto precipitacion.points:
precipitacion.points
## class : SpatialPointsDataFrame
## features : 83
## extent : -74.625, -72.075, 4.774999, 7.024999 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 1
## names : lluvia
## min values : 4.69908825556437
## max values : 93.0493799845378
Preparación adicional: En caso de que se pierda la conexión, podemos leer nuevamente el shapefile con el área de interés:
(aoi <- shapefile("C:/Users/Brian/Desktop/Daniela/Boyacá/Marco Geoestadístico Boyacá/15_BOYACA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))
## class : SpatialPolygonsDataFrame
## features : 123
## extent : -74.66496, -71.94885, 4.655196, 7.055557 (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 : 15, 15001, ALMEIDA, 1537, 25.35271459, 2017, BOYACÃ, 0.204497978002, 0.00206924119091
## max values : 15, 15897, ZETAQUIRÃ, Ordenanza 8 de Diciembre 4 de 1965, 1513.60104233, 2017, BOYACÃ, 2.40391520946, 0.123609862522
Necesitamos convertirlo a una característica espacial:
boyaca_sf <- sf::st_as_sf(aoi)
Ahora, disolver los límites internos:
(border_sf <-
boyaca_sf %>%
summarise(area = sum(MPIO_NAREA)))
## Simple feature collection with 1 feature and 1 field
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## CRS: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## area geometry
## 1 23138.05 POLYGON ((-73.03051 4.98230...
Convertir la característica espacial en un dataframe espacial:
(border <- as(border_sf, 'Spatial'))
## class : SpatialPolygonsDataFrame
## features : 1
## extent : -74.66496, -71.94885, 4.655196, 7.055557 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
## variables : 1
## names : area
## value : 23138.04812969
Otra conversión:
(boyaca.sf <- st_as_sf(aoi) %>% mutate(MUNIC = MPIO_CNMBR, CODIGO = MPIO_CCDGO) %>% select(MUNIC, CODIGO))
## Simple feature collection with 123 features and 2 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## CRS: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## First 10 features:
## MUNIC CODIGO geometry
## 1 TUNJA 15001 POLYGON ((-73.34014 5.58308...
## 2 ALMEIDA 15022 POLYGON ((-73.36793 5.01349...
## 3 AQUITANIA 15047 POLYGON ((-72.76242 5.63856...
## 4 ARCABUCO 15051 POLYGON ((-73.50487 5.84347...
## 5 BELÉN 15087 POLYGON ((-72.91692 6.08612...
## 6 BERBEO 15090 POLYGON ((-73.0677 5.27048,...
## 7 BETÉITIVA 15092 POLYGON ((-72.81796 5.97422...
## 8 BOAVITA 15097 POLYGON ((-72.64907 6.43640...
## 9 BOYACÃ\201 15104 POLYGON ((-73.34806 5.47411...
## 10 BRICEÑO 15106 POLYGON ((-73.89118 5.73749...
Ahora, intentaremos una st_intersection. Esto se hace usando la biblioteca sf:
# conversión
p.sf <- st_as_sf(precipitacion.points)
## intersectar polígonos con puntos, manteniendo la información de ambos
(precipitacion.sf = st_intersection(boyaca.sf, p.sf))
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
## Simple feature collection with 83 features and 3 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -74.625 ymin: 4.774999 xmax: -72.075 ymax: 7.024999
## CRS: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## First 10 features:
## MUNIC CODIGO lluvia geometry
## 122 CUBARÃ\201 15223 5.211703 POINT (-72.075 7.024999)
## 122.1 CUBARÃ\201 15223 4.699088 POINT (-72.375 6.874999)
## 122.2 CUBARÃ\201 15223 5.437004 POINT (-72.225 6.874999)
## 122.3 CUBARÃ\201 15223 7.714959 POINT (-72.075 6.874999)
## 121 CHISCAS 15180 6.492599 POINT (-72.375 6.724999)
## 123 GÃœICÃ\201N 15332 6.910823 POINT (-72.225 6.724999)
## 121.1 CHISCAS 15180 7.016574 POINT (-72.525 6.574999)
## 121.2 CHISCAS 15180 7.481505 POINT (-72.375 6.574999)
## 123.1 GÃœICÃ\201N 15332 7.693087 POINT (-72.225 6.574999)
## 112 TIPACOQUE 15810 10.633375 POINT (-72.675 6.424999)
Se hacen dos reproyecciones:
p.sf.magna <- st_transform(precipitacion.sf, crs=3116)
boyaca.sf.magna <- st_transform(boyaca.sf, crs=3116)
Una conversión:
(precipitacion2 <- as(p.sf.magna, 'Spatial'))
## class : SpatialPointsDataFrame
## features : 83
## extent : 939365.3, 1221370, 1019824, 1269064 (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 : ALMEIDA, 15022, 4.69908825556437
## max values : VILLA DE LEYVA, 15861, 93.0493799845378
Nuevamente, podemos escribir/guardar el conjunto de datos intermedio (por si acaso):
shapefile(precipitacion2, filename='C:/Users/Brian/Desktop/Daniela/geomatica/Informe final/precipitacion2.shp', overwrite=TRUE)
Si es necesario, lo leemos. De lo contrario, sigamos adelante:
precipitacion2$rainfall <- round(precipitacion2$lluvia, 1)
##precip2$rainfall <- round(precip2$chirps.v2.0.2020.03.6, 1)
Lo que tenemos
precipitacion2
## class : SpatialPointsDataFrame
## features : 83
## extent : 939365.3, 1221370, 1019824, 1269064 (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, rainfall
## min values : ALMEIDA, 15022, 4.69908825556437, 4.7
## max values : VILLA DE LEYVA, 15861, 93.0493799845378, 93
(boyaca2 <- as(boyaca.sf.magna, 'Spatial'))
## class : SpatialPolygonsDataFrame
## features : 123
## extent : 934932.9, 1235253, 1006592, 1272398 (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 : ALMEIDA, 15001
## max values : ZETAQUIRÃ, 15897
Es posible que queramos escribir/guardar el nuevo objeto (por si acaso):
shapefile(boyaca2, filename='C:/Users/Brian/Desktop/Daniela/geomatica/Informe final/boyaca2.shp', overwrite=TRUE)
Asegúrese de que las dos extensiones coincidan:
precipitacion2@bbox <- boyaca2@bbox
Trazar los datos con tmap: Este paquete genera mapas temáticos, los cuales son mapas geográficos en los que se visualizan las distribuciones de datos espaciales. La función tm_shape crea un elemento tmap que especifica un objeto de datos espaciales, al que nos referimos.
tm_shape(boyaca2) + tm_polygons() +
tm_shape(precipitacion2) +
tm_dots(col="rainfall", palette = "RdBu", midpoint = 50,
title="Precipitación muestreada \n (en mm)", size=0.25) +
tm_text("rainfall", just="center", xmod=.6, size = 0.52) +
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 8.1 x64 (build 9600)
##
## 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] RColorBrewer_1.1-2 leaflet_2.0.3 gstat_2.0-6 tmap_3.0
## [5] forcats_0.5.0 stringr_1.4.0 dplyr_0.8.5 purrr_0.3.3
## [9] readr_1.3.1 tidyr_1.0.2 tibble_2.1.3 ggplot2_3.3.0
## [13] tidyverse_1.3.0 sf_0.9-3 raster_3.0-12 rgdal_1.4-8
## [17] sp_1.4-1
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-144 fs_1.3.2 xts_0.12-0 lubridate_1.7.4
## [5] httr_1.4.1 tools_3.6.3 backports_1.1.5 R6_2.4.1
## [9] KernSmooth_2.23-16 lazyeval_0.2.2 rgeos_0.5-2 DBI_1.1.0
## [13] colorspace_1.4-1 withr_2.1.2 tidyselect_1.0.0 curl_4.3
## [17] compiler_3.6.3 leafem_0.1.1 cli_2.0.2 rvest_0.3.5
## [21] xml2_1.2.5 scales_1.1.0 classInt_0.4-3 digest_0.6.25
## [25] foreign_0.8-75 rmarkdown_2.1 base64enc_0.1-3 dichromat_2.0-0
## [29] pkgconfig_2.0.3 htmltools_0.4.0 dbplyr_1.4.2 htmlwidgets_1.5.1
## [33] rlang_0.4.5 readxl_1.3.1 httpcode_0.3.0 rstudioapi_0.11
## [37] FNN_1.1.3 farver_2.0.3 generics_0.0.2 zoo_1.8-7
## [41] jsonlite_1.6.1 crosstalk_1.1.0.1 magrittr_1.5 Rcpp_1.0.3
## [45] munsell_0.5.0 fansi_0.4.1 abind_1.4-5 lifecycle_0.2.0
## [49] stringi_1.4.6 leafsync_0.1.0 yaml_2.2.1 jqr_1.1.0
## [53] tmaptools_3.0 maptools_0.9-9 grid_3.6.3 parallel_3.6.3
## [57] geojsonio_0.9.2 crayon_1.3.4 lattice_0.20-38 stars_0.4-1
## [61] haven_2.2.0 geojson_0.3.2 hms_0.5.3 knitr_1.28
## [65] pillar_1.4.3 spacetime_1.2-3 codetools_0.2-16 crul_0.9.0
## [69] reprex_0.3.0 XML_3.99-0.3 glue_1.3.2 evaluate_0.14
## [73] V8_3.1.0 modelr_0.1.6 vctrs_0.2.4 png_0.1-7
## [77] cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1 xfun_0.12
## [81] lwgeom_0.2-1 broom_0.5.5 e1071_1.7-3 class_7.3-15
## [85] viridisLite_0.3.0 intervals_0.15.2 units_0.6-5
