Compute anomalies from monthly data and historic normals

[mikejohnson51]

Question from email:

---

What I'm trying to do is calculate anomalies between a "current" period  (mean values for 2004:2015) and a "baseline" period (1961-1990 normals), for Tmin coldest month, Tmax warmest month, climate water deficit (perhaps summed over the year), etc.

So for example, to calculate the anomaly for Tmin coldest month, for each site I need to calculate:

1. mean(Tmin.Jan), mean(Tmin.Feb), etc., across 2004-2015
2. the min. value from 1 (since we don't know which month is the coldest)
3. The min value of (Tmin.Jan., Tmin.Feb, etc.) from 1961-1990 normals
4. difference between 2 and 3

---

I think there are two approaches for this - a raster and a point based solution. The raster approach is simpler I find but will provide both for comparison:

library(raster)
library(AOI)
library(sf)
library(climateR)
library(dplyr)

Get Data for 100 random pts in Colorado ...

random_pts = AOI::aoi_get(state = "co") %>% 
  st_sample(100) %>% 
  st_cast("POINT") %>% 
  st_as_sf() %>% 
  mutate(ID = 1:n())

co = getTerraClim(random_pts, 
                  param = "tmin", 
                  startDate = "2004-01-01", 
                  endDate = "2014-12-01")

coNorm = getTerraClimNormals(random_pts, 
                             param = "tmin", 
                             period = "19611990")

Site based approach

ext = extract_sites(r = co, pts = random_pts, "ID")

ext_norm = extract_sites(r = coNorm, pts = random_pts, "ID")[[1]] %>% 
  tidyr::pivot_longer(-date, values_to = "norms") %>% 
  mutate(month = as.numeric(date), date = NULL) %>% 
  group_by(name) %>% 
  summarize(minNorm = min(norms)) %>% 
  ungroup()

out = ext$tmin %>% 
  tidyr::pivot_longer(-date) %>% 
  mutate(month = lubridate::month(date)) %>% 
  group_by(month, name) %>% 
  summarise(meanTmin = mean(value)) %>% 
  ungroup() %>% 
  group_by(name) %>% 
  summarise(minTmin = min(meanTmin)) %>%
  ungroup() %>% 
  left_join(ext_norm, by = c("name")) %>% 
  mutate(site_anom = minTmin - minNorm, ID = as.numeric(gsub("site_", "", name))) %>% 
  left_join(random_pts, by = "ID") %>% 
  st_as_sf() %>% 
  arrange(ID)
#> `summarise()` has grouped output by 'month'. You can override using the `.groups` argument.

plot(out['site_anom'], pch = 16)

Raster Based approach

indices      = rep(1:12, times = nlayers(co$tmin)/12)
co_tmin      = stackApply(co$tmin, indices = indices, mean)
diff_rast = min(co_tmin) - min(coNorm$`19611990_tmin`)
random_pts$raster_anom =  extract(diff_rast, random_pts)
plot(diff_rast)

Check for agreement

plot(random_pts$raster_anom, out$site_anom)
abline(0,1)

^{Created on 2021-03-10 by the reprex package (v0.3.0)}