Exploring the package {climateR}

Exploring the various datasets available publicly and freely to use and analyze with R, {terra} and {sf} with the {climateR} package

Maps

Geocomputation

{geodata}

Author

Aditya Dahiya

Published

January 30, 2025

Exploring the {climateR} package

Exploring {climateR} (by Mike Johnson) with {sf} (Pebesma and Bivand 2023) and {terra} (Hijmans 2024)

Code

# Install {cliamteR} package
# remotes::install_github("mikejohnson51/AOI") # suggested!
# remotes::install_github("mikejohnson51/climateR")

# Data Import and Wrangling Tools
library(sf)                   # Handling simple features in R
library(terra)                # Handling rasters in R
library(tidyterra)            # Rasters with ggplot2
library(tidyverse)            # All things tidy

# Final plot tools
library(scales)               # Nice Scales for ggplot2
library(fontawesome)          # Icons display in ggplot2
library(ggtext)               # Markdown text in ggplot2
library(showtext)             # Display fonts in ggplot2
library(colorspace)           # Lighten and Darken colours

# Package to explore
library(climateR)             # Climate Data
library(AOI)                  # Get area of interest

# Making tables in R
library(gt)                   # Beautiful Tables

bts = 12 # Base Text Size
sysfonts::font_add_google("Asap Condensed", "body_font")
showtext::showtext_auto()
theme_set(
  theme_minimal(
    base_size = bts,
    base_family = "body_font"
  ) +
    theme(
      text = element_text(
        colour = "grey20",
        lineheight = 0.3,
        margin = margin(0,0,0,0, "pt")
      )
    )
)

# Some basic caption stuff
# A base Colour
bg_col <- "white"
seecolor::print_color(bg_col)

# Colour for highlighted text
text_hil <- "grey30"
seecolor::print_color(text_hil)

# Colour for the text
text_col <- "grey20"
seecolor::print_color(text_col)


# Caption stuff for the plot
sysfonts::font_add(
  family = "Font Awesome 6 Brands",
  regular = here::here("docs", "Font Awesome 6 Brands-Regular-400.otf")
)
github <- "&#xf09b"
github_username <- "aditya-dahiya"
xtwitter <- "&#xe61b"
xtwitter_username <- "@adityadahiyaias"
social_caption_1 <- glue::glue("<span style='font-family:\"Font Awesome 6 Brands\";'>{github};</span> <span style='color: {text_hil}'>{github_username}  </span>")
social_caption_2 <- glue::glue("<span style='font-family:\"Font Awesome 6 Brands\";'>{xtwitter};</span> <span style='color: {text_hil}'>{xtwitter_username}</span>")
plot_caption <- paste0(
  "**Data**: {climateR} package  ",
  "  |  **Code:** ", 
  social_caption_1, 
  " |  **Graphics:** ", 
  social_caption_2
  )
rm(github, github_username, xtwitter, 
   xtwitter_username, social_caption_1, 
   social_caption_2)

`getTerraClim()`

The getTerraClim() function in the {climateR} package provides access to TerraClimate data, a high-resolution dataset that offers monthly climate and water balance variables from 1958 onward. These variables include key climate indicators such as precipitation, temperature, soil moisture, wind speed, and evaporation metrics.

Code

terraClim_vars <- tibble::tibble(
  Variable = c("aet", "def", "PDSI", "pet", "ppt", "q", "soil", "srad", 
               "swe", "tmax", "tmin", "vap", "vpd", "ws"),
  Units = c("mm", "mm", "unitless", "mm", "mm", "mm", "mm", "W/m^2", 
            "mm", "degC", "degC", "kPa", "kPa", "m/s"),
  Description = c("Water evaporation amount", 
                  "Potential evaporation minus actual evaporation",
                  "Palmer Drought Severity Index",
                  "Potential evaporation amount",
                  "Precipitation amount",
                  "Runoff amount",
                  "Soil moisture content",
                  "Downwelling shortwave flux in air",
                  "Liquid water content of surface snow",
                  "Maximum air temperature",
                  "Minimum air temperature",
                  "Water vapor partial pressure in air",
                  "Vapor pressure deficit",
                  "Wind speed")
)

terraClim_vars |> 
  gt::gt() |> 
  gtExtras::gt_theme_538()

Table 1: Summary of Climate and Hydrological Variables in getTerraClim()

Variable	Units	Description
aet	mm	Water evaporation amount
def	mm	Potential evaporation minus actual evaporation
PDSI	unitless	Palmer Drought Severity Index
pet	mm	Potential evaporation amount
ppt	mm	Precipitation amount
q	mm	Runoff amount
soil	mm	Soil moisture content
srad	W/m^2	Downwelling shortwave flux in air
swe	mm	Liquid water content of surface snow
tmax	degC	Maximum air temperature
tmin	degC	Minimum air temperature
vap	kPa	Water vapor partial pressure in air
vpd	kPa	Vapor pressure deficit
ws	m/s	Wind speed

British Isles: Monthly Rainfall Data

This code generates a faceted map of monthly rainfall in the British Isles from January 2018 to December 2022. It begins by defining the area of interest using aoi_get() and then retrieves precipitation data via getTerraClim(). The region is refined using ne_countries() and processed with st_union(), st_as_sf(), and st_crop() to match the raster data’s extent. The precipitation raster is cropped and masked using terra::crop() and terra::mask(). Monthly labels are extracted with str_remove() and format(). The plot is created using ggplot() with geom_spatraster() for the raster data and geom_sf() for country borders. Color scaling is handled by paletteer::scale_fill_paletteer_c(), and faceting is applied with facet_wrap(). The final visualization is styled using theme_map(), and the plot is saved with ggsave().

Code

# Try for a smaller area

uk_vec <- aoi_get(country = c("United Kingdom", "Ireland"))

uk_climate_raw <- getTerraClim(
  uk_vec,
  varname = "ppt",
  startDate = "2018-01-01",
  endDate = "2022-12-31"
)

uk_vec_detail <- rnaturalearth::ne_countries(
  country = c("United Kingdom", "Ireland"),
  scale = "medium",
  returnclass = "sf"
) |> 
  select(geometry) |> 
  st_union() |> 
  st_as_sf() |> 
  st_crop(
    st_bbox(
      c(
        ymin = 49.8,
        ymax = 59.5,
        xmin = -10.3,
        xmax = 2
      ),
      crs = st_crs("EPSG:4326")
    )
  ) |> 
  st_transform(st_crs(uk_climate_raw$ppt))

ggplot(uk_vec_detail) +
  geom_sf()

uk_rast <- uk_climate_raw$ppt |> 
  terra::crop(uk_vec_detail) |> 
  terra::mask(uk_vec_detail)


strip_labels <- uk_rast |> 
  names() |> 
  str_remove("ppt_") |> 
  str_remove("_total") |> 
  as_date() |> 
  format("%B\n%Y")

names(strip_labels) <- uk_rast |> 
  names()


g <- ggplot() +
  geom_spatraster(
    data = uk_rast[[c(1:60)]]
  ) +
  geom_sf(
    data = uk_vec_detail,
    fill = "transparent",
    linewidth = 0.25,
    colour = "grey20"
  ) +
  paletteer::scale_fill_paletteer_c(
    "grDevices::Blues 3",
    direction = -1,
    na.value = "transparent",
    limits = c(0, 450),
    oob = scales::squish
  ) +
  facet_wrap(
    ~lyr,
    labeller = labeller(lyr = strip_labels),
    nrow = 5,
    ncol = 12
  ) +
  coord_sf(clip = "off") +
  labs(
    title = "British Isles: Rainfall Pattern",
    fill = "Monthly Precipitation (mm)"
  ) +
  ggthemes::theme_map(
    base_size = 40,
    base_family = "body_font"
  ) +
  theme(
    legend.position = "bottom",
    legend.title.position = "top",
    strip.text = element_text(
      hjust = 0,
      margin = margin(6,0,-10,0,"pt"),
      size = 48,
      lineheight = 0.3
    ),
    strip.background = element_blank(),
    plot.title = element_text(
      margin = margin(0,0,10,0, "pt"),
      hjust = 0.5,
      size = 120
    ),
    panel.background = element_rect(
      fill = "transparent",
      colour = "transparent"
    ),
    legend.key.height = unit(10, "pt"),
    legend.key.width = unit(100, "pt"),
    legend.title = element_text(
      margin = margin(0,0,2,0, "mm"),
      size = 80,
      hjust = 0.5
    ),
    legend.text = element_text(
      margin = margin(2,0,0,0, "mm")
    ),
    plot.margin = margin(10,0,10,0, "pt"),
    legend.margin = margin(0,0,0,0, "pt"),
    legend.box.margin = margin(0,0,0,0, "pt"),
    legend.justification = c(0.5, 0)
  )

ggsave(
  filename = here::here(
    "geocomputation", "images",
    "climateR_package_1.png"
  ),
  plot = g,
  height = 4500,
  width = 6000,
  units = "px",
  bg = "white"
)

Figure 1: This faceted visualization illustrates monthly precipitation trends across the British Isles from January 2018 to December 2022. Each panel represents a different month, showcasing variations in rainfall intensity using a gradient color scale. Data is sourced from TerraClimate and spatially refined for accuracy.

Precipitation data graphic for India.

Code

# Get a nice map of India from Survey of India (simplified to
# save on computing time)
india_vec <- read_sf(
  here::here(
    "data", "india_map",
    "India_Country_Boundary.shp"
  )
) |> 
  # Simplify to save computing time
  st_simplify(dTolerance = 3000) |> 
  # Keep on bigger szed polygons and multiploygons
  filter(!st_is_empty(geometry)) |> 
  arrange(desc(Area)) |> 
  slice(1:2) |> 
  select(-Area) |> 
  st_transform("EPSG:4326")

# Quick check
india_vec |> 
  ggplot() +
  geom_sf()
#--------------------------------------------------------------
# Minor manual correction: Increase bounding box of India
# as the top of north is slightly getting cropped
# Get the current bounding box of India
india_bbox <- st_bbox(india_vec)

# Increase the northern boundary (ymax) by 2 degrees
india_bbox["ymax"] <- india_bbox["ymax"] + 2

# Convert the updated bounding box to an sf object
india_bbox_sf <- st_as_sfc(india_bbox, crs = st_crs(india_vec))

# Crop india_vec to match the expanded bounding box
india_vec <- st_crop(india_vec, india_bbox_sf)

# Remove temporary objects
rm(india_bbox, india_bbox_sf)
#-------------------------------------------------------------

# Download the TerraClimate dataset on PDSI from 
# https://www.climatologylab.org/terraclimate.html
india_rast_raw <- getTerraClim(
  india_vec,
  varname = "ppt",
  startDate = "2018-01-01",
  endDate = "2023-12-01"
)

# Match CRS of both: prefer to change CRS of vector
india_vec2 <- india_vec |> 
  st_transform(crs(india_rast_raw$ppt))

# Crop and Maskt he Raster to show only the data within 
# India's Administrative Borders
india_rast <- india_rast_raw$ppt |> 
  aggregate(fact = 4) |> 
  terra::crop(india_vec2) |> 
  terra::mask(india_vec2)

# Clean Names for the Panel Strip Text: Month and Year
strip_labels <- india_rast |> 
  names() |> 
  str_remove("ppt_") |> 
  str_remove("_total") |> 
  as_date() |> 
  format("%b\n%Y")

names(strip_labels) <- india_rast |> 
  names()

length(strip_labels)

# Temporarily save india_rast to save download time
# Check its range of values
writeRaster(
  india_rast,
  filename = "india_ppt_rast.tif",
  overwrite = TRUE
)

india_rast |>
  values() |>
  range(na.rm = T)

g2 <- ggplot() +
  geom_spatraster(
    data = india_rast
  ) +
  geom_sf(
    data = india_vec2,
    fill = "transparent",
    linewidth = 0.25,
    colour = "grey20"
  ) +
  paletteer::scale_fill_paletteer_c(
    "grDevices::Blues 3",
    direction = -1,
    na.value = "transparent",
    limits = c(0, 700),
    breaks = seq(0, 600, 200),
    oob = scales::squish
  ) +
  facet_wrap(
    ~lyr,
    labeller = labeller(lyr = strip_labels),
    ncol = 12
  ) +
  coord_sf(clip = "off") +
  labs(
    title = "India: Rainfall Patterns (2018 - 2023)",
    subtitle = str_wrap("Monthly rainfall data for India, from TerraClimate dataset, fetched as a Raster using {climateR}", 130),
    fill = "Precipitation (mm)",
    caption = plot_caption
  ) +
  ggthemes::theme_map(
    base_size = 40,
    base_family = "body_font"
  ) +
  theme(
    legend.position = "bottom",
    legend.title.position = "top",
    strip.text = element_text(
      hjust = 0.7,
      margin = margin(3,0,-20,0,"pt"),
      size = 48,
      lineheight = 0.3
    ),
    strip.background = element_blank(),
    plot.title = element_text(
      margin = margin(0,0,10,0, "pt"),
      hjust = 0.5,
      size = 160
    ),
    plot.subtitle = element_text(
      margin = margin(0,0,10,0, "pt"),
      hjust = 0.5,
      size = 70,
      lineheight = 0.3
    ),
    plot.caption = element_textbox(
      hjust = 0.5,
      margin = margin(30,0,0,0, "pt"),
      size = 60
    ),
    panel.background = element_rect(
      fill = "transparent",
      colour = "transparent"
    ),
    legend.key.height = unit(10, "pt"),
    legend.key.width = unit(100, "pt"),
    legend.title = element_text(
      margin = margin(0,0,2,0, "mm"),
      size = 80,
      hjust = 0.5
    ),
    legend.text = element_text(
      margin = margin(2,0,0,0, "mm"),
      size = 60
    ),
    plot.margin = margin(10,5,5,5, "pt"),
    legend.margin =  margin(-20,0,0,0, "pt"),
    legend.box.margin = margin(-20,0,0,0, "pt"),
    legend.justification = c(0.5, 0)
  )

ggsave(
  filename = here::here(
    "geocomputation", "images",
    "climateR_package_3.png"
  ),
  plot = g2,
  height = 3700,
  width = 4000,
  units = "px",
  bg = "white",
  limitsize = FALSE
)

This faceted visualization illustrates monthly precipitation trends across the India from January 2018 to December 2023. Each panel represents a different month, showcasing variations in rainfall intensity using a gradient color scale. Data is sourced from TerraClimate and spatially refined for accuracy.

India: Drought Severity Trends (2010 - 2023)

The Palmer Drought Severity Index (PDSI) is a widely used metric for measuring long-term drought conditions, based on temperature, precipitation, and soil moisture balance. It provides a standardized scale where negative values indicate drought (with -4 or lower signifying extreme drought) and positive values represent wet conditions.

The analysis shown in Figure 2 uses the TerraClimate dataset from the Climatology Lab, a high-resolution global dataset that offers monthly climate and hydrological variables at a ~4 km resolution. The dataset is accessed via the climateR package, which facilitates retrieval of climate data from TerraClimate. The geographic boundaries for India are obtained using Survey of India, and spatial data operations are conducted with sf and terra. The final visualization employs ggplot2 adn tidyterra to create a faceted graph displaying PDSI trends from 2010 to 2024, offering insights into regional drought patterns over time.

Code

# Get a nice map of India from Survey of India (simplified to
# save on computing time)
india_vec <- read_sf(
  here::here(
    "data", "india_map",
    "India_Country_Boundary.shp"
  )
) |> 
  # Simplify to save computing time
  st_simplify(dTolerance = 3000) |> 
  # Keep on bigger szed polygons and multiploygons
  filter(!st_is_empty(geometry)) |> 
  arrange(desc(Area)) |> 
  slice(1:2) |> 
  select(-Area) |> 
  st_transform("EPSG:4326")

# Quick check
india_vec |> 
  ggplot() +
  geom_sf()

# Minor manual correction: Increase bounding box of India
# as the top of north is slightly getting cropped
# Get the current bounding box of India
india_bbox <- st_bbox(india_vec)

# Increase the northern boundary (ymax) by 2 degrees
india_bbox["ymax"] <- india_bbox["ymax"] + 2

# Convert the updated bounding box to an sf object
india_bbox_sf <- st_as_sfc(india_bbox, crs = st_crs(india_vec))

# Crop india_vec to match the expanded bounding box
india_vec <- st_crop(india_vec, india_bbox_sf)

# Remove temporary objects
rm(india_bbox, india_bbox_sf)

# Download the TerraClimate dataset on PDSI from 
# https://www.climatologylab.org/terraclimate.html
# india_rast_raw <- getTerraClim(
#   india_vec,
#   varname = "PDSI",
#   startDate = "2010-01-01",
#   endDate = "2023-12-01"
# )
# Use Downloaded Data
india_rast_raw <- rast("india_rast.tif")

# Match CRS of both: prefer to change CRS of vector
# Code for using directly downloaded data
# india_vec2 <- india_vec |> 
#   st_transform(crs(india_rast_raw$PDSI))
# Code for using downloaded data
india_vec2 <- india_vec |> 
  st_transform(crs(india_rast_raw))

# Crop and Maskt he Raster to show only the data within 
# India's Administrative Borders
# india_rast <- india_rast_raw$PDSI |> 
#   aggregate(fact = 4) |> 
#   terra::crop(india_vec2) |> 
#   terra::mask(india_vec2)

# Code for using downloaded data
india_rast <- india_rast_raw |> 
  aggregate(fact = 4) |> 
  terra::crop(india_vec2) |> 
  terra::mask(india_vec2)

# Clean Names for the Panel Strip Text: Month and Year
strip_labels <- india_rast |> 
  names() |> 
  str_remove("PDSI_") |> 
  str_remove("_total") |> 
  as_date() |> 
  format("%b\n%Y")

names(strip_labels) <- india_rast |> 
  names()

length(strip_labels)

# Temporarily save india_rast to save download time
# Check its range of values
# writeRaster(
#   india_rast,
#   filename = "india_rast.tif"
# )

# india_rast |> 
#   values() |> 
#   range(na.rm = T)

# A dummy tibble for writing labels along with maps
df_labels <- tibble(strip_label = strip_labels) |> 
  mutate(
    strip_label = str_replace_all(strip_label, "\n", " "),
    month = rep(1:12, 14),
    month_label = month(month, label = TRUE, abbr = TRUE),
    year_label  = as.character(rep(2010:2023, each = 12)),
    lyr = names(strip_labels)
  ) |> 
  mutate(
    month_label = if_else(year_label == "2010", month_label, ""),
    year_label = if_else(month == 1, year_label, "")
  )

# Some testing code to create facets with text only in few of them
# g1 <- ggplot(
#   data = df_labels
# ) +
#   geom_text(
#     data = df_labels |> filter(str_detect(lyr, "2010")),
#     mapping = aes(
#       label = month_label,
#       x = 82.5, y = 39
#     ),
#     size = 20,
#     family = "body_font"
#   ) +
#   geom_text(
#     data = df_labels |> filter(month == 1),
#     mapping = aes(
#       label = year_label,
#       x = 64, y = 22
#     ),
#     size = 20,
#     family = "body_font"
#   ) +
#   coord_sf(
#     default_crs = "EPSG:4326",
#     clip = "off",
#     xlim = c(68, 97),
#     ylim = c(7, 37)
#   ) +
#   facet_wrap(
#     ~lyr,
#     ncol = 12
#   ) +
#   theme_minimal() +
#   theme(
#     strip.text = element_blank(),
#     axis.text = element_blank(),
#     axis.title = element_blank(),
#     panel.grid = element_blank()
#   )
# 
# ggsave(
#   filename = here::here(
#     "geocomputation", "images",
#     "climateR_package_temp.png"
#   ),
#   plot = g1,
#   height = 2500 * 4.1,
#   width = 6000,
#   units = "px",
#   bg = "white",
#   limitsize = FALSE
# )
# 
# ggplot() +
#   geom_spatraster(
#     data = india_rast[[1]]
#   ) +
#   coord_sf(
#     default_crs = "EPSG:4326",
#     clip = "off",
#     xlim = c(68, 97),
#     ylim = c(7, 37)
#   )

g2 <- ggplot() +
  geom_spatraster(
    data = india_rast
  ) +
  geom_sf(
    data = india_vec2,
    fill = "transparent",
    linewidth = 0.25,
    colour = "grey20"
  ) +
  
  # Strip Labels along the months at top
  geom_text(
    data = df_labels |> filter(str_detect(lyr, "2010")),
    mapping = aes(
      label = month_label,
      x = 82.5, y = 45
    ),
    size = 50,
    family = "body_font"
  ) +
  
  # Strip labels for Year at the left side
  geom_text(
    data = df_labels |> filter(month == 1),
    mapping = aes(
      label = year_label,
      x = 52, y = 22
    ),
    size = 50,
    family = "body_font"
  ) +
  coord_sf(
    default_crs = "EPSG:4326",
    clip = "off",
    xlim = c(68, 97),
    ylim = c(7, 37)
  ) +
  paletteer::scale_fill_paletteer_c(
    "ggthemes::Orange-Blue-White Diverging",
    direction = 1,
    na.value = "transparent",
    limits = c(-10, 10),
    oob = scales::squish
  ) +
  facet_wrap(
    ~lyr,
    labeller = labeller(lyr = strip_labels),
    ncol = 12
  ) +
  labs(
    title = "India: Drought Severity (2010 - 2023)",
    subtitle = str_wrap("Palmer Drought Severity Index (PDSI) measures long-term drought and wetness conditions based on temperature, precipitation, and soil moisture balance. Negative values indicate drought severity, while positive values signify wetter-than-normal conditions.", 130),
    fill = "Palmer Drought Severity Index (PDSI)",
    caption = plot_caption
  ) +
  ggthemes::theme_map(
    base_size = 40,
    base_family = "body_font"
  ) +
  theme(
    legend.position = "bottom",
    legend.title.position = "top",
    # strip.text = element_text(
    #   hjust = 0.7,
    #   margin = margin(3,0,-20,0,"pt"),
    #   size = 48,
    #   lineheight = 0.3
    # ),
    strip.text = element_blank(),
    strip.background = element_blank(),
    plot.title.position = "plot",
    plot.title = element_text(
      margin = margin(0,0,10,0, "pt"),
      hjust = 0.5,
      size = 210
    ),
    plot.subtitle = element_text(
      margin = margin(0,0,50,0, "pt"),
      hjust = 0.5,
      size = 80,
      lineheight = 0.3
    ),
    plot.caption = element_textbox(
      hjust = 0.5,
      margin = margin(20,0,0,0, "pt"),
      size = 80
    ),
    panel.background = element_rect(
      fill = "transparent",
      colour = "transparent"
    ),
    legend.key.height = unit(10, "pt"),
    legend.key.width = unit(100, "pt"),
    legend.title = element_text(
      margin = margin(0,0,2,0, "mm"),
      size = 80,
      hjust = 0.5
    ),
    legend.text = element_text(
      margin = margin(2,0,0,0, "mm"),
      size = 60
    ),
    plot.margin = margin(10,5,10,90, "pt"),
    legend.margin =  margin(-20,0,0,0, "pt"),
    legend.box.margin = margin(-20,0,0,0, "pt"),
    legend.justification = c(0.5, 0)
  )

ggsave(
  filename = here::here(
    "geocomputation", "images",
    "climateR_package_2.png"
  ),
  plot = g2,
  height = 9000,
  width = 6200,
  units = "px",
  bg = "white",
  limitsize = FALSE
)

Figure 2: This faceted visualization displays the monthly PDSI values across India from 2010 to 2024. Each horizontal row represents an year, while each column represents a month, highlighting spatial variations in drought and wetness conditions. Negative values (shades of red) indicate drought severity, while positive values (shades of blue) reflect wetter-than-normal conditions. Data sourced from TerraClimate via the Climatology Lab.

References

Hijmans, Robert J. 2024. “Terra: Spatial Data Analysis.” https://CRAN.R-project.org/package=terra.

Pebesma, Edzer, and Roger Bivand. 2023. “Spatial Data Science: With Applications in r.” https://doi.org/10.1201/9780429459016.