Administrative-Level Aggregation

library(xaci)

Besides the national scalar mode (area = TRUE) and the grid-cell mode (area = FALSE, see vignette("xaci-visualization")), xaci can aggregate components and the ACI itself at the level of administrative units — e.g. French departments or regions — by passing admin_level to calculate_aci() or to any individual *_component() function.

How it works

Administrative aggregation relies on two pieces of geographic data, both fetched over the network the first time they are needed:

Both of these require network access and are therefore not run in this vignette. In an interactive session, the call looks like this:

admin_mask_L1 <- build_admin_mask(
  lon = grid_aci$lon, lat = grid_aci$lat,
  country_abbrev = "FRA", admin_level = 1, crs_metric = 2154   # Lambert-93
)

dept_assignment <- assign_sealevel_to_admin(
  country_abbrev = "FRA", admin_level = 1, crs_metric = 2154
)

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

dept_aci_FRA <- calculate_aci(
  country_abbrev      = "FRA",
  study_period          = c("2011-01-01", "2015-12-31"),
  reference_period       = c("2011-01-01", "2013-12-31"),
  years                    = 2011:2015,
  granularity               = "month",
  admin_level                = 1,
  crs_metric                  = 2154,
  load_dir                     = results_dir,
  computed_components           = TRUE
)

dim(dept_aci_FRA)   # months x (7 variables x number of departments)
plot_aci_map(dept_aci_FRA, variable = "ACI", time_index = "mean")

Passing admin_mask (and, for sea level, admin_assignment) explicitly to an individual component avoids rebuilding it on every call — useful when computing several components for the same country/level:

prec_admin_L1 <- 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,
  admin_mask                    = admin_mask_L1,
  computed_components            = TRUE,
  load_dir                        = results_dir
)

The aggregation mechanics, without the network

The two GADM/rnaturalearth downloads above only produce two simple data structures — an admin_mask (area-fraction weights per grid cell per administrative unit) and an admin_assignment (station-to-unit mapping and coastal fractions). Everything downstream of them is ordinary R and can be demonstrated with hand-built versions of these structures, exactly as the package’s own test suite does to keep tests network-free.

A synthetic two-region country

lon <- c(-1, 0)
lat <- c(43, 44)

# A minimal admin_mask: 2 grid cells, cell (i, j=1) entirely in region "A",
# cell (i, j=2) entirely in region "B" (weight 1 -- no partial overlap here,
# to keep the arithmetic easy to follow).
synthetic_admin_mask <- list(
  units   = c("A", "B"),
  lon     = lon,
  lat     = lat,
  weights = list(
    `1` = c(A = 1, B = 0),   # cell (lon[1], lat[1])
    `2` = c(A = 0, B = 1),   # cell (lon[1], lat[2])
    `3` = c(A = 1, B = 0),   # cell (lon[2], lat[1])
    `4` = c(A = 0, B = 1)    # cell (lon[2], lat[2])
  )
)

We reuse a small synthetic temperature component (built as in vignette("xaci-components")) and aggregate it per region with reduce_dataarray_to_dataframe() — the same function calculate_aci() calls internally once admin_mask is known:

build_synthetic_t2m <- function(path, lon, lat, time_vec, origin) {
  time_hours <- as.numeric(difftime(time_vec, origin, units = "hours"))
  nlo <- length(lon); nla <- length(lat); nt <- length(time_vec)
  set.seed(5)
  trend    <- seq_len(nt) / nt   # mild warming trend, see note below
  seasonal <- 288 + 10 * sin(2 * pi * seq_len(nt) / (24 * 365)) + 0.6 * trend
  vals <- array(NA_real_, dim = c(nlo, nla, nt))
  for (i in seq_len(nlo)) {
    for (j in seq_len(nla)) {
      # A clear north/south offset so the two regions differ visibly
      vals[i, j, ] <- seasonal + 4 * j + rnorm(nt, sd = 1.5)
    }
  }
  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
  )
  var_t2m <- ncdf4::ncvar_def("t2m", "K", list(dim_lon, dim_lat, dim_time),
                              missval = NA, prec = "double")
  nc <- ncdf4::nc_create(path, list(var_t2m))
  ncdf4::ncvar_put(nc, var_t2m, vals)
  ncdf4::nc_close(nc)
  invisible(path)
}

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")
t2m_file <- tempfile(fileext = ".nc")
build_synthetic_t2m(t2m_file, lon, lat, time_vec, origin)

# reference_period spans 3 years (not fewer) for the same reason as the
# other vignettes: with only 1, every month has a single reference sample
# and standardises to exactly 0; with exactly 2, standardising always
# collapses to +-1/sqrt(2) (a mathematical identity, see
# vignette("xaci-components")) -- informative for neither region here. 3
# reference years (1 degree of freedom per month) is a build-speed
# compromise that still lets regions A and B genuinely differ, though a
# longer reference would be more realistic. study_period extends one year
# beyond it (2014), which isn't constrained this way and shows
# unambiguously genuine anomalies, driven by the warming trend in the
# synthetic data.
reference_period <- c("2011-01-01", "2013-12-31")
study_period      <- c("2011-01-01", "2014-12-31")

t90_grid <- temperature_component(
  temperature_data_path = t2m_file,
  country_abbrev          = "XXX",
  reference_period         = reference_period,
  study_period             = study_period,
  percentile               = 90, extremum = "max", above_thresholds = TRUE,
  area                     = FALSE,     # keep the spatial dimension
  admin_mask                = synthetic_admin_mask
)

head(t90_grid, 3)   # inside reference_period (2011-2013)
#>                 t90_A      t90_B
#> 2011-01-01 -0.7823697 -0.9558920
#> 2011-02-01 -0.6112760 -0.7125328
#> 2011-03-01 -0.7296229 -0.5299002
tail(t90_grid, 3)   # 2014, outside reference_period: genuine anomalies
#>               t90_A     t90_B
#> 2014-10-01 2.672668 0.7958107
#> 2014-11-01 4.458087 0.7962436
#> 2014-12-01 2.198950 1.5060550

temperature_component() returned one column per region directly, because we passed admin_mask in. Internally, each grid cell is first standardised against its own reference-period distribution, and only then averaged (weighted by area fraction) into region A and region B — so the two columns can differ, but only insofar as the two regions’ underlying cells behave differently, not because of a raw, unstandardised temperature gap between them.

The sea-level equivalent

assign_sealevel_to_admin()’s output — list(station_ids, factors) — can be built by hand the same way, using real PSMSL station IDs (see vignette("xaci-components") for how to generate synthetic tide-gauge .txt files for these):

synthetic_admin_assignment <- list(
  station_ids = list(A = 1, B = 61),      # PSMSL IDs: Brest -> A, Marseille -> B
  factors     = c(A = 0.6, B = 0.3)       # coastal fraction per region
)

Passed as admin_assignment to sealevel_component(..., area = TRUE), this drives reduce_sealevel_over_region() to produce one sealevel_A / sealevel_B column instead of a single national series — see ?reduce_sealevel_over_region for the exact aggregation logic, which mirrors reduce_dataarray_to_dataframe() above but on tide-gauge station data instead of a spatial grid.