Use quantile delta mapping (QDM) to bias adjust forecast, reanalysis

CHANGELOG.md

@@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Changed
+
+- BREAKING: Rewrite forecast download, climate data parsing, adding QDM bias correction and "noleap" calendar. ([@brews](https://github.com/brews), [PR#17](https://github.com/ClimateImpactLab/poreallas/pull/17))
+
 ## [0.3.0] - 2026-08-20
 
 ### Added

README.md

@@ -20,13 +20,16 @@ Key configurations are set through environment variables or a .env file (see `ex
 
 The current configurations are:
 
-* POREALLAS_TAS_FORECAST_URI: URI to the cleaned ECMWF S51 ensemble air temperature Zarr Store.
-* POREALLAS_ERA5_URI: URI to the Zarr Store of cleaned daily ERA5 dataset used for historical climate and impacts analysis.
+* POREALLAS_TAS_FORECAST_URI: URI to the parsed and bias-adjusted ECMWF S51 ensemble air temperature Zarr Store.
+* POREALLAS_ERA5_URI: URI to the Zarr Store of parsed and bias-adjusted daily ERA5 dataset used for historical climate and impacts analysis.
 * POREALLAS_GAMMA_URI: URI to the Zarr Store of "gamma" parameters used when calculating to calculate a mortality response function.
 * POREALLAS_REGIONS_URI: URI to the Zarr Store of region and grid weights or "segment weights".
-* POREALLAS_REGIONS_POLYGONS_URI: URI to geoparquet file with polygons for each region. Used for mapping.
 * POREALLAS_SOCIOECONOMICS_URI: URI to file with each region's GDP per capita (gdppc).
 * POREALLAS_EFFECTS_URI: Optional URI to write Zarr store of projected mortality effects. Will not write output if unset.
+* POREALLAS_REGIONS_POLYGONS_URI: URI to geoparquet file with polygons for each region. Used for mapping.
+* POREALLAS_PARSED_GMFD_URI: URI to the Zarr Store of parsed GMFD air temperature data.
+* POREALLAS_PARSED_ERA5_URI: URI to the Zarr Store of parsed daily ERA5 air temperature data.
+* POREALLAS_PARSED_FORECAST_URI: URI to the Zarr Store of parsed ECMWF S51 ensemble air temperature.
 
 These are used to run the prototype in `scripts/` for downloads, parsing/cleaning, and projecting.
 
@@ -44,7 +47,6 @@ uv run scripts/05-project_effects.py
 
 from the root of this repository.
 
-Alternatively, running `uv run scripts/05-project_effects_marimo.py` will project with a prototype GMFD bias-adjustment.
 
 ### Data and parsing
 

example.env

@@ -1,5 +1,8 @@
-POREALLAS_TAS_FORECAST_URI = "./data/parsed/s51_tas.zarr"
-POREALLAS_ERA5_URI = "./data/era5_daily_tas_1995_2025_regrid.zarr"
+POREALLAS_PARSED_ERA5_URI = "./data/parsed/era5.zarr"
+POREALLAS_PARSED_GMFD_URI = "./data/parsed/gmfd.zarr"
+POREALLAS_PARSED_FORECAST_URI = "./data/parsed/forecast.zarr"
+POREALLAS_TAS_FORECAST_URI = "./data/parsed/forecast_adj.zarr"
+POREALLAS_ERA5_URI = "./data/parsed/era5_adj.zarr"
 POREALLAS_GAMMA_URI = "./data/parsed/gamma.zarr"
 POREALLAS_REGIONS_URI = "./data/parsed/segment_weights.zarr"
 POREALLAS_REGIONS_POLYGONS_URI = "./data/parsed/impact_region.parquet"

pyproject.toml

@@ -23,6 +23,7 @@ dependencies = [
     "xarray[complete]>=2026.4.0",
     "xclim>=0.61.1",
     "xhistogram>=0.3.2",
+    "xsdba>=0.7.0",
 ]
 
 [build-system]

scripts/00-download_s51_tasmin_tasmax.py → scripts/00-download_forecast.py

@@ -1,29 +1,19 @@
-# Notes in prep for mortality projection based on seasonal ENSO forecasts.
+# Download ECMWF S51 minimum and maximum daily temperatures for a forecast beginning in a TARGET_MONTH from 1981 through 2025.
+# The tasmin and tasmax data are written to separate files.
 
-# See https://www.ecmwf.int/en/forecasts/documentation-and-support/seasonal
-# https://cds.climate.copernicus.eu/datasets/seasonal-original-single-levels?tab=download
-# https://iri.columbia.edu/our-expertise/climate/forecasts/seasonal-climate-forecasts/
+import cdsapi
 
-# Daily ERA5:
-# https://cds.climate.copernicus.eu/datasets/derived-era5-single-levels-daily-statistics
-# # ARCO ERA5:
-# https://github.com/google-research/arco-era5
-# Seasonal forecast daily + subdaily
-# https://cds.climate.copernicus.eu/datasets/seasonal-original-single-levels
 
-import cdsapi
+TARGET_MONTH = 5
+START_YEAR = 1981
+STOP_YEAR = 2026
+OUT_DIRECTORY = "./data/raw/s51_hist_tasmin_tasmax/"
 
-client = cdsapi.Client()
 
-# Trying to download -daily forecast data.
-# https://cds.climate.copernicus.eu/datasets/seasonal-original-single-levels?tab=download
 dataset = "seasonal-original-single-levels"
 request = {
     "originating_centre": "ecmwf",
     "system": "51",
-    "variable": ["maximum_2m_temperature_in_the_last_24_hours"],
-    "year": ["2026"],
-    "month": ["05"],
     "day": ["01"],
     "leadtime_hour": [
         "24",
@@ -245,20 +235,28 @@
     "data_format": "netcdf",
 }
 
-client.retrieve(dataset, request, "./data/raw/s51_tasmax.nc")
-# ds_tasmax = xr.open_dataset("download_s51_tasmax.nc")
+client = cdsapi.Client()
+
+# Stuff month string with a leading "0" if there is only a single character.
+target_month = str(TARGET_MONTH).zfill(2)
+request["month"] = [target_month]
+
+for yr in range(START_YEAR, STOP_YEAR + 1):
+    request["year"] = [str(yr)]
+
+    request["variable"] = ["minimum_2m_temperature_in_the_last_24_hours"]
+    out_path = OUT_DIRECTORY + f"tasmin-{yr}-{target_month}.nc"
+
+    print(f"Beginning download to {out_path}")
+    client.retrieve(dataset, request, out_path)
+    print(f"Downloaded to {out_path}")
+
+    request["variable"] = ["maximum_2m_temperature_in_the_last_24_hours"]
 
-request["variable"] = ["minimum_2m_temperature_in_the_last_24_hours"]
-client.retrieve(dataset, request, "./data/raw/s51_tasmin.nc")
-# ds_tasmin = xr.open_dataset("download_s51_tasmin.nc")
+    out_path = OUT_DIRECTORY + f"tasmax-{yr}-{target_month}.nc"
 
-# ds_tas = xr.merge(
-#     [
-#         xr.open_dataset("download_s51_tasmax.nc"),
-#         xr.open_dataset("download_s51_tasmin.nc"),
-#     ],
-#     compat="no_conflicts",
-# )
+    print(f"Beginning download to {out_path}")
+    client.retrieve(dataset, request, out_path)
+    print(f"Downloaded to {out_path}")
 
-# # Estimate daily tas from daily tasmax and daily tasmin.
-# ds_tas["tas"] = (ds_tas["mx2t24"] + ds_tas["mn2t24"]) / 2
+print("All forecast downloads complete")

scripts/01-parse_era5.py

@@ -0,0 +1,139 @@
+# Parse ERA5 data store to prepare for analysis.
+#
+# Run on notebooks.cilresearch.org with container image pangeo/pangeo-notebook:2026.06.04.
+#
+# This script loads and parses ERA5 data. It is run on the cluster
+# because it loads from a petabyte-scale dataset co-located with this cluster.
+# Data regridding also uses a compiled library which can be difficult to install on
+# some platforms, but is readily available on the cluster.
+#
+# More information on the ERA5 data store hosted on GCP:
+# https://console.cloud.google.com/marketplace/product/bigquery-public-data/arco-era5
+# https://github.com/google-research/arco-era5/
+
+import datetime
+import os
+import uuid
+
+import dask
+from dask_gateway import GatewayCluster  # type: ignore[ty:unresolved-import]
+from dotenv import load_dotenv
+import xarray as xr
+import xesmf as xe  # type: ignore[ty:unresolved-import]
+
+load_dotenv()
+
+OUT_ZARR = os.environ["POREALLAS_PARSED_ERA5_URI"]
+START_YEAR = 1981
+STOP_YEAR = 2025
+TARGET_REGRID_URI = "s51_hcm.nc"
+JUPYTER_IMAGE = os.environ.get("JUPYTER_IMAGE")
+UID = str(uuid.uuid4())
+START_TIME = datetime.datetime.now(datetime.UTC).isoformat()
+
+print(
+    f"""
+        {JUPYTER_IMAGE=}
+        {START_TIME=}
+        {UID=}
+    """
+)
+
+
+def open_regrid_target(uri: str) -> xr.Dataset:
+    """Open/clean a dataset to use as a regridding target"""
+    # Using the S51 seasonal monthly seasonal hindcast ensemble mean from copernicus as the target grid for our regrid...
+    # Selecting so only have coords for latitude and longitude for regridding.
+    target = xr.open_dataset(uri).isel(
+        {"forecast_reference_time": 0, "forecastMonth": 0}, drop=True
+    )
+    return target
+
+
+def open_era5(
+    start_year: int | str,
+    stop_year: int | str,
+    uri="gs://gcp-public-data-arco-era5/ar/full_37-1h-0p25deg-chunk-1.zarr-v3",
+) -> xr.Dataset:
+    """Opens and parses Googe's ARCO ERA5 store, returning rechunked daily tas dataset
+
+    This can be very heavy and data-IO intensive.
+    """
+    ds = xr.open_zarr(
+        uri,
+        chunks=None,
+        storage_options=dict(token="anon"),
+    )
+
+    # Grab only valid periods
+    ar_full_37_1h = ds.sel(
+        time=slice(ds.attrs["valid_time_start"], ds.attrs["valid_time_stop"])
+    )
+
+    # This is huge so only get what we need. It also needs to be chunked so
+    # it isn't read all into memory at once.
+    clipped_window = (
+        ar_full_37_1h["2m_temperature"]
+        .sel(time=slice(str(start_year), str(stop_year)))
+        .chunk({"time": "auto", "latitude": -1, "longitude": -1})
+    )
+
+    # Collect the subdaily data into daily means and rechunk again.
+    daily = clipped_window.resample(time="D").mean()
+    clipped_window_daily = daily.chunk(
+        {"time": "auto", "latitude": -1, "longitude": -1}
+    )
+
+    # We made it a DataArray but let's make it "tas" in a Dataset.
+    clipped_window_daily.name = "tas"
+    out_ds = clipped_window_daily.to_dataset()
+
+    # Add metadata from the full-sized data.
+    out_ds.attrs |= ds.attrs
+    return out_ds
+
+
+dask.config.set({"distributed.comm.timeouts.connect": "60s"})
+cluster = GatewayCluster(worker_image=JUPYTER_IMAGE, scheduler_image=JUPYTER_IMAGE)
+client = cluster.get_client()
+print(client.dashboard_link)
+cluster.scale(50)
+
+regrid_target = open_regrid_target(TARGET_REGRID_URI)
+
+era5 = open_era5(
+    start_year=START_YEAR,
+    stop_year=STOP_YEAR,
+)
+
+# Cannot have leap years in QDM bias adjustment so convert to a no-leapyear calendar.
+era5 = era5.convert_calendar("noleap", dim="time")
+
+regridder = xe.Regridder(era5, regrid_target, method="bilinear", periodic=True)
+era5_regrid = regridder(era5)
+era5_regrid.attrs |= era5.attrs
+
+# Metadata on units is required later in the workflow.
+era5_regrid["tas"].attrs["units"] = "K"
+
+# Add additional general metadata.
+era5_regrid.attrs |= {
+    "poreallas_created_at": START_TIME,
+    "poreallas_uid": UID,
+    "poreallas_description": "Parsed ERA5 climate fields",
+}
+era5_regrid["tas"].attrs |= {
+    "poreallas_created_at": START_TIME,
+    "poreallas_uid": UID,
+    "poreallas_description": "Parsed ERA5 tas field",
+}
+
+# All of time needs to be in a single chunk for QDM bias adjustment.
+# This generally gets you ~110 MiB chunks.
+era5_regrid = era5_regrid.chunk({"time": -1, "latitude": 30, "longitude": 60})
+
+era5_regrid.to_zarr(OUT_ZARR, consolidated=True)
+print(f"Output written to {OUT_ZARR}")
+
+cluster.scale(0)
+cluster.shutdown()