Our sixth notebook: Climate data in R
Manuel Reyes & Mauren Daza
1. Introduction
In this R notebook written with [R Markdown] (http://rmarkdown.rstudio.com), you can see the procedure that is done with climatic data to illustrate, visualize and distinguish geospatial capabilities.
1.1 What is climateR?
There are several packages in R that contain climate data, but, in this notebook we use * ClimateR *, which fulfills an all-in-one function, which simplifies all development. In this package only three databases have global coverage.
It began by installing the package and it is suggested to do the following:
#install the packages from github, do it in the console not from here
#remotes::install_github("mikejohnson51/AOI") # suggested!
#remotes::install_github("mikejohnson51/climateR")library(AOI)
library(climateR)
library(sf)
library(raster)
library(rasterVis)
library(dplyr)#2. Get Climate Data climateR offers support for sf, sfc, and bbox objects. Here we will be requesting data for the boundary box of the Sucre department.
Sucre <- read_sf("C:/Users/Mauren Daza/Desktop/Geo/70_SUCRE/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"); select(Sucre, MPIO_CCDGO, MPIO_CNMBR,)2.1 Obtaining meteorological TerraClimate data
The TerraClimate dataset contains meteorological and hydraulic data from 1958 to the present, which are uploaded monthly through [TerraClimate website] (https://climatedataguide.ucar.edu/climate-data) with its different variables and we can see here how these are derived and how their spatial resolution changes.
For this, the data is obtained:
# precipitation
tc_prcp = getTerraClim(Sucre, param = "prcp", startDate = "2019-02-01")## Spherical geometry (s2) switched off
## Spherical geometry (s2) switched on(tc_tmp <- tc_prcp[[1]])## class : RasterStack
## dimensions : 46, 29, 1334, 1 (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -75.70833, -74.5, 8.25, 10.16667 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X2019.02
## min values : 0.1
## max values : 11.6library(leaflet)
library(RColorBrewer)
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]")We can check what variables are provided by TerraClimate by getting the metadata:
head(param_meta$terraclim)Now, let’s get the Palmer Drought Severity Index (PDSI). More info on PSDI [Data PDSI] (https://climatedataguide.ucar.edu/climate-data/palmer-drought-severity-index-pdsi)
# PDSI
tc_palmer = getTerraClim(Sucre, param = "palmer", startDate = "2019-02-01")## Spherical geometry (s2) switched off
## Spherical geometry (s2) switched on(tc_tmp <- tc_palmer[[1]])## class : RasterStack
## dimensions : 46, 29, 1334, 1 (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -75.70833, -74.5, 8.25, 10.16667 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X2019.02
## min values : -4
## max values : -1.6pal <- 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")2.2 Obtaining TerraClimate normals
Let’s get the Climate Water Deficit data in February (period: 1981- 2010):
# CWD
(wat_def = getTerraClimNormals(Sucre, param = "water_deficit", period = "19812010", month=2))## Spherical geometry (s2) switched off
## Spherical geometry (s2) switched on## $terraclim_19812010_water_deficit
## class : RasterStack
## dimensions : 46, 29, 1334, 1 (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -75.70833, -74.5, 8.25, 10.16667 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X02
## min values : 38.1
## max values : 98.7tc_tmp <- wat_def[[1]]Make sure to check what are the units of the CWD variable:
pal <- colorNumeric(c("green", "#9acd32","yellow", "orange", "#fc7300"),
values(tc_tmp$X02),na.color = "transparent")
#Create map
leaflet() %>% addTiles() %>%
addRasterImage(tc_tmp$X02, colors = pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(tc_tmp$X02),
title = "WaterDeficit-February")2.3 Obtaining CHIRPS data
We can also obtain pentad rainfall data provided by Climate Hazards Center at University of California in Santa Barbara. The dataset is known as Climate Hazards Center InfraRed Precipitation with Station data (CHIRPS). Please go to the CHIRPS data website. Check the main characteristics of the data and take relevant notes.
Let’s get the rainfall corresponding to November 2020. In a month, there are six pentads (i.e. periods of five consecutive days).
(chirps = getCHIRPS(Sucre, startDate = "2020-11-01", endDate = "2020-11-06" ))## Spherical geometry (s2) switched off
## Spherical geometry (s2) switched on## class : RasterStack
## dimensions : 38, 25, 950, 6 (nrow, ncol, ncell, nlayers)
## resolution : 0.05, 0.05 (x, y)
## extent : -75.75001, -74.50001, 8.25, 10.15 (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, 255names(chirps) <- c("P1","P2","P3", "P4", "P5", "P6");chirps## class : RasterStack
## dimensions : 38, 25, 950, 6 (nrow, ncol, ncell, nlayers)
## resolution : 0.05, 0.05 (x, y)
## extent : -75.75001, -74.50001, 8.25, 10.15 (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, 255Notice a few things about this object:
-dimensions: the “size” of the file in pixels -nrow, ncol: the number of rows and columns in the data (imagine a spreadsheet or a matrix). -ncells: the total number of pixels or cells that make up the raster. -resolution: the size of each pixel (in decimal degrees in this case). Let’s remind that 1 decimal degree represents about 111.11 km at the Equator. -extent: the spatial extent of the raster. This value will be in the same coordinate units as the coordinate reference system of the raster. -crs: the coordinate reference system string for the raster. This raster is in geographic coordinates with a datum of WGS 84.
Note also that each layer represents rainfall accumulated in five days.
(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
## 16.701053 2.129474 6.537895 2.500000 15.610526 21.070526for (penta in capas) {
tmp <- chirps[[penta]]
print(cellStats(tmp, max))
}## [1] 62
## [1] 21
## [1] 27
## [1] 20
## [1] 76
## [1] 65(valores <- seq(from=1,to=180,by=5))## [1] 1 6 11 16 21 26 31 36 41 46 51 56 61 66 71 76 81 86 91
## [20] 96 101 106 111 116 121 126 131 136 141 146 151 156 161 166 171 176pal <- 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")
m3. Basic statistics of climate data
A quick exploration of the CHIRPS statistics:
# histograma
par(mfrow=c(2,3))
for (penta in capas) {
hist(chirps[[penta]], col= hcl.colors(5, palette= "Temps"))
}
# 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= 16.7"
## [1] "sd= 11.76"
## [1] "pentad= P2"
## [1] "mean= 2.13"
## [1] "sd= 1.48"
## [1] "pentad= P3"
## [1] "mean= 6.54"
## [1] "sd= 5.59"
## [1] "pentad= P4"
## [1] "mean= 2.5"
## [1] "sd= 2.15"
## [1] "pentad= P5"
## [1] "mean= 15.61"
## [1] "sd= 12.29"
## [1] "pentad= P6"
## [1] "mean= 21.07"
## [1] "sd= 13.54"4. Reproducibility
sessionInfo()## R version 4.0.5 (2021-03-31)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Spanish_Mexico.1252 LC_CTYPE=Spanish_Mexico.1252
## [3] LC_MONETARY=Spanish_Mexico.1252 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Mexico.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] RColorBrewer_1.1-2 leaflet_2.0.4.1 dplyr_1.0.7
## [4] rasterVis_0.50.2 latticeExtra_0.6-29 lattice_0.20-41
## [7] terra_1.3-4 raster_3.4-13 sp_1.4-5
## [10] sf_1.0-1 climateR_0.1.0 AOI_0.2.0.9000
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.2 rnaturalearth_0.1.0 sass_0.4.0
## [4] jsonlite_1.7.2 viridisLite_0.4.0 foreach_1.5.1
## [7] bslib_0.2.5.1 shiny_1.6.0 assertthat_0.2.1
## [10] highr_0.9 yaml_2.2.1 pillar_1.6.1
## [13] glue_1.4.2 digest_0.6.27 promises_1.2.0.1
## [16] rvest_1.0.0 colorspace_2.0-2 htmltools_0.5.1.1
## [19] httpuv_1.6.1 pkgconfig_2.0.3 s2_1.0.6
## [22] purrr_0.3.4 xtable_1.8-4 scales_1.1.1
## [25] jpeg_0.1-8.1 later_1.2.0 tibble_3.1.2
## [28] proxy_0.4-26 leaflet.extras_1.0.0 farver_2.1.0
## [31] generics_0.1.0 ellipsis_0.3.2 hexbin_1.28.2
## [34] magrittr_2.0.1 crayon_1.4.1 mime_0.11
## [37] evaluate_0.14 fansi_0.5.0 doParallel_1.0.16
## [40] xml2_1.3.2 class_7.3-18 tools_4.0.5
## [43] USAboundaries_0.3.1 lifecycle_1.0.0 stringr_1.4.0
## [46] munsell_0.5.0 compiler_4.0.5 jquerylib_0.1.4
## [49] e1071_1.7-7 RNetCDF_2.4-2 rlang_0.4.11
## [52] classInt_0.4-3 units_0.7-2 grid_4.0.5
## [55] iterators_1.0.13 htmlwidgets_1.5.3 crosstalk_1.1.1
## [58] base64enc_0.1-3 rmarkdown_2.9 wk_0.4.1
## [61] codetools_0.2-18 DBI_1.1.1 R6_2.5.0
## [64] zoo_1.8-9 knitr_1.33 rgdal_1.5-23
## [67] fastmap_1.1.0 utf8_1.2.1 KernSmooth_2.23-18
## [70] stringi_1.6.2 parallel_4.0.5 Rcpp_1.0.6
## [73] vctrs_0.3.8 png_0.1-7 tidyselect_1.1.1
## [76] xfun_0.24