Computing the Six ACI Components

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Each of the six ACI components can be computed independently of calculate_aci() — useful to inspect a single component, to recompute one after a data update, or simply to understand what each function does before combining them. This vignette walks through all of them on a small synthetic dataset built with the same pattern as vignette("xaci-intro").

Building a synthetic ERA5-like extract

We generate four years of hourly data (t2m, tp, u10, v10) on a tiny 3x3 grid, plus a land-mask NetCDF, entirely in memory / temp files — no network access required. We use the first three years (2011-2013) as the reference period and all four years (2011-2014) as the study period: this matters for what the standardised output looks like, see the note right after the data are generated below.

build_synthetic_nc <- function(path, var, unit, lon, lat, time_vec, origin, vals) {
  time_hours <- as.numeric(difftime(time_vec, origin, units = "hours"))
  dim_lon  <- ncdf4::ncdim_def("longitude", "degrees_east", lon)
  dim_lat  <- ncdf4::ncdim_def("latitude", "degrees_north", lat)
  dim_time <- ncdf4::ncdim_def(
    "time", paste0("hours since ", format(origin, "%Y-%m-%d %H:%M:%S")),
    time_hours, unlim = TRUE
  )
  ncvar <- ncdf4::ncvar_def(var, unit, list(dim_lon, dim_lat, dim_time),
                            missval = NA, prec = "double")
  nc <- ncdf4::nc_create(path, list(ncvar))
  ncdf4::ncvar_put(nc, ncvar, vals)
  ncdf4::nc_close(nc)
  invisible(path)
}

build_synthetic_mask <- function(path, lon, lat) {
  dim_lon <- ncdf4::ncdim_def("longitude", "degrees_east", lon)
  dim_lat <- ncdf4::ncdim_def("latitude", "degrees_north", lat)
  var_mask <- ncdf4::ncvar_def("country", "1", list(dim_lon, dim_lat),
                               missval = NA, prec = "double")
  nc <- ncdf4::nc_create(path, list(var_mask))
  ncdf4::ncvar_put(nc, var_mask, matrix(1, length(lon), length(lat)))
  ncdf4::nc_close(nc)
  invisible(path)
}

set.seed(42)
lon      <- c(-1, 0, 1)
lat      <- c(43, 44, 45)
origin   <- as.POSIXct("1900-01-01 00:00:00", tz = "UTC")
time_vec <- seq(as.POSIXct("2011-01-01 00:00", tz = "UTC"),
                as.POSIXct("2014-12-31 23:00", tz = "UTC"), by = "hour")
nlo <- length(lon); nla <- length(lat); nt <- length(time_vec)

# A mild warming/drying trend across the series, so that 2014 (the one
# out-of-reference year) shows genuine anomalies rather than pure noise.
trend <- seq_len(nt) / nt

# --- t2m: seasonal cycle in Kelvin, plus warming trend ---
seasonal_t <- 288 + 10 * sin(2 * pi * seq_len(nt) / (24 * 365)) + 0.6 * trend
t2m_vals <- array(NA_real_, c(nlo, nla, nt))
for (i in seq_len(nlo)) for (j in seq_len(nla))
  t2m_vals[i, j, ] <- seasonal_t + (i + j) + rnorm(nt, sd = 1.5)

# --- tp: precipitation (m/hour), mostly dry with occasional rain events,
# with rain becoming rarer over time (drying trend) ---
tp_vals <- array(0, c(nlo, nla, nt))
for (i in seq_len(nlo)) for (j in seq_len(nla)) {
  rain_hours <- rbinom(nt, 1, 0.08 * (1 - 0.3 * trend))
  tp_vals[i, j, ] <- rain_hours * rexp(nt, rate = 800)
}

# --- u10 / v10: wind components (m/s) ---
u10_vals <- array(rnorm(nlo * nla * nt, mean = 3, sd = 4), c(nlo, nla, nt))
v10_vals <- array(rnorm(nlo * nla * nt, mean = 1, sd = 4), c(nlo, nla, nt))

t2m_file  <- tempfile(fileext = ".nc")
tp_file   <- tempfile(fileext = ".nc")
u10_file  <- tempfile(fileext = ".nc")
v10_file  <- tempfile(fileext = ".nc")
mask_file <- tempfile(fileext = ".nc")

build_synthetic_nc(t2m_file, "t2m", "K",     lon, lat, time_vec, origin, t2m_vals)
build_synthetic_nc(tp_file,  "tp",  "m",     lon, lat, time_vec, origin, tp_vals)
build_synthetic_nc(u10_file, "u10", "m s-1", lon, lat, time_vec, origin, u10_vals)
build_synthetic_nc(v10_file, "v10", "m s-1", lon, lat, time_vec, origin, v10_vals)
build_synthetic_mask(mask_file, lon, lat)

reference_period <- c("2011-01-01", "2013-12-31")  # 3 years
study_period      <- c("2011-01-01", "2014-12-31")  # 4 years
country_abbrev    <- "XXX"   # fictitious code for this synthetic example

Three things worth knowing before reading the output below, all consequences of using a small, synthetic reference period rather than the multi-decade history you’d have with real data:

reference_period must span at least 2 years. Standardisation divides by a reference-period standard deviation computed separately for each calendar month; with only 1 year, every month has exactly 1 reference sample, an undefined (NA) standard deviation, and (for the sea-level component specifically) an empty result.
A handful of reference years makes every component look artificially alike. Standardising n reference observations against their own mean/sd always satisfies two constraints: they sum to 0, and their squares sum to n - 1. That leaves only n - 2 “degrees of freedom” for the actual data to express itself through. With n = 2 there are zero degrees of freedom left: every month’s anomaly is forced to exactly +1/sqrt(2) or -1/sqrt(2) (≈ ±0.71), for every component, regardless of what the underlying data actually looked like — which is why, with a 2-year reference, temperature, precipitation, drought and wind all “look the same” during the reference period: they’re not similar, they’re mathematically forced into the same shape. The reference period here spans 3 years (1 degree of freedom per month) rather than 2 or 5+: enough to break that exact-magnitude artifact and let each component’s standardised values differ from one another, while keeping this vignette quick to build. It’s a deliberate compromise for build speed, not a fully realistic reference sample — with real, multi-decade data this is a non-issue either way.
Months in study_period but outside reference_period (2014 here) aren’t constrained this way at all and show unambiguously genuine anomalies, driven by the warming/drying trend built into the synthetic data.

Temperature: `temperature_component()`

Both the hot-day frequency (\(T_{90}\)) and the cold-night frequency (\(T_{10}\)) come from the same function, called twice with different percentile / extremum / above_thresholds arguments:

t90 <- temperature_component(
  temperature_data_path = t2m_file,
  country_abbrev         = country_abbrev,
  reference_period        = reference_period,
  study_period            = study_period,
  mask_path               = mask_file,
  percentile              = 90, extremum = "max", above_thresholds = TRUE,
  area                    = TRUE
)

t10 <- temperature_component(
  temperature_data_path = t2m_file,
  country_abbrev         = country_abbrev,
  reference_period        = reference_period,
  study_period            = study_period,
  mask_path               = mask_file,
  percentile              = 10, extremum = "min", above_thresholds = FALSE,
  area                    = TRUE
)

head(t90, 3)   # inside reference_period (2011-2013): reflects the data
#> 2011-01-01 2011-02-01 2011-03-01 
#> -0.9149914 -0.9743106 -0.8320503
tail(t90, 3)   # 2014, outside reference_period: genuine anomalies
#> 2014-10-01 2014-11-01 2014-12-01 
#>   2.044340   1.236330   1.585021
head(t10, 3)
#> 2011-01-01 2011-02-01 2011-03-01 
#>  0.8992058  0.9424617  0.9635368
tail(t10, 3)
#> 2014-10-01 2014-11-01 2014-12-01 
#>  -1.999621  -1.826386  -1.225086

Internally, temperature_component() averages a “day” half and a “night” half of each day, computes the reference-period threshold at the requested percentile, then counts (monthly) how often each day/night crosses it.

Precipitation: `precipitation_component()`

The precipitation component is the maximum 5-day sliding sum of precipitation, standardised against the reference period:

prec <- precipitation_component(
  precipitation_data_path = tp_file,
  country_abbrev           = country_abbrev,
  reference_period          = reference_period,
  study_period              = study_period,
  mask_path                 = mask_file,
  var_name                  = "tp",
  window_size                = 5L,
  area                       = TRUE
)

head(prec, 3)
#> 2011-01-01 2011-02-01 2011-03-01 
#>  1.1428345  0.9358717  1.1215291
tail(prec, 3)
#> 2014-10-01 2014-11-01 2014-12-01 
#>  -2.496143  -1.857892  -1.149866

Drought: `drought_component()`

The drought component is based on the maximum number of consecutive dry days (CDD) per year, interpolated to monthly resolution and standardised:

drought <- drought_component(
  precipitation_data_path = tp_file,
  country_abbrev           = country_abbrev,
  reference_period          = reference_period,
  study_period              = study_period,
  mask_path                 = mask_file,
  area                       = TRUE
)

head(drought, 3)
#> 2011-01-01 2011-02-01 2011-03-01 
#> -0.9506522 -1.0762085 -1.1529024
tail(drought, 3)
#> 2014-10-01 2014-11-01 2014-12-01 
#>   1.529558   1.279779   1.080634

Wind: `wind_component()`

The wind component measures how often wind power (derived from the u10 and v10 components) exceeds its reference-period 90th percentile:

wind <- wind_component(
  wind_u10_data_path = u10_file,
  wind_v10_data_path = v10_file,
  country_abbrev       = country_abbrev,
  reference_period      = reference_period,
  study_period          = study_period,
  mask_path             = mask_file,
  area                  = TRUE
)

head(wind, 3)
#> 2011-01-01 2011-02-01 2011-03-01 
#>          0          0          0
tail(wind, 3)
#> 2014-10-01 2014-11-01 2014-12-01 
#>  0.1971326  0.2370370  0.2365591

Sea level: `sealevel_component()`

Unlike the other components, sea level comes from PSMSL tide-gauge records, not gridded NetCDF data. xaci bundles the PSMSL station metadata (names, coordinates, country) in inst/extdata/psmsl_data.csv, but not the tide-gauge measurement history itself (downloaded on demand by request_sealevel_data(), which requires network access).

To keep this vignette self-contained, we generate synthetic monthly measurements for two real French stations (Brest and Marseille, PSMSL IDs 1 and 61) in the format sealevel_load_data() expects, instead of downloading real records.

See the note above on why reference_period spans 2 full years here rather than 1: with a single reference year, sea-level standardisation would divide by an undefined (NA) per-month standard deviation, and — after that NA propagates through every standardised value — return an empty result.

psmsl_meta <- load_psmsl_data()
brest     <- psmsl_meta[psmsl_meta$ID == 1, ]
marseille <- psmsl_meta[psmsl_meta$ID == 61, ]
brest[, c("Station Name", "ID", "lat", "lon", "Country")]
#> # A tibble: 1 × 5
#>   `Station Name`    ID   lat   lon Country
#>   <chr>          <dbl> <dbl> <dbl> <chr>  
#> 1 BREST              1  48.4 -4.50 FRA
marseille[, c("Station Name", "ID", "lat", "lon", "Country")]
#> # A tibble: 1 × 5
#>   `Station Name`    ID   lat   lon Country
#>   <chr>          <dbl> <dbl> <dbl> <chr>  
#> 1 MARSEILLE         61  43.3  5.35 FRA

# PSMSL encodes the month as a year fraction (see ?sealevel_correct_date_format)
month_frac <- c("0417", "125", "2083", "2917", "375", "4583",
                "5417", "625", "7083", "7917", "875", "9583")

build_synthetic_psmsl_station <- function(dir, station_id, years,
                                          base_level, trend_mm_per_year) {
  lines <- character(0)
  for (y in years) {
    for (m in seq_len(12)) {
      level <- base_level + trend_mm_per_year * (y - years[1]) + rnorm(1, sd = 15)
      lines <- c(lines, sprintf("%d.%s;%.1f;0;000", y, month_frac[m], level))
    }
  }
  writeLines(lines, file.path(dir, paste0(station_id, ".txt")))
}

psmsl_dir <- tempfile("psmsl_")
dir.create(psmsl_dir)
set.seed(123)
build_synthetic_psmsl_station(psmsl_dir, 1,  2011:2014, base_level = 7020, trend_mm_per_year = 3)
build_synthetic_psmsl_station(psmsl_dir, 61, 2011:2014, base_level = 6980, trend_mm_per_year = 4)

sealevel_national <- sealevel_component(
  country_abbrev    = "FRA",     # must match the metadata's Country column
  study_period       = study_period,
  reference_period   = reference_period,
  area                = TRUE,
  sealevel_dir        = psmsl_dir
)

head(sealevel_national, 3)   # inside reference_period (2011-2013)
#>              sealevel
#> 2011-01-01 -1.1529915
#> 2011-02-01  0.3086872
#> 2011-03-01  0.8661490
tail(sealevel_national, 3)   # 2014, outside reference_period: genuine anomalies
#>               sealevel
#> 2014-10-01 -1.55288664
#> 2014-11-01  1.28824360
#> 2014-12-01  0.09876605

sealevel_component() also has a grid-cell mode (area = FALSE), which inverse-distance-weight-interpolates the station anomalies onto an ERA5 grid — this is what calculate_aci() uses internally when producing grid-cell (mapping) output:

sealevel_grid <- sealevel_component(
  country_abbrev  = "FRA",
  study_period     = study_period,
  reference_period = reference_period,
  area              = FALSE,
  mask_path         = mask_file,      # supplies the target lon/lat grid
  # Our fictitious grid (lon -1/0/1, lat 43/44/45) sits ~350-650 km from
  # Brest and Marseille; the default max_dist_km = 500 would leave a couple
  # of cells with no station in range (NA sea level -> NA ACI there). We
  # widen it here purely because this toy grid's placement is arbitrary --
  # with a real ERA5 grid this cutoff is meaningful and NA cells far from
  # the coast are expected, see ?sealevel_component.
  max_dist_km       = 800,
  sealevel_dir      = psmsl_dir
)

dim(sealevel_grid$data)
#> [1]  3  3 48

Caching heavy computations (steps 3-4)

With real ERA5 data, each *_component() call above can take a long time. Pass save = TRUE once to cache the result to disk, then computed_components = TRUE on subsequent calls to reload it instantly — this is exactly what calculate_aci() does internally.

save_dir/load_dir default to NULL (a sub-directory of tempdir(), cleared at the end of the session); pass your own persistent directory, as below, to actually benefit from the cache across sessions:

results_dir <- tools::R_user_dir("xaci", which = "data")

precipitation_component(
  precipitation_data_path = "data/era5/FRA/tp_2011_2015.nc",
  country_abbrev           = "FRA",
  reference_period          = c("2011-01-01", "2013-12-31"),
  study_period              = c("2011-01-01", "2015-12-31"),
  mask_path                 = "data/era5/FRA/mask_FRA.nc",
  area                       = FALSE,
  save                       = TRUE,
  save_dir                   = results_dir
)

# Later, or in a different session:
prec_national <- precipitation_component(
  precipitation_data_path = "data/era5/FRA/tp_2011_2015.nc",
  country_abbrev           = "FRA",
  reference_period          = c("2011-01-01", "2013-12-31"),
  study_period              = c("2011-01-01", "2015-12-31"),
  mask_path                 = "data/era5/FRA/mask_FRA.nc",
  area                       = TRUE,
  computed_components        = TRUE,
  load_dir                   = results_dir
)

The next vignette, vignette("xaci-full-pipeline"), shows how calculate_aci() combines all six components above into the index itself.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.

Computing the Six ACI Components

Building a synthetic ERA5-like extract

Temperature: temperature_component()

Precipitation: precipitation_component()

Drought: drought_component()

Wind: wind_component()

Sea level: sealevel_component()

Caching heavy computations (steps 3-4)

Temperature: `temperature_component()`

Precipitation: `precipitation_component()`

Drought: `drought_component()`

Wind: `wind_component()`

Sea level: `sealevel_component()`