Package 'xaci'

Description

Once data has been reduced to daily (or coarser) resolution, it is small enough to hand off to the rest of the existing base-R pipeline (resample_monthly(), standardize_metric(), .compute_aci_grid(), etc.) unchanged. This is the bridge between the terra-based loading/reduction steps and that pipeline.

Usage

.spatraster_to_list(r)

Arguments

Value

A list with elements data ([lon x lat x time] array), lon, lat, time – the same structure produced by load_netcdf().

Description

Computes the full Actuarial Climate Index by combining all five components: temperature (T90 and T10), precipitation, drought, wind, and sea level.

Description

Aggregate a monthly data frame to a given temporal granularity

Usage

aggregate_granularity(df, granularity = "month")

Arguments

Value

A data.frame with aggregated values and appropriate row names.

Description

Produces a GIF animation (via gganimate) showing the evolution of a spatial variable at each time step. Requires gganimate and gifski (or png) to be installed.

Usage

animate_aci_map(
  data,
  variable = "ACI",
  country_abbrev = NULL,
  admin_level = NULL,
  crs_metric = 4326,
  var_label = "Standardised anomaly",
  title = "ACI over time",
  palette = "RdBu",
  reverse = TRUE,
  n_breaks = 9,
  borders = TRUE,
  fps = 4,
  width = 800,
  height = 600,
  save_path = NULL
)

Arguments

Value

A gganimate animation object (invisibly).

Why dontrun

This example is wrapped in dontrun because it depends on the optional gganimate and gifski packages (declared in Suggests, not guaranteed to be installed on every machine running R CMD check), and because rendering a multi-frame GIF takes well over 5 seconds, making it unsuitable for donttest as well.

Examples

## Not run: 
lon <- seq(-5, 8, length.out = 10)
lat <- seq(42, 51, length.out = 8)
arr <- array(stats::rnorm(10 * 8 * 5), dim = c(10, 8, 5))
grid <- list(ACI = arr, lon = lon, lat = lat,
            time = seq_len(5))

animate_aci_map(grid, variable = "ACI", country_abbrev = "FRA",
                fps = 2, save_path = file.path(tempdir(), "aci_animation.gif"))

## End(Not run)

Description

Sets grid cells to NA where the mask value is below threshold.

Usage

apply_mask(dataset, mask_path, var_name, threshold = 0.8)

Arguments

Value

The dataset list with data masked in-place.

Description

Equivalent of apply_mask(): sets cells to NA where the mask value is below threshold. Processed block-by-block by terra.

Usage

apply_mask_terra(r, mask_path, threshold = 0.8, chunk_size = 20000)

Arguments

Details

Note sur la limite de 65535 couches : le format interne utilise par terra::mask() pour ecrire son resultat sur disque (fichier temporaire) ne supporte pas plus de 65535 couches en sortie. Au-dela (cas frequent en donnees horaires sur plusieurs decennies : ~300k couches pour 35 ans), on ne peut pas masquer l'objet en un seul appel. Cette fonction bascule donc automatiquement sur un traitement par blocs temporels ('chunk_size' couches a la fois), ecrits en GeoTIFF (BigTIFF) puis recombines en une seule source multi-fichiers, sans jamais materialiser l'ensemble en RAM. Le masque etant purement spatial (identique a chaque pas de temps), le decoupage en blocs ne change pas le resultat : un pixel exclu par le masque l'est de la meme facon dans chaque bloc.

Precision numerique : terra ecrit sur disque en simple precision (datatype = "FLT4S") PAR DEFAUT (voir terra::terraOptions()), des que le resultat ne tient pas en RAM – ce qui inclut la branche par blocs ci-dessous (TOUJOURS ecrite sur disque). Sans le datatype = "FLT8S" explicite ci-dessous, un masquage qui declenche une ecriture disque perdrait silencieusement de la precision par rapport a un masquage qui reste en RAM – incoherence potentielle non liee au chunking lui-meme, mais au meme risque general (voir la note similaire dans resample_daily_terra()).

Value

The masked terra::SpatRaster.

Description

Performs a spatial join between PSMSL station coordinates and administrative unit polygons, and computes for each unit the fraction of its perimeter that is coastline. The coastline layer is cropped to the country bounding box before intersection to speed up computation. All geometries are projected to crs_metric before length calculations to ensure results are in metres.

Usage

assign_sealevel_to_admin(country_abbrev, admin_level = 1, crs_metric = 4326)

Arguments

Value

A list with two elements:

station_ids

Named list: keys are administrative unit names, values are integer vectors of PSMSL station IDs within that unit.

factors

Named numeric vector: keys are administrative unit names, values are the coastal fraction (coastline length / total perimeter) in [0, 1]. Zero for landlocked units.

Description

For each ERA5 grid cell, computes the fraction of its surface area that falls within each administrative unit polygon (using sf geometry intersection). All geometries are projected to crs_metric before area calculations to ensure results are in metres squared. The result can be computationally expensive for fine grids and large countries — consider caching it with saveRDS().

Usage

build_admin_mask(
  lon,
  lat,
  country_abbrev,
  admin_level = 1,
  resolution = 0.25,
  crs_metric = 4326,
  cache_dir = .gadm_cache_dir()
)

Arguments

Value

A list with elements:

weights

A named list of length length(lon) * length(lat). Each element is a named numeric vector of (unit -> fractional area weight) pairs.

units

Character vector of all administrative unit names.

lon

Input longitudes.

lat

Input latitudes.

Description

This is the main entry point of the package. It instantiates all five climate components, standardises them over the reference period, and combines them into a single ACI time series:

Usage

calculate_aci(
  country_abbrev,
  study_period,
  reference_period,
  years = NULL,
  temperature_data_path = NULL,
  precipitation_data_path = NULL,
  wind_u10_data_path = NULL,
  wind_v10_data_path = NULL,
  mask_data_path = NULL,
  sealevel_dir = NULL,
  percentile_high = 90,
  percentile_low = 10,
  granularity = "month",
  area = TRUE,
  factor = 1/5,
  max_dist_km = 500,
  admin_level = NULL,
  crs_metric = 4326,
  save = FALSE,
  save_dir = NULL,
  load_dir = NULL,
  computed_components = FALSE,
  engine = c("base", "terra"),
  cores = 1L
)

Arguments

Details

$ACI = \frac{T_{high} - T_{low} + P + D + \alpha \cdot SL + W}{5 + \alpha}$

where $\alpha$ is the sea-level erosion factor (coastal fraction, default 1/5). For administrative units without coastal stations, $\alpha = 0$ and the denominator becomes 5. At grid-cell level, $\alpha \in \{0, 1\}$ depending on whether the cell lies within max_dist_km of a tide-gauge station.

Value

National scalar (area = TRUE, admin_level = NULL)

A data.frame with columns t<percentile_high>, t<percentile_low>, precipitation, drought, wind, sealevel, and ACI, indexed by dates at the chosen granularity.

Grid-cell (area = FALSE, admin_level = NULL)

A named list with one array [lon x lat x nt] per component (ACI, t<percentile_high>, t<percentile_low>, precipitation, drought, wind, sealevel), plus lon, lat, and time (character vector of period labels at the chosen granularity).

Administrative (admin_level integer)

A data.frame indexed by dates at the chosen granularity, with one <component>_<unit> column per administrative unit FOR EACH component (ACI, t<percentile_high>, t<percentile_low>, precipitation, drought, wind, sealevel) – e.g. ACI_<unit>, t90_<unit>, precipitation_<unit>, etc. This mirrors the grid-cell mode's per-component structure, and lets plot_aci_map() plot any individual component at admin level, not just ACI itself.

Examples

The four examples below all require real ERA5 NetCDF files (temperature, precipitation, wind, country mask), obtained beforehand via download_era5/download_mask and a Copernicus CDS API key (see cds_set_key). They cannot be run automatically by R CMD check and are therefore wrapped in dontrun: (1) national ACI at annual granularity; (2) grid-cell ACI (full spatial output); (3) ACI with custom percentiles (T95/T5); and (4) ACI by department (admin level 2), using Lambert-93 for France.

Examples

## Not run: 
result <- calculate_aci(
  temperature_data_path   = "t2m_1960-2020.nc",
  precipitation_data_path = "tp_1960-2020.nc",
  wind_u10_data_path      = "u10_1960-2020.nc",
  wind_v10_data_path      = "v10_1960-2020.nc",
  country_abbrev          = "FRA",
  mask_data_path          = "mask_FRA.nc",
  study_period            = c("1980-01-01", "2020-12-31"),
  reference_period        = c("1961-01-01", "1990-12-31"),
  granularity             = "year",
  area                    = TRUE,
  factor                  = 1/5
)

grid <- calculate_aci(
  temperature_data_path   = "t2m_1960-2020.nc",
  precipitation_data_path = "tp_1960-2020.nc",
  wind_u10_data_path      = "u10_1960-2020.nc",
  wind_v10_data_path      = "v10_1960-2020.nc",
  country_abbrev          = "FRA",
  mask_data_path          = "mask_FRA.nc",
  study_period            = c("1980-01-01", "2020-12-31"),
  reference_period        = c("1961-01-01", "1990-12-31"),
  granularity             = "year",
  area                    = FALSE,
  max_dist_km             = 500
)

result_custom <- calculate_aci(
  temperature_data_path   = "t2m_1960-2020.nc",
  precipitation_data_path = "tp_1960-2020.nc",
  wind_u10_data_path      = "u10_1960-2020.nc",
  wind_v10_data_path      = "v10_1960-2020.nc",
  country_abbrev          = "FRA",
  mask_data_path          = "mask_FRA.nc",
  study_period            = c("1980-01-01", "2020-12-31"),
  reference_period        = c("1961-01-01", "1990-12-31"),
  granularity             = "year",
  percentile_high         = 95,
  percentile_low          = 5
)

result_dept <- calculate_aci(
  temperature_data_path   = "t2m_1960-2020.nc",
  precipitation_data_path = "tp_1960-2020.nc",
  wind_u10_data_path      = "u10_1960-2020.nc",
  wind_v10_data_path      = "v10_1960-2020.nc",
  country_abbrev          = "FRA",
  mask_data_path          = "mask_FRA.nc",
  study_period            = c("1980-01-01", "2020-12-31"),
  reference_period        = c("1961-01-01", "1990-12-31"),
  granularity             = "year",
  area                    = FALSE,
  admin_level             = 2,
  crs_metric              = 2154
)

## End(Not run)

Description

Calculate the half-day (day or night) temperature component

Usage

calculate_halfday_component(
  dataset,
  reference_period,
  part_of_day,
  extremum,
  percentile,
  above_thresholds,
  window_days = 5L
)

Arguments

Value

A list with data [lon × lat × months] and monthly time.

Description

Calculate the half-day (day or night) temperature component (terra version)

Usage

calculate_halfday_component_terra(
  r,
  reference_period,
  part_of_day,
  extremum,
  percentile,
  above_thresholds,
  mask_path = NULL,
  threshold = 0.8,
  cores = 1L,
  target_chunk_gb = 1,
  window_days = 5L
)

Arguments

Value

Same structure as calculate_halfday_component(): list(data, time, lon, lat).

Description

Computes the rolling window_size-day sum of precipitation, then takes the monthly maximum.

Usage

calculate_maximum_precipitation_over_window(
  dataset,
  var_name = "tp",
  window_size = 5L
)

Arguments

Value

A list with data [lon x lat x months] and time.

Description

Calculate maximum precipitation over a rolling window (terra version)

Usage

calculate_maximum_precipitation_over_window_terra(
  r,
  var_name = "tp",
  window_size = 5L,
  mask_path = NULL,
  threshold = 0.8,
  target_chunk_gb = 1
)

Arguments

Value

Same structure as calculate_maximum_precipitation_over_window().

Description

Uses a rolling window (in days) over the reference period, followed by a group-by-day-of-year percentile. Matches the CLIMDEX / ACI methodology (see ACI (2018), Appendix A of Garrido et al.): the threshold for a given calendar day is the n-th percentile of the daily TX/TN values falling within a window of window_days days (default 5, i.e. $\pm$2 days) around that calendar day, pooled across all years of the reference period.

Usage

calculate_percentiles(daily_dataset, n, reference_period, window_days = 5L)

Arguments

Details

Bug fix #1 (see NEWS): earlier versions of this function computed the rolling window/percentile on hourly t2m values (with window_size expressed in hours: 80 for day, 40 for night), then compared the resulting threshold to a daily extremum in .crossing_frequency(). Since the daily max (or min) of ~16 (or 8) hourly values almost always exceeds (or falls below) the 90th (or 10th) percentile of the marginal hourly distribution, this made the exceedance frequency during the reference period itself massively higher than the intended ~10 the reference frequency approximately by construction. calculate_percentiles() now takes the already-reduced daily extremum series (the output of temp_extremum()) as input, so the threshold and the tested variable are computed from the same daily quantity.

Bug fix #2 (see NEWS): even after bug fix #1, the threshold was still computed as a percentile of rolling percentiles: a window_days-wide rolling n-th percentile was first applied along the daily series, and a second n-th percentile was then taken of those already-extreme rolled values, grouped by day-of-year across reference years. This two-stage procedure is systematically biased: taking the n-th percentile of a set of local n-th percentiles pushes the final threshold further into the tail than a single, direct percentile would, so it under-counts the true exceedance frequency (observed ~5.6 even with bug fix #1 alone applied). The standard CLIMDEX/ACI methodology instead pools the raw daily values that fall within window_days of a given calendar day, across all years of the reference period, into a single sample, and takes one percentile of that pooled sample. calculate_percentiles() now follows this pooling approach (no more rolling window / zoo dependency for this function), which empirically recovers the intended ~10

Value

A numeric array [lon x lat x 366] (day-of-year 1-366).

Description

Equivalent of calculate_percentiles().

Usage

calculate_percentiles_terra(
  r_daily,
  n,
  reference_period,
  window_days = 5L,
  filename = "",
  cores = 1L
)

Arguments

Details

Bug fix #1 (see NEWS): this function used to take the raw hourly raster, filter it to day/night hours, and compute the rolling-window percentile directly on those hourly values (with window_size expressed in hours: 80 for day, 40 for night) – while calculate_halfday_component_terra() compared the resulting threshold to a daily max/min. Because the daily extremum of ~16 (or 8) hourly values almost always crosses the 90th (or 10th) percentile of the hourly distribution, this made the exceedance frequency during the reference period itself massively higher than the intended ~10 already-reduced daily extremum raster (the output of temp_extremum_terra()) as input, with window_days expressed in days (default 5), so the threshold and the tested variable are computed from the same daily quantity.

Bug fix #2 (see NEWS): even after bug fix #1, the threshold was still computed as a rolling-window n-th percentile (terra::roll()) followed by a second n-th percentile of those already-extreme rolled values, grouped by day-of-year (terra::tapp()). This two-stage procedure systematically under-counts the true exceedance frequency (observed ~5.6 applied) – see calculate_percentiles() (base-R engine) for the full explanation and numeric illustration, which applies identically here. calculate_percentiles_terra() now instead pools the RAW daily values within window_days of a given calendar day, across ALL years of the reference period, and takes a SINGLE terra::app() quantile of that pooled sample per cell, per calendar day – the standard CLIMDEX/ACI methodology. This also removes the terra::roll()/terra::tapp() dependency for this function entirely (simpler, and no longer sensitive to terra::roll()'s argument names changing across terra versions).

Performance AND memory note: looping over 366 calendar days and calling terra::app() on a (typically ~100-200 layer) subset of the reference period for each is comparable in total workload to the single whole-series terra::roll() call used before bug fix #2, but touches far fewer layers per call. A single such call over a whole country's grid could still, in principle, use substantial memory for a large grid – as a precaution the raster is still ALWAYS split into small spatial tiles sized to a conservative, fixed memory target (see target_tile_gb in .calculate_percentiles_terra_tiled()), regardless of cores – cores only controls how many of those already-memory-safe tiles run concurrently (default 1: one at a time).

Value

A terra::SpatRaster with 366 layers (day-of-year 1-366).

Description

Calculate monthly wind exceedance frequency

Usage

calculate_period_wind_exceedance_frequency(
  u10_dataset,
  v10_dataset,
  reference_period
)

Arguments

Value

A list with data [lon x lat x periods] and time.

Description

Drop-in, memory-safe replacement for calculate_period_wind_exceedance_frequency(). Only the initial wind_power_terra() call differs from the base-R version; the threshold, crossing, and monthly-aggregation steps reuse the exact same helpers, since they already operate on small, daily-resolution data.

Usage

calculate_period_wind_exceedance_frequency_terra(
  u10_r,
  v10_r,
  reference_period,
  mask_path = NULL,
  threshold = 0.8,
  target_chunk_gb = 1
)

Arguments

Value

Same as calculate_period_wind_exceedance_frequency().

Description

Computes a rolling (sliding-window) sum along the time dimension of a 3-D climate array.

Usage

calculate_rolling_sum(dataset, var_name, window_size)

Arguments

Value

A list with the same structure as dataset but with data replaced by the rolling sum (leading incomplete windows are NA).

Description

Wrapper around wf_set_key for the new CDS API (post-October 2024). The token is stored securely in the system keyring and never needs to appear in your scripts again.

Usage

cds_set_key(token)

Arguments

Details

To obtain your token:

1. Create a free account at https://cds.climate.copernicus.eu.

2. Go to your profile → Personal Access Token.

3. Copy the token (a UUID string) and pass it to this function once.

Value

Invisibly NULL. The token is saved in the system keyring.

Examples

## Not run: 
cds_set_key("xxxxxxxx-xxxx-xxxx-xxxx-xxxxxxxxxxxx")

## End(Not run)

Description

Low-level helpers shared by all ACI component functions.

Description

Memory-safe alternatives to load_netcdf()/apply_mask() and to the hourly-resolution steps of the temperature, precipitation and wind pipelines (temp_extremum(), calculate_percentiles(), the daily resampling inside wind_power()).

Why this file exists: for a whole country at hourly resolution over 40+ years, ncdf4::ncvar_get() materialises the full [lon x lat x time] array in RAM in one shot (tens of GB). terra reads NetCDF lazily (via GDAL) and processes computations block-by-block, writing results to disk instead of accumulating them in memory.

Only the hourly-scale steps are ported here. Once data has been reduced to DAILY resolution (via resample_daily_terra(), temp_extremum_terra(), ...), the resulting object is small enough (~40 years x 365 days) to convert back to a plain array with .spatraster_to_list() and hand off, unchanged, to the rest of the existing (already tested) pipeline – resample_monthly(), standardize_metric(), .compute_aci_grid(), etc.

Description

Calculate days with wind power above the 90th-percentile threshold

Usage

days_above_wind_thresholds(u10_dataset, v10_dataset, reference_period)

Arguments

Value

A list with binary data [lon x lat x days] and time.

Description

Functions to download ERA5 climate variables and country mask from the Copernicus Climate Data Store (CDS) via the ecmwfr package (v2.0+, new CDS API — Personal Access Token only, no UID required).

Description

Downloads hourly ERA5 single-level data for one variable and one or several years, for a given geographical bounding box. Files are downloaded year by year and then optionally merged into a single NetCDF using the ncdf4 package.

Usage

download_era5(
  variable = c("t2m", "tp", "u10", "v10"),
  years,
  area = c(51.5, -5.5, 41, 10),
  country_abbrev = NULL,
  dest_dir = NULL,
  merge = TRUE,
  overwrite = FALSE
)

Arguments

Details

The CDS Personal Access Token must be stored beforehand with cds_set_key.

Value

Invisibly, the path to the final NetCDF file (merged if merge = TRUE, otherwise the directory containing yearly files).

Examples

## Not run: 
# Store your token once
cds_set_key("xxxxxxxx-xxxx-xxxx-xxxx-xxxxxxxxxxxx")

# Automatic path: <tempdir()>/xaci_era5/FRA/tp_1961_1990.nc
download_era5(
  variable       = "tp",
  years          = 1961:1990,
  area           = c(51.5, -5.5, 41.0, 10.0),
  country_abbrev = "FRA"
)

# Or explicit path:
download_era5(
  variable = "tp",
  years    = 1961:1990,
  area     = c(51.5, -5.5, 41.0, 10.0),
  dest_dir = "my/custom/path"
)

## End(Not run)

Description

Convenience wrapper that calls download_era5 four times (one per variable: t2m, tp, u10, v10).

Usage

download_era5_all(
  years,
  area = c(51.5, -5.5, 41, 10),
  country_abbrev = NULL,
  dest_dir = NULL,
  merge = TRUE,
  overwrite = FALSE
)

Arguments

Value

Invisibly, a named character vector of output paths (one per variable).

Examples

## Not run: 
cds_set_key("xxxxxxxx-xxxx-xxxx-xxxx-xxxxxxxxxxxx")

# Produces (by default): <tempdir()>/xaci_era5/FRA/t2m_1961_1990.nc, tp_..., u10_..., v10_...
download_era5_all(
  years          = 1961:1990,
  area           = c(51.5, -5.5, 41.0, 10.0),
  country_abbrev = "FRA"
)

## End(Not run)

Description

Downloads the ERA5 land-sea mask for the given bounding box, then rescales it to [0, 1] and saves it as mask_<country_abbrev>.nc in dest_dir. The mask variable is renamed to country so it is directly compatible with apply_mask.

Usage

download_mask(
  country_abbrev = "FRA",
  area = c(51.5, -5.5, 41, 10),
  dest_dir = NULL,
  overwrite = FALSE
)

Arguments

Value

Invisibly, the path to the mask NetCDF file.

Examples

## Not run: 
cds_set_key("xxxxxxxx-xxxx-xxxx-xxxx-xxxxxxxxxxxx")

# Produces (by default): <tempdir()>/xaci_era5/FRA/mask_FRA.nc
download_mask(
  country_abbrev = "FRA",
  area      = c(51.5, -5.5, 41.0, 10.0)
)

## End(Not run)

Description

Computes the drought component of the Actuarial Climate Index via maximum consecutive dry days (CDD).

Description

Computes standardised consecutive dry days.

Usage

drought_component(
  precipitation_data_path,
  country_abbrev,
  reference_period,
  study_period,
  mask_path = NULL,
  area = FALSE,
  admin_level = NULL,
  admin_mask = NULL,
  crs_metric = 4326,
  computed_components = FALSE,
  save = FALSE,
  save_dir = NULL,
  load_dir = NULL
)

Arguments

Value

Named numeric vector (area = TRUE), standardised list (area = FALSE), or data.frame per admin unit.

Description

Drop-in, memory-safe replacement for drought_component().

Usage

drought_component_terra(
  precipitation_data_path,
  country_abbrev,
  reference_period,
  study_period,
  mask_path = NULL,
  area = FALSE,
  admin_level = NULL,
  admin_mask = NULL,
  crs_metric = 4326,
  computed_components = FALSE,
  save = FALSE,
  save_dir = NULL,
  load_dir = NULL,
  target_chunk_gb = 1
)

Arguments

Value

Same as drought_component().

Description

For month $j \in 1..11$ of year $k$, the monthly CDD interpolates between the ANNUAL values of year $k-1$ and year $k$ (each annual value being anchored at December 31 of its own year):

$CDD_j(k) = \frac{12-j}{12} \cdot CDD_{k-1} + \frac{j}{12} \cdot CDD_k$

and $CDD_{12}(k) = CDD_k$ exactly (December of year $k$ coincides with that year's own annual anchor point). This matches Eq. (A.6) of Garrido et al. The first year, having no $k-1$ to interpolate from, repeats its own annual value for all 12 months.

Usage

drought_interpolate(cdd_annual)

Arguments

Details

Bug fix (see NEWS): earlier versions of this function computed this same interpolation formula correctly, but assigned the result to the WRONG calendar year – labelling month $j$ of year $\text{years}[k]$ as w1 * data[,,k] + w2 * data[,,k+1], which is actually month $j$ of year $\text{years}[k+1]$ per Eq. (A.6) (whose December value must equal that year's OWN annual value, not the following year's). Verified numerically: with annual CDD values 10/20/30/40 for 2000-2003, December 2000 came out as 20 (year 2001's annual value) instead of 10 (year 2000's own value). The interpolation arithmetic itself was correct; only the year label was off by one. As a consequence, the "repeat own value for all 12 months" edge-case handling now applies to the FIRST year (which has no preceding year to interpolate from) instead of the last (which, under the corrected labelling, has no such issue: Eq. A.6 only ever needs year $k-1$ to compute year $k$'s months, so the last year is computed normally).

Value

A list with data [lon × lat × months] and monthly time.

Description

Memory-safe equivalents of max_consecutive_dry_days() and drought_component(). Only the hourly-to-daily reduction runs through terra; drought_interpolate() and the consecutive-dry-day scan already operate on small, annual/daily-resolution data, so they are reused as-is via the shared helper .max_consecutive_dry_days_from_daily() (defined in drought.R) – no logic is duplicated.

Description

This function is a constructor. It reads a NetCDF variable and, if mask_path is provided, zeroes out (sets to NA) all grid cells where the mask value is below threshold.

Usage

load_component(data_path, var_name, mask_path = NULL, threshold = 0.8)

Arguments

Value

A list with elements data (3-D array [lon × lat × time]), lon, lat, time (POSIXct vector), and var_name.

Description

Drop-in, memory-safe replacement for load_component() intended for full-resolution, grid-cell-level ("area = FALSE") workflows over long historical periods.

Usage

load_component_terra(
  data_path,
  var_name,
  mask_path = NULL,
  threshold = 0.8,
  chunk_size = 20000
)

Arguments

Value

A terra::SpatRaster.

Description

Equivalent of load_netcdf(), but never materialises pixel values in RAM: only metadata (dimensions, CRS, time) is read up front.

Usage

load_netcdf_terra(path, var_name)

Arguments

Value

A terra::SpatRaster with one layer per time step and terra::time() attached.

Description

Load the bundled PSMSL station metadata

Usage

load_psmsl_data()

Value

A data.frame of PSMSL tide-gauge station information.

Description

A dry day is defined as total daily precipitation < 1 mm. The function returns the annual maximum number of consecutive dry days for each grid cell.

Usage

max_consecutive_dry_days(dataset)

Arguments

Value

A list with data [lon × lat × years] and annual time.

Description

Calculate maximum consecutive dry days (terra version)

Usage

max_consecutive_dry_days_terra(
  r,
  mask_path = NULL,
  threshold = 0.8,
  target_chunk_gb = 1
)

Arguments

Value

Same structure as max_consecutive_dry_days().

Description

Produces either a bar chart (one panel per component, type = "bar") or a stacked area chart (type = "stacked") showing the six standardised components that make up the ACI.

Usage

plot_aci_components(
  aci_df,
  type = c("stacked", "line", "bar"),
  components = NULL,
  title = "ACI Component Contributions",
  subtitle = NULL
)

Arguments

Value

A ggplot object.

Examples

# Toy monthly ACI data.frame (see plot_aci_timeseries() for details)
n     <- 60L
dates <- format(seq(as.Date("2015-01-01"), by = "month", length.out = n),
                "%Y-%m")
aci_df <- data.frame(
  t90           = 20 + stats::rnorm(n, sd = 2),
  t10           = 5  + stats::rnorm(n, sd = 2),
  precipitation = stats::rnorm(n, sd = 1),
  drought       = stats::rnorm(n, sd = 1),
  wind          = stats::rnorm(n, sd = 1),
  sealevel      = stats::rnorm(n, sd = 1),
  row.names     = dates
)
aci_df$ACI <- with(aci_df, (t90 - t10 + precipitation + drought +
                               wind + sealevel) / 5)

plot_aci_components(aci_df, type = "bar")
plot_aci_components(aci_df, type = "stacked", components = c("t90", "t10"))

Description

Arranges the four core plots (time series, component facets, boxplots and density) into a single patchwork figure. Requires the patchwork package to be installed.

Usage

plot_aci_dashboard(aci_df, title = "ACI Summary Dashboard")

Arguments

Value

A patchwork object (inherits from ggplot).

Examples

# Toy monthly ACI data.frame (see plot_aci_timeseries() for details)
n     <- 60L
dates <- format(seq(as.Date("2015-01-01"), by = "month", length.out = n),
                "%Y-%m")
aci_df <- data.frame(
  t90           = 20 + stats::rnorm(n, sd = 2),
  t10           = 5  + stats::rnorm(n, sd = 2),
  precipitation = stats::rnorm(n, sd = 1),
  drought       = stats::rnorm(n, sd = 1),
  wind          = stats::rnorm(n, sd = 1),
  sealevel      = stats::rnorm(n, sd = 1),
  row.names     = dates
)
aci_df$ACI <- with(aci_df, (t90 - t10 + precipitation + drought +
                               wind + sealevel) / 5)
plot_aci_dashboard(aci_df)

Description

Displays the statistical distribution of each standardised component (and optionally the global ACI) as boxplots (type = "boxplot"), violin plots (type = "violin"), or density curves (type = "density").

Usage

plot_aci_distribution(
  aci_df,
  components = NULL,
  include_aci = FALSE,
  type = c("boxplot", "violin", "density"),
  title = "ACI Component Distributions",
  subtitle = NULL
)

Arguments

Value

A ggplot object.

Examples

# Toy monthly ACI data.frame (see plot_aci_timeseries() for details)
n     <- 60L
dates <- format(seq(as.Date("2015-01-01"), by = "month", length.out = n),
                "%Y-%m")
aci_df <- data.frame(
  t90           = 20 + stats::rnorm(n, sd = 2),
  t10           = 5  + stats::rnorm(n, sd = 2),
  precipitation = stats::rnorm(n, sd = 1),
  drought       = stats::rnorm(n, sd = 1),
  wind          = stats::rnorm(n, sd = 1),
  sealevel      = stats::rnorm(n, sd = 1),
  row.names     = dates
)
aci_df$ACI <- with(aci_df, (t90 - t10 + precipitation + drought +
                               wind + sealevel) / 5)
plot_aci_distribution(aci_df, type = "violin")
plot_aci_distribution(aci_df, type = "density")

Description

Handles three types of input, at any spatial aggregation level produced by calculate_aci() (grid-cell, admin level 1, admin level 2, ...):

• Bare array (e.g. grid$t90, grid$ACI): a single component array as stored in the grid-cell list returned by calculate_aci(area = FALSE). Self-describing via its attached lon/lat/time/country_abbrev attributes, so it can be passed on its own. Renders a raster map. Since ACI and its components share the exact same array structure, this works identically whether variable is a component or the ACI itself.

• Grid-cell list (the full output of calculate_aci(area = FALSE), or of load_component()): variable selects which field of the list to map. Renders a raster map.

• Administrative data.frame (output of calculate_aci(admin_level = N)): renders a choropleth map. country_abbrev/admin_level/crs_metric are read from the data.frame's attributes when not supplied explicitly.

Usage

plot_aci_map(
  data,
  variable = "ACI",
  time_index = "mean",
  country_abbrev = NULL,
  admin_level = NULL,
  crs_metric = NULL,
  var_label = "Standardised anomaly",
  title = "ACI Spatial Distribution",
  palette = "RdBu",
  reverse = TRUE,
  n_breaks = 9,
  borders = TRUE
)

Arguments

Value

A ggplot object.

Offline vs. online examples

The first two examples below use a small synthetic 10x8x5 (lon x lat x time) grid and run offline (borders = FALSE). Overlaying administrative borders, or plotting an admin-level choropleth, additionally requires downloading boundary polygons via geodata::gadm() and therefore a working internet connection; the corresponding examples are wrapped in dontrun and not run automatically by R CMD check.

Examples

lon <- seq(-5, 8, length.out = 10)
lat <- seq(42, 51, length.out = 8)
arr <- array(stats::rnorm(10 * 8 * 5), dim = c(10, 8, 5))
arr <- structure(arr, lon = lon, lat = lat,
                  time = seq_len(5), country_abbrev = "FRA")

# Bare array, self-describing thanks to its attributes
plot_aci_map(arr, time_index = "mean", borders = FALSE)

# Equivalent, via the parent list and `variable`
grid <- list(t90 = arr, lon = lon, lat = lat, time = seq_len(5))
plot_aci_map(grid, variable = "t90", time_index = 1, borders = FALSE)

## Not run: 
plot_aci_map(arr, time_index = "mean", country_abbrev = "FRA")

admin <- calculate_aci(
  temperature_data_path   = "t2m_1960-2020.nc",
  precipitation_data_path = "tp_1960-2020.nc",
  wind_u10_data_path      = "u10_1960-2020.nc",
  wind_v10_data_path      = "v10_1960-2020.nc",
  country_abbrev          = "FRA",
  mask_data_path          = "mask_FRA.nc",
  study_period            = c("1961-01-01", "2020-12-31"),
  reference_period        = c("1961-01-01", "1990-12-31"),
  granularity             = "year",
  area                    = FALSE,
  admin_level             = 1
)
plot_aci_map(admin, variable = "ACI", time_index = 1)

## End(Not run)

Description

Convenience wrapper around plot_aci_map that computes the mean over a user-defined date range rather than a single time index.

Usage

plot_aci_map_mean(
  data,
  variable = "ACI",
  period = NULL,
  country_abbrev = NULL,
  admin_level = NULL,
  crs_metric = NULL,
  var_label = "Standardised anomaly",
  title = NULL,
  palette = "RdBu",
  reverse = TRUE,
  n_breaks = 9,
  borders = TRUE
)

Arguments

Value

A ggplot object.

Description

Draws a line chart of the ACI values returned by calculate_aci, with an optional smoothed trend and a horizontal zero reference line.

Usage

plot_aci_timeseries(
  aci_df,
  smooth = TRUE,
  span = 0.2,
  title = "Actuarial Climate Index (ACI)",
  colour = "#1F77B4",
  fill_area = TRUE
)

Arguments

Value

A ggplot object.

Examples

# Toy monthly ACI data.frame (not from calculate_aci(), just for
# illustration), with the columns and "YYYY-MM" row names that
# calculate_aci(area = TRUE) would return.
n     <- 60L
dates <- format(seq(as.Date("2015-01-01"), by = "month", length.out = n),
                "%Y-%m")
aci_df <- data.frame(
  t90           = 20 + stats::rnorm(n, sd = 2),
  t10           = 5  + stats::rnorm(n, sd = 2),
  precipitation = stats::rnorm(n, sd = 1),
  drought       = stats::rnorm(n, sd = 1),
  wind          = stats::rnorm(n, sd = 1),
  sealevel      = stats::rnorm(n, sd = 1),
  row.names     = dates
)
aci_df$ACI <- with(aci_df, (t90 - t10 + precipitation + drought +
                               wind + sealevel) / 5)
plot_aci_timeseries(aci_df)

Description

Computes the precipitation component of the Actuarial Climate Index.

Description

Computes the standardised anomaly of maximum monthly precipitation over a rolling 5-day window.

Usage

precipitation_component(
  precipitation_data_path,
  country_abbrev,
  reference_period,
  study_period,
  mask_path = NULL,
  var_name = "tp",
  window_size = 5L,
  area = FALSE,
  admin_level = NULL,
  admin_mask = NULL,
  crs_metric = 4326,
  computed_components = FALSE,
  save = FALSE,
  save_dir = NULL,
  load_dir = NULL
)

Arguments

Value

Named numeric vector, standardised list, or data.frame per admin unit.

Description

Drop-in, memory-safe replacement for precipitation_component().

Usage

precipitation_component_terra(
  precipitation_data_path,
  country_abbrev,
  reference_period,
  study_period,
  mask_path = NULL,
  var_name = "tp",
  window_size = 5L,
  area = FALSE,
  admin_level = NULL,
  admin_mask = NULL,
  crs_metric = 4326,
  computed_components = FALSE,
  save = FALSE,
  save_dir = NULL,
  load_dir = NULL,
  target_chunk_gb = 1
)

Arguments

Value

Same as precipitation_component().

Description

Memory-safe equivalents of calculate_maximum_precipitation_over_window() and precipitation_component(). Only the hourly-to-daily reduction runs through terra; the rolling-sum and monthly-max logic already operates on small, daily-resolution data, reused as-is via the shared helper .max_precipitation_from_daily() (defined in precipitation.R) – no logic is duplicated.

Description

If admin_mask is NULL, averages over all non-NA cells (national mean). Otherwise computes a weighted mean per administrative unit using surface-fraction weights.

Usage

reduce_dataarray_to_dataframe(metric, column_name = "value", admin_mask = NULL)

Arguments

Value

A data.frame with one column per administrative unit (or one column for the national mean), indexed by month-start dates.

Description

Reads the bundled psmsl_data.csv to identify stations for the given country abbreviation, downloads the corresponding data files from the PSMSL website, and stores them locally.

Usage

request_sealevel_data(country_abbrev, dest_dir = NULL)

Arguments

Value

Invisibly, the path to the destination directory.

Description

Equivalent of resample_daily(). Groups layers by calendar day and applies fun via terra::tapp().

Usage

resample_daily_terra(r, fun = "mean", filename = "", target_chunk_gb = 1)

Arguments

Details

Memory note : pour agreger par jour, terra::tapp() doit lire, PAR BLOC SPATIAL, la serie temporelle COMPLETE de r (readValues() recupere nrows x ncol x nlyr(r) valeurs d'un coup, puis matrix() en fait une copie contigue). Sur une grille pourtant petite (quelques milliers de cellules), nlyr(r) horaire sur plusieurs decennies (des centaines de milliers de couches) suffit a lui seul a depasser la limite memoire d'un processus R (constate empiriquement : plafond de 16 Go de R sous macOS, "vector memory limit ... reached", pour ~65 ans horaires sur la France entiere) – MEME quand le resultat final (agrege par jour) est minuscule. C'est exactement la meme classe de probleme que celle documentee pour terra::roll() dans .calculate_percentiles_terra_tiled() : le pic memoire depend de la taille de L'ENTREE lue par bloc, pas de la sortie. On applique donc ici le meme principe, en decoupant TEMPORELLEMENT (et non spatialement, puisque c'est le nombre de couches qui explose ici, pas le nombre de cellules) : r est scinde en blocs de jours CONSECUTIFS (donc de couches horaires contigues – l'ordre chronologique de r est requis), chaque bloc est agrege independamment et ecrit sur disque, puis les blocs journaliers (petits) sont recombines. Le resultat est identique, bloc par bloc ou en un seul appel : chaque jour est agrege a partir de ses propres couches horaires, jamais a cheval sur deux blocs.

Precision numerique : les blocs intermediaires sont ecrits sur disque en datatype = "FLT8S" (double precision, 64 bits) – explicitement force pour eviter que terra ne bascule par defaut sur du float 32 bits, ce qui romprait l'invariance du resultat par rapport a target_chunk_gb (chemin non-decoupe = calcul garde en RAM, en double precision R). Le resultat est ainsi, a l'arithmetique flottante pres (memes operations, meme ordre par groupe), rigoureusement identique quel que soit le nombre de blocs.

Value

A terra::SpatRaster with one layer per day.

Description

Processes PSMSL tide-gauge data and computes the sea-level component of the Actuarial Climate Index.

Description

Replace PSMSL sentinel values (-99999) with NA

Usage

sealevel_clean_data(df)

Arguments

Value

Cleaned data.frame.

Description

Calculate the sea-level component of the ACI

Usage

sealevel_component(
  country_abbrev,
  study_period,
  reference_period,
  mask_path = NULL,
  area = TRUE,
  max_dist_km = 500,
  sealevel_dir = NULL,
  admin_level = NULL,
  admin_assignment = NULL,
  crs_metric = 4326,
  computed_components = FALSE,
  save = FALSE,
  save_dir = NULL,
  load_dir = NULL
)

Arguments

Value

If area = FALSE: a list with a [lon x lat x t] array, lon, lat, time — same structure as other grid-cell components. If area = TRUE: a data.frame of station anomalies (national or per admin unit).

Description

Compute monthly reference statistics for sea-level data

Usage

sealevel_compute_monthly_stats(df, reference_period, stats)

Arguments

Value

A numeric matrix, 12 rows (calendar months "1"-"12" as row names) x one column per station (station names as returned by colnames(df)). Each station is standardised against its own monthly reference statistics, independently of the other stations – consistent with how the other ACI components (ERA5 grid cells, administrative units) are each standardised against their own reference, not a value pooled across the whole spatial domain.

Description

PSMSL encodes dates as YYYY.fraction where the fraction identifies the month.

Usage

sealevel_correct_date_format(df)

Arguments

Value

The same data.frame with Date row names ("YYYY-MM-01"), rows with unrecognised dates removed.

Description

Load sea-level txt files from a directory

Usage

sealevel_load_data(directory)

Arguments

Value

A data.frame with one column per station and a numeric date as row names.

Description

Full sea-level processing pipeline

Usage

sealevel_process(directory, study_period, reference_period)

Arguments

Value

A named list with:

data

Standardised data.frame of anomalies [time x stations], row names "YYYY-MM-DD". Stations listed in excluded_stations are entirely absent from this data.frame (all-NA columns are dropped), not merely NA.

coords

A data.frame with columns station_id, lon, lat, one row per station present in data (excluded stations are also dropped here).

excluded_stations

Character vector of station column names (possibly empty) that have ZERO observations in reference_period for at least one calendar month, and are therefore impossible to standardise against Eq. A.12 (no reference mean is definable for that station/month) – see the "Bug fix" note below. Use this to know which stations (e.g. ones installed after 1990) silently contribute nothing, so this isn't found out by surprise downstream.

Bug fix (see NEWS)

Stations with literally no data during reference_period (e.g. a tide gauge installed in 2018, for a 1961-1990 reference period) used to get a reference mean of NaN (not NA), which propagated through sealevel_standardize_data() to make the entire standardised series of that station NaN – including during years where the station DOES have perfectly good data. This had no visible effect only because reduce_sealevel_over_region()'s rowMeans(..., na.rm = TRUE) happens to also drop NaN (R treats NaN as a kind of NA), so the station was silently and permanently excluded from every national/regional average it could have contributed to, with no warning. For the French stations listed in Garrido et al.'s Table B.2, 12 of the 40 stations (30%) have zero overlap with 1961-1990 and were affected. This function now (a) produces clean NA instead of NaN for this case, and more importantly (b) explicitly detects and reports these stations via warning() and the excluded_stations return value, instead of relying on an implicit, easy-to-miss floating-point coincidence. Excluding such stations is not itself a choice this function makes – there is no principled reference baseline to compute for a station that did not exist during the reference period – but that exclusion is now visible rather than silent.

Description

Standardise sea-level data over the study period

Usage

sealevel_standardize_data(df, monthly_means, monthly_std_devs, study_period)

Arguments

Value

A data.frame of standardised anomalies for the study period, with rows containing all-NA removed. Each station's values are standardised against its own monthly reference (see sealevel_compute_monthly_stats()), not a value pooled across stations.

Description

For each calendar month, computes the mean and standard deviation over the reference period, then returns (x - mean) / sd. Optionally averages over the spatial dimensions first.

Usage

standardize_metric(metric, reference_period, area = FALSE)

Arguments

Value

If area = TRUE: a named numeric vector (names = dates). Otherwise: a list with the same structure as metric but standardised values.

Description

Compute daily temperature extremum (min or max) for day or night hours

Usage

temp_extremum(dataset, extremum, period)

Arguments

Value

A list with data [lon × lat × days] and daily time.

Description

Equivalent of temp_extremum(). Filters layers by hour-of-day (a cheap, lazy operation on a SpatRaster – no data is read), then reduces to daily resolution via resample_daily_terra().

Usage

temp_extremum_terra(r, extremum, period, filename = "", target_chunk_gb = 1)

Arguments

Value

A terra::SpatRaster with one layer per day.

Description

Computes the temperature component of the Actuarial Climate Index.

Description

Combines day and night half-day components (equal weighting), then standardises relative to the reference period.

Usage

temperature_component(
  temperature_data_path,
  country_abbrev,
  reference_period,
  study_period,
  mask_path = NULL,
  percentile = 90,
  extremum = "max",
  above_thresholds = TRUE,
  area = FALSE,
  admin_level = NULL,
  admin_mask = NULL,
  crs_metric = 4326,
  computed_components = FALSE,
  save = FALSE,
  save_dir = NULL,
  load_dir = NULL,
  window_days = 5L
)

Arguments

Value

Named numeric vector, standardised list, or data.frame per admin unit.

Description

Drop-in, memory-safe replacement for temperature_component(), intended for full-resolution, grid-cell-level ("area = FALSE") workflows over long historical periods (40+ years hourly).

Usage

temperature_component_terra(
  temperature_data_path,
  country_abbrev,
  reference_period,
  study_period,
  mask_path = NULL,
  percentile = 90,
  extremum = "max",
  above_thresholds = TRUE,
  area = FALSE,
  admin_level = NULL,
  admin_mask = NULL,
  crs_metric = 4326,
  cores = 1L,
  computed_components = FALSE,
  save = FALSE,
  save_dir = NULL,
  load_dir = NULL,
  target_chunk_gb = 1,
  window_days = 5L
)

Arguments

Value

Same as temperature_component(): a standardized metric list, or a data frame if admin_mask/admin_level is provided.

Description

Memory-safe equivalents of calculate_halfday_component() and temperature_component(), for full-resolution, grid-cell-level workflows over long historical periods (see R/component_terra.R for the underlying loading/reduction primitives).

Only the hourly-scale steps (temp_extremum_terra(), calculate_percentiles_terra()) run through terra. Once reduced to daily/monthly resolution, the exact same .crossing_frequency() helper used by the base-R pipeline takes over – no logic is duplicated.

Description

Helper functions to manipulate and merge climate data.

Description

A collection of ggplot2-based plotting functions to visualize the Actuarial Climate Index (ACI) and its components. Four families of plots are provided:

1. Time series of the global ACI (plot_aci_timeseries)

2. Component decomposition stacked / faceted chart (plot_aci_components)

3. Geographic maps from a NetCDF-derived array (plot_aci_map)

4. Distribution plots - boxplots and density ridges (plot_aci_distribution)

Description

Computes the wind component of the Actuarial Climate Index.

Description

Calculate the wind component of the ACI

Usage

wind_component(
  wind_u10_data_path,
  wind_v10_data_path,
  country_abbrev,
  reference_period,
  study_period,
  mask_path = NULL,
  area = FALSE,
  admin_level = NULL,
  admin_mask = NULL,
  crs_metric = 4326,
  computed_components = FALSE,
  save = FALSE,
  save_dir = NULL,
  load_dir = NULL
)

Arguments

Value

Named numeric vector, standardised list, or data.frame per admin unit.

Description

Drop-in, memory-safe replacement for wind_component(), intended for full-resolution, grid-cell-level ("area = FALSE") workflows over long historical periods (40+ years hourly).

Usage

wind_component_terra(
  wind_u10_data_path,
  wind_v10_data_path,
  country_abbrev,
  reference_period,
  study_period,
  mask_path = NULL,
  area = FALSE,
  admin_level = NULL,
  admin_mask = NULL,
  crs_metric = 4326,
  computed_components = FALSE,
  save = FALSE,
  save_dir = NULL,
  load_dir = NULL,
  target_chunk_gb = 1
)

Arguments

Value

Same as wind_component().

Description

Wind speed = sqrt(u10^2 + v10^2), wind power = 0.5 * rho * ws^3.

Usage

wind_power(u10_dataset, v10_dataset, reference_period = NULL)

Arguments

Details

Bug fix (see NEWS): earlier versions averaged u10 and v10 to daily means SEPARATELY first, then took the magnitude of that daily-mean vector (sqrt(mean(u)^2 + mean(v)^2)) as the "daily mean wind speed". This "vector averaging" systematically UNDER-estimates the true daily mean wind speed whenever wind direction varies over the day, because opposing (or merely differing) hourly components partially cancel out before the magnitude is taken – e.g. a wind blowing at a constant 8 m/s whose direction oscillates by only +/-60 degrees over the day averages out to a vector magnitude of ~6 m/s, a ~58% error on the resulting wind power (0.5*rho*w^3), verified numerically. The standard meteorological convention for a "daily mean wind speed" is instead to average the INSTANTANEOUS scalar speed sqrt(u^2+v^2) computed at each timestep, which is always positive and does not suffer from this cancellation; that is what this function now does.

Value

A list with data [lon x lat x days] (wind power in W/m²) and daily time.

Description

Terra-based equivalent of wind_power(): computes daily mean u/v via resample_daily_terra() (block-by-block, never loading the full hourly cube in RAM), then wind speed/power arithmetic directly on the (now daily, small) SpatRasters before converting back to the package's list format.

Usage

wind_power_terra(
  u10_r,
  v10_r,
  reference_period = NULL,
  mask_path = NULL,
  threshold = 0.8,
  target_chunk_gb = 1
)

Arguments

Value

Same structure as wind_power(): list(data, time, lon, lat).

Description

Memory-safe equivalents of wind_power() and wind_component(). Only wind_power_terra() touches raw hourly data (via terra); wind_thresholds()'s internal logic, threshold-crossing, and monthly aggregation already operate on DAILY resolution data (small even for 40+ years), so they are reused as-is via the shared helpers .wind_thresholds_from_wp(), .days_above_from_wp(), .monthly_frequency_from_binary() defined in wind.R – no logic is duplicated.

Description

Calculate wind power thresholds (90th percentile per day-of-year)

Usage

wind_thresholds(u10_dataset, v10_dataset, reference_period)

Arguments

Value

A list with data [lon x lat x all_days] giving the threshold for each day of the full series, and time.