Skip to contents

Model-ready meteorological forcing from gridded reanalysis and climate projections, bias-corrected against local observations and temporally disaggregated. metscale covers the meteorological side of process-model forcing without depending on any model: extract, correct, disaggregate, and a set of clean meteorological conversions. It integrates with the Aquatic Ecosystem Model Ensemble (AEME) when that is installed, but does not require it.

What it does

Task Function
Download ERA5 forcing: global daily point series (ISIMIP3a, no key), hourly GRIB from the Copernicus Data Store, or aggregate downloaded netCDF / GRIB to daily download_era5_isimip_point(), download_era5_cds() / read_era5_grib_point(), convert_era5_netcdf()
Extract hourly ERA5-Land at a point or averaged over a lake polygon, in standard units and time zone extract_era5_hourly_met(), extract_era5_lake_met()
Standardise a table of measured meteorology (names, units, timezone, resampling) prepare_obs_met(), standardise_met(), guess_met_vars()
Translate a met table between the AEME MET_* scheme and CF / CMIP short names and units (tas, pr, sfcWind, …) met_to_cf(), cf_to_met()
Fit a bias correction ERA5 -> observations (monthly scaling, variance scaling, regression, empirical / trend-preserving quantile mapping) with leave-one-year-out cross-validation fit_met_bias_correction()
Apply it to the full record and aggregate to daily apply_met_bias_correction(), met_to_daily()
Build a bias-corrected daily baseline for a climate-scenario (delta-change) workflow bias_correct_daily_baseline(), ?scenario_workflow
Read daily CMIP6 / CCAM projections at a point (model-calendar aware) for the delta-change step extract_climate_point()
Disaggregate daily meteorology to hourly / 3-hourly (method of fragments or mean diurnal cycle) disaggregate_met_to_hourly(), build_diurnal_climatology()
Fill a minimal met set out to everything a lake model needs expand_met()
Solar geometry, clear-sky shortwave, hourly shortwave from daily solar_zenith_angle(), clear_sky_swr(), estimate_hourly_swr()
Adjust wind speed between measurement heights (e.g. a 2 m buoy anemometer to the 10 m reanalysis / lake-model convention) wind_at_height(), met_wind_at_height()
Humidity, wind-vector, pressure and longwave conversions rh_to_dewpoint(), uv2ds() / ds2uv(), station_pressure(), calc_in_lwr(), calc_humidity_vars(), calc_cc()

Standard variables

Columns follow the AEME MET_* scheme (met_vars()): MET_radswd W/m2, MET_radlwd W/m2, MET_cldcvr fraction, MET_tmpair degC, MET_humrel %, MET_tmpdew degC, MET_prvapr hPa, MET_prsttn / MET_prmslp Pa, MET_wndspd m/s, MET_wnddir degrees (from), MET_wnduvu / MET_wnduvv m/s, MET_pprain / MET_ppsnow mm per timestep.

The default time zone is "UTC", matching the clock ERA5, ERA5-Land and CMIP6 are published on, so the extractors apply no hidden shift. Every tz argument resolves the same way: an explicit value first, then the tz attribute carried on the input, then UTC. Pass a fixed-offset zone such as "Etc/GMT-12" (NZST, UTC+12, no daylight saving) to work in local time while keeping a gap-free regular sub-daily series.

met_to_cf() maps this scheme onto the CF / CMIP short names and units (MET_tmpair degC → tas K, MET_pprain mm/step → pr kg m-2 s-1, MET_wndspd → sfcWind, MET_radswd → rsds, …) for handing data to tools that expect that vocabulary; cf_to_met() is the inverse.

Relationship to AEME

metscale holds the canonical implementation of the meteorology helpers that previously lived inside AEME (expand_met, standardise_met, calc_cc, calc_in_lwr, calc_swr / estimate_hourly_swr, the humidity, wind and pressure conversions). Three long-standing defects were fixed on the way across:

  • Time zone. expand_met() / calc_cc() / estimate_hourly_swr() take an explicit tz. The AEME versions force the session time zone to UTC and read the wall clock, which shifts the solar-geometry curve by the UTC offset (12 h for New Zealand) whenever the input is a local POSIXct.
  • Pressure units. expand_met() returns MET_prsttn / MET_prmslp in Pa. AEME’s station-pressure helper returned hPa despite the documented unit, so an expanded set with no supplied pressure came out a factor of 100 low.
  • Wind components. ds2uv() / uv2ds() are mutual inverses using the meteorological “direction from” convention. AEME::expand_met() swapped u and v when deriving components from speed and direction.

Example

Runs on the bundled Lake Rotorua sample data (three years of hourly ERA5-Land, the matching buoy record, and daily CMIP6 / CCAM projections).

library(metscale)
ex <- system.file("extdata", package = "metscale")

era5 <- read.csv(file.path(ex, "rotorua_era5_hourly_met.csv.gz"), check.names = FALSE)
era5$Date <- as.POSIXct(era5$Date, tz = "Etc/GMT-12")
obs  <- prepare_obs_met(file.path(ex, "rotorua_buoy_met_aeme_hr.csv.gz"),
                        resample = "hour", tz = "Etc/GMT-12",
                        wind_height = 2)   # buoy anemometer at 2 m -> 10 m

## 1. bias-correct ERA5 against the buoy, apply, aggregate to daily
bc        <- fit_met_bias_correction(era5, obs)   # leave-one-year-out CV
bc$skill
corrected <- apply_met_bias_correction(era5, bc)
daily     <- met_to_daily(corrected)

## 2. daily CMIP6 projections at the lake (input to the delta-change baseline)
cmip <- extract_climate_point(file.path(ex, "rotorua_cmip6"),
                            lon = 176.2717, lat = -38.0790)

## 3. disaggregate a daily series back to hourly for a sub-daily model
hourly <- disaggregate_met_to_hourly(daily, donor = corrected, method = "fragments",
                                     lon = 176.2717, lat = -38.0790)

## 4. hand off in CF / CMIP names and units: tas [K], pr [kg m-2 s-1], ...
cf <- met_to_cf(hourly)

The full bias-correction → delta-change → disaggregation workflow is walked through in vignette("scenario-workflow").

Installation

# install.packages("pak")
pak::pak("limnotrack/metscale")