My sixth notebook: Climate data in R
Ivan Lizarazo
13.04.2021
1. Why this notebook?
This is an R Notebook written to illustrate how to obtain, process and visualize climate data in R. It aims to help Geomatica Basica students at Universidad Nacional to get started with R geospatial capabilities. Note that, here, I focus on explaining the technical procedure rather than on interpreting or analyzing the output.
A few tips for writing your own notebook:
- Write and run every chunk step-by-step
- In case of errors, read the message, interpret it and try to fix it
- In case the error persists, write a message on slack (be as precise and descriptive as possible)
- Once your html file is produced, it is time to publish it at RPubs
- After publishing the notebook, review it, fix poor quality wording, and republish it using the update option
In case of trouble, please review the rmarkdown documentation as well as the Leaflet for R documentation.
Do not forget to look for help at community.rstudio.com/.
2. Introduction to climateR
There are several R packages that can be used to get access to climate data. Here, we will use climateR which seeks to simplifiy the steps needed to it. It currently provides access to the following gridded climate sources using a single parameter:
Note that not all the available climate data have global coverage. Many of them cover only the USA.
Let’s start by installing the package:
#install the packages from the console not from here
#remotes::install_github("mikejohnson51/AOI") # suggested!
#remotes::install_github("mikejohnson51/climateR")rm(list=ls())library(AOI)
library(climateR)
library(sf)
library(raster)
library(rasterVis)
library(dplyr)3. Get Climate Data
climateR offers support for \(sf\), \(sfc\), and \(bbox\) objects. Here we will be requesting data for the boundary box of the Boyaca department.
(boyaca <- st_read("../../datos/shapefiles/15_BOYACA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))## Reading layer `MGN_MPIO_POLITICO' from data source `/Users/ivanlizarazo/Documents/ivan/GB/datos/shapefiles/15_BOYACA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 123 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## Geodetic CRS: WGS 843.1 Obtaining meteorological TerraClimate data
The TerraClimate dataset comprises meteorological and water balance variables for 1958-present, available on a monthly timestep. Please go to the [TerraClimate website] (https://climatedataguide.ucar.edu/climate-data) to know more about these dataset. Take notes on the available variables, the way they derive them, and their spatial resolution.
Let’s get the data:
# precipitation
tc_prcp = getTerraClim(boyaca, param = "prcp", startDate = "2019-02-01")tc_tmp <- tc_prcp[[1]]tc_tmp## class : RasterStack
## dimensions : 59, 66, 3894, 1 (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -74.66667, -71.91667, 4.625, 7.083333 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X2019.02
## min values : 0
## max values : 212library(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 here.
# PDSI
tc_palmer = getTerraClim(boyaca, param = "palmer", startDate = "2019-02-01")tc_tmp <- tc_palmer[[1]]tc_tmp## class : RasterStack
## dimensions : 59, 66, 3894, 1 (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -74.66667, -71.91667, 4.625, 7.083333 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X2019.02
## min values : -3.4
## max values : 1.4pal <- 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")3.2 Obtaining TerraClimate normals
Let’s get the Climate Water Deficit data in February (period: 1981- 2010):
# CWD
wat_def = getTerraClimNormals(boyaca, param = "water_deficit", period = "19812010", month=2)wat_def## $`19812010_water_deficit`
## class : RasterStack
## dimensions : 59, 66, 3894, 1 (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -74.66667, -71.91667, 4.625, 7.083333 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X02
## min values : 3.4
## max values : 65.2tc_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")leaflet() %>% addTiles() %>%
addRasterImage(tc_tmp$X02, colors = pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(tc_tmp$X02),
title = "WaterDeficit-February")In your Informe No. 2 analyze the eventual link between topography and climate.
3.2 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(boyaca, startDate = "2020-11-01", endDate = "2020-11-06" )chirps## class : RasterStack
## dimensions : 49, 56, 2744, 6 (nrow, ncol, ncell, nlayers)
## resolution : 0.05, 0.05 (x, y)
## extent : -74.7, -71.9, 4.65, 7.1 (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 : 49, 56, 2744, 6 (nrow, ncol, ncell, nlayers)
## resolution : 0.05, 0.05 (x, y)
## extent : -74.7, -71.9, 4.65, 7.1 (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
## 52.147959 11.902697 9.064869 16.631560 29.383382 14.997813for (penta in capas) {
tmp <- chirps[[penta]]
print(cellStats(tmp, max))
}## [1] 172
## [1] 118
## [1] 100
## [1] 124
## [1] 154
## [1] 218(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")# Print the map
m4. 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]])
}
# 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= 52.15"
## [1] "sd= 29.3"
## [1] "pentad= P2"
## [1] "mean= 11.9"
## [1] "sd= 20.67"
## [1] "pentad= P3"
## [1] "mean= 9.06"
## [1] "sd= 15.15"
## [1] "pentad= P4"
## [1] "mean= 16.63"
## [1] "sd= 16.67"
## [1] "pentad= P5"
## [1] "mean= 29.38"
## [1] "sd= 26.97"
## [1] "pentad= P6"
## [1] "mean= 15"
## [1] "sd= 24.96"5. Reproducibility
Cite as: Lizarazo, I., 2021. A tutorial on getting global TerraClimate and CHIRPS data from R. Available at https://rpubs.com/ials2un/trrclmt.
sessionInfo()## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## 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.4
## [4] rasterVis_0.49 latticeExtra_0.6-29 lattice_0.20-41
## [7] raster_3.4-5 sp_1.4-5 sf_0.9-8
## [10] climateR_0.0.4 AOI_0.2.0.9000
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.2 rnaturalearth_0.1.0 sass_0.3.1
## [4] jsonlite_1.7.2 viridisLite_0.3.0 foreach_1.5.1
## [7] bslib_0.2.4 shiny_1.6.0 assertthat_0.2.1
## [10] highr_0.8 yaml_2.2.1 pillar_1.5.1
## [13] glue_1.4.2 digest_0.6.27 promises_1.2.0.1
## [16] rvest_0.3.6 colorspace_2.0-0 htmltools_0.5.1.1
## [19] httpuv_1.5.5 pkgconfig_2.0.3 purrr_0.3.4
## [22] xtable_1.8-4 scales_1.1.1 jpeg_0.1-8.1
## [25] later_1.1.0.1 tibble_3.1.0 proxy_0.4-25
## [28] leaflet.extras_1.0.0 generics_0.1.0 farver_2.0.3
## [31] ellipsis_0.3.1 hexbin_1.28.2 magrittr_2.0.1
## [34] crayon_1.4.1 mime_0.10 evaluate_0.14
## [37] fansi_0.4.2 doParallel_1.0.16 xml2_1.3.2
## [40] class_7.3-18 tools_4.0.4 USAboundaries_0.3.1
## [43] lifecycle_1.0.0 stringr_1.4.0 munsell_0.5.0
## [46] compiler_4.0.4 jquerylib_0.1.3 e1071_1.7-6
## [49] RNetCDF_2.4-2 rlang_0.4.10 classInt_0.4-3
## [52] units_0.7-1 grid_4.0.4 iterators_1.0.13
## [55] htmlwidgets_1.5.3 crosstalk_1.1.1 base64enc_0.1-3
## [58] rmarkdown_2.7 codetools_0.2-18 DBI_1.1.1
## [61] R6_2.5.0 zoo_1.8-8 knitr_1.31
## [64] rgdal_1.5-23 fastmap_1.1.0 utf8_1.1.4
## [67] KernSmooth_2.23-18 stringi_1.5.3 parallel_4.0.4
## [70] Rcpp_1.0.6 vctrs_0.3.6 png_0.1-7
## [73] tidyselect_1.1.0 xfun_0.21