Datos climáticos de Casanare en R

1. Introducción

CHIRPS que incorpora imágenes de satélite de 5 km × 5 km de resolución con datos de estaciones in situ para crear series temporales de precipitaciones cuadriculadas para el análisis de tendencias, eventos severos y seguimiento de sequías estacionales (Tymvios et al., 2016).

2. Librerias

# remotes::install_github("mikejohnson51/AOI")
# remotes::install_github("mikejohnson51/climateR")

rm(list=ls())

library(AOI)
library(climateR)
library(sf)

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

library(raster)

Loading required package: sp

library(rasterVis)

Loading required package: lattice

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

library(leaflet)

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(RColorBrewer)

3. Datos climáticos

(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
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
First 10 features:
   DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR                                MPIO_CRSLC
1          85        001      YOPAL     Ordenanza 0038 del 8 de Julio de 1942
2          85        010    AGUAZUL                                      1950
3          85        015    CHÁMEZA Ordenanza 0021 del 3 de Diciembre de 1959
4          85        136  LA SALINA                                      1780
5          85        139       MANÍ                                      1953
6          85        162  MONTERREY                                      1960
7          85        225    NUNCHÍA                                      1748
8          85        230     OROCUÉ                                      1845
9          85        263       PORE                                      1799
10         85        279    RECETOR                                      1798
   MPIO_NAREA MPIO_CCNCT MPIO_NANO DPTO_CNMBR SHAPE_AREA SHAPE_LEN ORIG_FID
1   2482.9386      85001      2020   CASANARE 0.20233677 3.1592326       55
2   1442.5073      85010      2020   CASANARE 0.11756477 2.2242959       56
3    313.5325      85015      2020   CASANARE 0.02556381 1.0678483       57
4    199.8943      85136      2020   CASANARE 0.01631711 0.8413073       58
5   3754.5485      85139      2020   CASANARE 0.30568793 4.1536144       59
6    779.1725      85162      2020   CASANARE 0.06349198 1.5420451       60
7   1101.8132      85225      2020   CASANARE 0.08981176 2.0246932       61
8   4753.4841      85230      2020   CASANARE 0.38683367 4.1890475       62
9    780.4712      85263      2020   CASANARE 0.06361918 1.7216899       63
10   181.2268      85279      2020   CASANARE 0.01477602 0.6745774       64
                         geometry
1  POLYGON ((-72.39513 5.56853...
2  POLYGON ((-72.56545 5.36972...
3  POLYGON ((-72.81017 5.36659...
4  POLYGON ((-72.33885 6.34470...
5  POLYGON ((-72.34155 5.06495...
6  POLYGON ((-72.89989 5.03465...
7  POLYGON ((-72.19558 5.71924...
8  POLYGON ((-71.50965 5.21730...
9  POLYGON ((-72.04587 5.83819...
10 POLYGON ((-72.80501 5.38092...

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

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_CCNCT 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...

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...

3.1 Obtención de datos meteorológicos de TerraClimate

tc_prcp = getTerraClim(casanare, param = "prcp", startDate = "2019-02-01")

Spherical geometry (s2) switched off
Spherical geometry (s2) switched on

tc_tmp <- tc_prcp[[1]]

tc_tmp

class      : RasterStack 
dimensions : 51, 78, 3978, 1  (nrow, ncol, ncell, nlayers)
resolution : 0.04166667, 0.04166667  (x, y)
extent     : -73.08333, -69.83333, 4.25, 6.375  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      : X2019.02 
min values :      0.2 
max values :    161.1 

pal <- colorNumeric(c("red", "orange", "#fcc000","yellow", "cyan", "blue", "#3240cd"), values(tc_tmp$X2019.02), 
                    na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(tc_tmp$X2019.02 , colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(tc_tmp$X2019.02),
    title = "Rainfall-Feb.2019 [mm]")

Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded ellps WGS 84 in Proj4 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_defs +type=crs
Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum World Geodetic System 1984 in Proj4 definition

writeRaster(tc_tmp, filename="precip", format="GTiff", overwrite=TRUE)

head(param_meta$terraclim)

** Índice de Severidad de la Sequía de Palmer (PDSI) **

tc_palmer = getTerraClim(casanare, param = "palmer", startDate = "2019-02-01")

Spherical geometry (s2) switched off
Spherical geometry (s2) switched on

tc_tmp <- tc_palmer[[1]]
tc_tmp

class      : RasterStack 
dimensions : 51, 78, 3978, 1  (nrow, ncol, ncell, nlayers)
resolution : 0.04166667, 0.04166667  (x, y)
extent     : -73.08333, -69.83333, 4.25, 6.375  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      : X2019.02 
min values :     -3.5 
max values :      1.9 

pal <- colorNumeric(c("#fc7300","orange", "yellow","#9acd32", "green"), values(tc_tmp$X2019.02), 
                    na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(tc_tmp$X2019.02, colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(tc_tmp$X2019.02),
    title = "PDSI-Feb.2019")

Warning in colors(.) :
  Some values were outside the color scale and will be treated as NA

3.2 Obtención de las normales de TerraClimate

Datos del déficit hídrico climático en febrero (período: 1981- 2010):

# CWD
wat_def = getTerraClimNormals(casanare, param = "water_deficit", period = "19812010", month=2)

Spherical geometry (s2) switched off
Spherical geometry (s2) switched on

wat_def

$terraclim_19812010_water_deficit
class      : RasterStack 
dimensions : 51, 78, 3978, 1  (nrow, ncol, ncell, nlayers)
resolution : 0.04166667, 0.04166667  (x, y)
extent     : -73.08333, -69.83333, 4.25, 6.375  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      :  X02 
min values :  4.8 
max values : 83.4 

tc_tmp <- wat_def[[1]]

pal <- colorNumeric(c("green", "#9acd32","yellow", "orange", 
                    "#fc7300"), values(tc_tmp$X02), 
                    na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(tc_tmp$X02, colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(tc_tmp$X02),
    title = "WaterDeficit-February")

Warning in colors(.) :
  Some values were outside the color scale and will be treated as NA

En su Informe Nº 2 analice el eventual vínculo entre la topografía y el clima.

3.3 Obtención de los datos de CHIRPS

chirps = getCHIRPS(casanare, startDate = "2020-11-01", endDate = "2020-11-06" )

Spherical geometry (s2) switched off
Spherical geometry (s2) switched on

chirps

class      : RasterStack 
dimensions : 42, 66, 2772, 6  (nrow, ncol, ncell, nlayers)
resolution : 0.05, 0.05  (x, y)
extent     : -73.10001, -69.80001, 4.25, 6.35  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      : prcp_20201101, prcp_20201102, prcp_20201103, prcp_20201104, prcp_20201105, prcp_20201106 
min values :             0,             0,             0,             0,             0,             0 
max values :           255,           255,           255,           255,           255,           255 

names(chirps) <- c("P1","P2","P3", "P4", "P5", "P6")

chirps

class      : RasterStack 
dimensions : 42, 66, 2772, 6  (nrow, ncol, ncell, nlayers)
resolution : 0.05, 0.05  (x, y)
extent     : -73.10001, -69.80001, 4.25, 6.35  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      :  P1,  P2,  P3,  P4,  P5,  P6 
min values :   0,   0,   0,   0,   0,   0 
max values : 255, 255, 255, 255, 255, 255 

dimensions: el “tamaño” del archivo en píxeles nrow, ncol: el número de filas y columnas de los datos (imagine una hoja de cálculo o una matriz). ncells: el número total de píxeles o celdas que componen el raster. resolution: el tamaño de cada píxel (en grados decimales en este caso). Recordemos que 1 grado decimal representa aproximadamente 111,11 km en el Ecuador. extent: la extensión espacial de la trama. Este valor estará en las mismas unidades de coordenadas que el sistema de referencia de coordenadas del raster. crs: la cadena del sistema de referencia de coordenadas del ráster. Este ráster está en coordenadas geográficas con un datum de WGS 84. Obsérvese también que cada capa representa las precipitaciones acumuladas en cinco días.

(capas <- names(chirps))

[1] "P1" "P2" "P3" "P4" "P5" "P6"

pentads <-c("P1", "P2", "P3", "P4", "P5", "P6")

cellStats(chirps, mean)

       P1        P2        P3        P4        P5        P6 
44.117605 20.787879 17.773449 24.879509  7.705988  7.511544 

for (penta in capas) {
  tmp <- chirps[[penta]]
  print(cellStats(tmp, max))
}

[1] 171
[1] 115
[1] 100
[1] 122
[1] 77
[1] 218

(valores <- seq(from=1,to=230,by=4))

 [1]   1   5   9  13  17  21  25  29  33  37  41  45  49  53  57  61  65  69  73
[20]  77  81  85  89  93  97 101 105 109 113 117 121 125 129 133 137 141 145 149
[39] 153 157 161 165 169 173 177 181 185 189 193 197 201 205 209 213 217 221 225
[58] 229

pal <- colorNumeric(c("red", "orange", "#fcc000","yellow", 
                        "cyan", "blue", "#3240cd"), 
                         valores, na.color = "transparent")

# Create leaflet widget --------------------------------------------------------
m <- leaflet() %>%
  addTiles(group = "OSM (default)") 
# Add multiple layers with a loop ----------------------------------------------
 for (penta in capas) {
  tmp <- chirps[[penta]]
  m <- m %>% 
  addRasterImage(tmp, colors = pal, opacity = 0.8, group= penta) 
 }

#
# Additional leaflet options 
m <- m %>%
  # Add layers controls
  addLayersControl(
    baseGroups = pentads,
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  # Add common legend
  addLegend(pal = pal, values = valores,
   title = "CHIRPS - Nov.2020")

m

Estadistica basica de los datos climaticos

# histograma
par(mfrow=c(2,3))
for (penta in capas) {
  hist(chirps[[penta]])
}

# works for large files
for (penta in capas) {
   tmp <-  chirps[[penta]]
   media= cellStats(tmp, 'mean')
   desv= cellStats(tmp, 'sd')
   print(paste("pentad=", penta))
   print(paste("mean=", round(media,digits = 2)))
   print(paste("sd=", round(desv,digits=2)))
}

[1] "pentad= P1"
[1] "mean= 44.12"
[1] "sd= 27.74"
[1] "pentad= P2"
[1] "mean= 20.79"
[1] "sd= 22.61"
[1] "pentad= P3"
[1] "mean= 17.77"
[1] "sd= 17.07"
[1] "pentad= P4"
[1] "mean= 24.88"
[1] "sd= 18.05"
[1] "pentad= P5"
[1] "mean= 7.71"
[1] "sd= 12.23"
[1] "pentad= P6"
[1] "mean= 7.51"
[1] "sd= 17.48"

writeRaster(tc_tmp, filename="precip", format="GTiff", overwrite=TRUE)