TerraclimateR

The terraclimate data set¹ provides data on 14 climatic variables through time since 1958 at ~4km resolution in space and monthly resolution in time. TerraclimateR is a lightweight package that implements a single function (get_terraclim) to access these data and import them directly into R in long format.

General information:

https://www.climatologylab.org/terraclimate.html

Usage

Code

library(TerraclimateR)
suppressPackageStartupMessages(library(tidyverse))

Input data (df) is required in the following format: a data.frame (or matrix) containing at least three variables (columns): a (unique) identifier, latitude, and longitude:

Code

df <- tidyr::tribble(
  ~id,            ~lat,       ~lon,
  "Vilcabamba",    -4.260641, -79.22274,
  "Paracas",      -13.833881, -76.25046,
  "Puerto Viejo",   9.656562, -82.75600,
  "Trinidad",      21.808230, -79.98102,
  "Cartagena",     10.172151, -75.74658
)

To retrieve climatic data for these observations (rows), a character containing selected terraclimate variables is passed to the clim_vars argument (e.g. c("aet", "def")). The follwing variables are available:

Available terraclimate variables:

Code	Variable	Detail	Unit
aet	Actual Evapotranspiration	monthly total	\(mm\)
def	Climate Water Deficit	monthly total	\(mm\)
pet	Potential evapotranspiration	monthly total	\(mm\)
ppt	Precipitation	monthly total	\(mm\)
q	Runoff	monthly total	\(mm\)
soil	Soil Moisture	total column at end of month	\(mm\)
srad	Downward surface shortwave radiation		\(W/m^2\)
swe	Snow water equivalent	at end of month	\(mm\)
tmax	Max Temperature	average for month	\(C\)
tmin	Min Temperature	average for month	\(C\)
vap	Vapor pressure	average for month	\(kPa\)
ws	Wind speed	average for month	\(m/s\)
vpd	Vapor Pressure Deficit	average for month	\(kPa\)
PDSI	Palmer Drought Severity Index	at end of month	unitless

To retrieve all of these variables, pass climvars = "ALL".

Example:

Retrieving data for minimum temperature per month (tmin) and monthly precipitation (ppt):

Code

nc_df <- get_terraclim(df        = df,
                       id_var    = id,
                       lon_var   = lon,
                       lat_var   = lat,
                       clim_vars = c("tmin", "ppt"),
                       show_prog = TRUE,
                       conserve  = TRUE)

The column id is retained and still contains our sample identifiers. The column year contains values ranging from 1958, the start of the terraclimate data set, to the latest year for which data are available. month specifies the month, and variable indicates the terraclimate variable. value is the value of the respective variable at the given month in the given year at the given location. Because conserve = TRUE was set in the function call, all the columns of the initial data.frame are still present. If the data.frame contains a lot of variables and hundreds of observations, it is advisable to set conserve = FALSE since this results in a lot of redundant information and may slow down computation. Also, if column names year, month, variable, or value are in original data they will be renamed (e.g. year to year.x and year.y) to guarantee unique variable names. Returning the data in long format facilitates further processing in a tidy style (see tidy data).

Code

str(nc_df)

## 'data.frame':    7800 obs. of  7 variables:
##  $ id      : chr  "Cartagena" "Cartagena" "Cartagena" "Cartagena" ...
##  $ year    : int  2022 2022 2022 2022 2022 2022 2022 2022 2022 2022 ...
##  $ month   : Factor w/ 12 levels "Jan","Feb","Mar",..: 4 3 8 5 6 7 12 9 10 11 ...
##  $ variable: Factor w/ 2 levels "ppt","tmin": 1 1 1 1 1 1 1 1 1 1 ...
##  $ value   : num  200.5 4.6 217 213.7 281.6 ...
##  $ lat     : num  10.2 10.2 10.2 10.2 10.2 ...
##  $ lon     : num  -75.7 -75.7 -75.7 -75.7 -75.7 ...

Retrieving data for only one year:

Code

get_terraclim(df        = df[1, ],
              id_var    = id,
              lon_var   = lon,
              lat_var   = lat,
              clim_vars = "ppt",
              show_prog = FALSE,
              conserve  = FALSE,
              year      = 2021)

##            id year month variable value
## 1  Vilcabamba 2021   Jan      ppt 102.9
## 2  Vilcabamba 2021   Feb      ppt 125.8
## 3  Vilcabamba 2021   Mar      ppt 336.7
## 4  Vilcabamba 2021   Apr      ppt 162.9
## 5  Vilcabamba 2021   May      ppt  54.8
## 6  Vilcabamba 2021   Jun      ppt  41.4
## 7  Vilcabamba 2021   Jul      ppt  20.3
## 8  Vilcabamba 2021   Aug      ppt  14.9
## 9  Vilcabamba 2021   Sep      ppt  20.8
## 10 Vilcabamba 2021   Oct      ppt 101.2
## 11 Vilcabamba 2021   Nov      ppt  72.3
## 12 Vilcabamba 2021   Dec      ppt  70.4

Retrieving data for a range of years:

Code

nc_df_1990_2000 <- get_terraclim(df        = df[1, ],
                                 id_var    = id,
                                 lon_var   = lon,
                                 lat_var   = lat,
                                 clim_vars = "ppt",
                                 show_prog = FALSE,
                                 conserve  = FALSE) %>% 
  dplyr::filter(year %in% 1990:2000)

str(nc_df_1990_2000)

## 'data.frame':    132 obs. of  5 variables:
##  $ id      : chr  "Vilcabamba" "Vilcabamba" "Vilcabamba" "Vilcabamba" ...
##  $ year    : int  1990 1990 1990 1990 1990 1990 1990 1990 1990 1990 ...
##  $ month   : Factor w/ 12 levels "Jan","Feb","Mar",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ variable: Factor w/ 1 level "ppt": 1 1 1 1 1 1 1 1 1 1 ...
##  $ value   : num  82.3 129.9 115.9 107.6 44.9 ...

NA Handling

Missing values, i.e., NA’s, are not discarded by get_terraclim.

Code

incomplete_df <- tribble(
  ~id,            ~lat,       ~lon,
  "Vilcabamba",    -4.260641, -79.22274,
  NA,             -13.833881, -76.25046,
  "Puerto Viejo",         NA, -82.75600,
  "Trinidad",             NA,        NA,
  NA,                     NA,        NA
)

Code

incomplete_nc_df <- get_terraclim(df        = incomplete_df,
                                  id_var    = id,
                                  lon_var   = lon,
                                  lat_var   = lat,
                                  clim_vars = c("tmin"),
                                  conserve  = TRUE)

incomplete_nc_df %>% group_by(id, lat, lon) %>% reframe(data_length = n())

## # A tibble: 5 × 4
##   id              lat   lon data_length
##   <chr>         <dbl> <dbl>       <int>
## 1 NA.2         -13.8  -76.3         780
## 2 NA.5          NA     NA             1
## 3 Puerto Viejo  NA    -82.8           1
## 4 Trinidad      NA     NA             1
## 5 Vilcabamba    -4.26 -79.2         780

Code

summary(incomplete_nc_df)

##       id                 year          month     variable        value      
##  Length:1563        Min.   :1958   Jan    :130   tmin:1560   Min.   :10.84  
##  Class :character   1st Qu.:1974   Feb    :130   NA's:   3   1st Qu.:12.77  
##  Mode  :character   Median :1990   Mar    :130               Median :13.51  
##                     Mean   :1990   Apr    :130               Mean   :14.18  
##                     3rd Qu.:2006   May    :130               3rd Qu.:14.88  
##                     Max.   :2022   (Other):910               Max.   :21.47  
##                     NA's   :3      NA's   :  3               NA's   :3      
##       lat               lon        
##  Min.   :-13.834   Min.   :-82.76  
##  1st Qu.:-13.834   1st Qu.:-79.22  
##  Median : -9.047   Median :-79.22  
##  Mean   : -9.047   Mean   :-77.74  
##  3rd Qu.: -4.261   3rd Qu.:-76.25  
##  Max.   : -4.261   Max.   :-76.25  
##  NA's   :3         NA's   :2

But exclusion of missing values is easy:

Code

incomplete_nc_df %>% drop_na() %>% summary()

##       id                 year          month     variable        value      
##  Length:1560        Min.   :1958   Jan    :130   tmin:1560   Min.   :10.84  
##  Class :character   1st Qu.:1974   Feb    :130               1st Qu.:12.77  
##  Mode  :character   Median :1990   Mar    :130               Median :13.51  
##                     Mean   :1990   Apr    :130               Mean   :14.18  
##                     3rd Qu.:2006   May    :130               3rd Qu.:14.88  
##                     Max.   :2022   Jun    :130               Max.   :21.47  
##                                    (Other):780                              
##       lat               lon        
##  Min.   :-13.834   Min.   :-79.22  
##  1st Qu.:-13.834   1st Qu.:-79.22  
##  Median : -9.047   Median :-77.74  
##  Mean   : -9.047   Mean   :-77.74  
##  3rd Qu.: -4.261   3rd Qu.:-76.25  
##  Max.   : -4.261   Max.   :-76.25  
##

Integration into Tidy Workflow

Transformations (Examples)

Code

nc_df %>% 
  group_by(id, month, variable) %>% 
  summarise(monthly_mean = mean(value), .groups = "drop") %>% 
  pivot_wider(values_from = monthly_mean, names_from = variable)

## # A tibble: 60 × 4
##    id        month     ppt  tmin
##    <chr>     <fct>   <dbl> <dbl>
##  1 Cartagena Jan     3.22   22.7
##  2 Cartagena Feb     0.778  22.9
##  3 Cartagena Mar     5.83   23.5
##  4 Cartagena Apr   103.     24.1
##  5 Cartagena May   160.     24.3
##  6 Cartagena Jun   160.     24.4
##  7 Cartagena Jul   117.     24.2
##  8 Cartagena Aug   184.     24.4
##  9 Cartagena Sep   149.     24.3
## 10 Cartagena Oct   240.     23.9
## # ℹ 50 more rows

Code

nc_df %>% group_by(id, year, variable) %>% summarize(yearly_sum = sum(value), .groups = "drop")

## # A tibble: 650 × 4
##    id         year variable yearly_sum
##    <chr>     <int> <fct>         <dbl>
##  1 Cartagena  1958 ppt           1105.
##  2 Cartagena  1958 tmin           280.
##  3 Cartagena  1959 ppt            960.
##  4 Cartagena  1959 tmin           270.
##  5 Cartagena  1960 ppt           1525.
##  6 Cartagena  1960 tmin           276.
##  7 Cartagena  1961 ppt           1256.
##  8 Cartagena  1961 tmin           273.
##  9 Cartagena  1962 ppt           1262 
## 10 Cartagena  1962 tmin           270.
## # ℹ 640 more rows

Visualisation

It’s always convenient to plot the data in order to quickly identify incorrect coordinates, e.g., inverted values.

Code

world_map <- ggplot2::map_data("world")

dplyr::distinct(world_map, region) %>%
  ggplot() +
  geom_map(map = world_map, aes(map_id = region), 
           fill = "gray90", color = "gray20", linewidth = 0.2) +
  geom_point(data = df, mapping = aes(lon, lat), fill = "#3366ff", size = 3, shape = 21) +
  geom_label(data = df %>% mutate(lon = case_when(id == "Puerto Viejo" ~ lon - 3, 
                                                  id == "Cartagena" ~ lon + 3,
                                                  .default = lon)), 
             mapping = aes(lon, lat - 2.5, label = id)) +
  expand_limits(x = c(min(df$lon) - 15, max(df$lon) + 15), 
                y = c(min(df$lat) - 5, max(df$lat) + 5)) +
  coord_equal(expand = FALSE) +
  labs(x = "Longitude", y = "Latitude", title = "Geographic Locations of Observations") +
  theme_linedraw()

Code

nc_df %>% filter(id == "Cartagena" & variable == "tmin") %>% 
ggplot() +
  aes(group = year, x = month, y = value, color = year) +
  geom_line(linewidth = 1.5) +
  scale_color_viridis_c(option = "rocket", name = "Year") +
  labs(x = "Month", y = expression(Temp[min]~"[°C]"), title = "Yearly Course of Minimum Temperature in Cartagena, CO") +
  theme_dark() +
  theme(legend.key.height = unit(0.8, "in"))

Code

nc_df %>% 
  filter(id != "Paracas" & variable == "ppt") %>% 
  group_by(id, variable, month) %>% 
  summarise(monthly_mean = mean(value), se = sd(value)/sqrt(n()), .groups = "drop") %>% 
  ggplot() +
  geom_col(aes(x = month, y = monthly_mean, fill = id), color = "gray20") +
  geom_errorbar(aes(x = month, ymin = monthly_mean - se, ymax = monthly_mean + se), 
                color = "gray20", width = 0.3) +
  facet_wrap(. ~ id) +
  labs(x = "Month", y = "Precipitation [mm]", title = "Average Precipitation per Month") +
  scale_fill_brewer(palette = "Set1") +
  theme_classic() +
  theme(legend.position = "none", strip.background = element_blank())

Code

nc_df %>% 
  filter(variable == "ppt") %>% 
  group_by(id, variable, year) %>% 
  summarise(yearly_sum = sum(value), .groups = "drop") %>% 
  filter(id == "Vilcabamba") %>% 
  ggplot() +
  geom_smooth(aes(x = year, y = yearly_sum), 
              color = "#3366ff", linewidth = 1.5, method = lm, formula = "y ~ x") +
  geom_point(aes(x = year, y = yearly_sum), color = "gray20") +
  geom_line(aes(x = year, y = yearly_sum), color = "gray20") +
  labs(x = "Year", y = "Precipitation per Year [mm]", title = "Yearly Precipitation in Vilcabamba") +
  theme_classic()

Additional Information

There are currently two helper functions used internally in TerraclimateR. If you are interested you get further information:

Code

?TerraclimateR:::get_terraclim_single()

?TerraclimateR:::show_progress()

Abatzoglou, J., Dobrowski, S., Parks, S. et al. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. Sci Data 5, 170191 (2018). https://doi.org/10.1038/sdata.2017.191 ↩︎