From f7cbae56637e5e3d89bb77aac08c431e954d0dc7 Mon Sep 17 00:00:00 2001 From: Brewster Malevich Date: Mon, 3 Aug 2026 17:47:50 -0700 Subject: [PATCH 1/5] Use quantile delta mapping (QDM) to bias adjust forecast, reanalysis --- README.md | 10 +- example.env | 7 +- pyproject.toml | 1 + ...smin_tasmax.py => 00-download_forecast.py} | 60 +++--- scripts/01-download_and_parse_era5.py | 103 ---------- scripts/01-parse_era5.py | 139 +++++++++++++ scripts/01-parse_forecast.py | 87 ++++++++ scripts/01-parse_gmfd.py | 105 ++++++++++ ...repare_GMFD_ERA5_bias_adjustment_climos.py | 188 ------------------ scripts/02-bias_adjust_era5.py | 69 +++++++ scripts/02-bias_adjust_forecast.py | 163 +++++++++++++++ scripts/02-clean_estimate_s51_tas.py | 92 --------- scripts/05-project_effects.py | 85 +++++++- uv.lock | 140 +++++++++++++ 14 files changed, 828 insertions(+), 421 deletions(-) rename scripts/{00-download_s51_tasmin_tasmax.py => 00-download_forecast.py} (71%) delete mode 100644 scripts/01-download_and_parse_era5.py create mode 100644 scripts/01-parse_era5.py create mode 100644 scripts/01-parse_forecast.py create mode 100644 scripts/01-parse_gmfd.py delete mode 100644 scripts/01-prepare_GMFD_ERA5_bias_adjustment_climos.py create mode 100644 scripts/02-bias_adjust_era5.py create mode 100644 scripts/02-bias_adjust_forecast.py delete mode 100644 scripts/02-clean_estimate_s51_tas.py diff --git a/README.md b/README.md index b31cc5f..e2882a6 100644 --- a/README.md +++ b/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 diff --git a/example.env b/example.env index 753421e..dff7f4e 100644 --- a/example.env +++ b/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" diff --git a/pyproject.toml b/pyproject.toml index 7ce7981..5925c6e 100644 --- a/pyproject.toml +++ b/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] diff --git a/scripts/00-download_s51_tasmin_tasmax.py b/scripts/00-download_forecast.py similarity index 71% rename from scripts/00-download_s51_tasmin_tasmax.py rename to scripts/00-download_forecast.py index 6411693..5ef0d15 100644 --- a/scripts/00-download_s51_tasmin_tasmax.py +++ b/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") diff --git a/scripts/01-download_and_parse_era5.py b/scripts/01-download_and_parse_era5.py deleted file mode 100644 index 656097d..0000000 --- a/scripts/01-download_and_parse_era5.py +++ /dev/null @@ -1,103 +0,0 @@ -# Running this on notebooks.cilresearch.org with pangeo/pangeo-notebook:2026.04.29 -# -# 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. -# -# If you run this on a remote cluster you'll need to ensure the process has access to the -# required data in the ./data directory. Be sure to download the processed data -# to your ./data/parsed/ directory, also. - -import datetime -import os -import uuid - -import dask -from dotenv import load_dotenv -import xarray as xr -import xesmf as xe -from dask_gateway import GatewayCluster - -load_dotenv() - -START_YEAR = 1993 # Need at least 1993-2016 to match Copernicus CDS hindcast period for seasonal forecasts. -STOP_YEAR = 2025 -TARGET_REGRID_URI = "s51_hcm.nc" -OUT_ZARR = os.environ["POREALLAS_ERA5_URI"] - - -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=} - """ -) - -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) - -# Get ERA5 from Google -# https://console.cloud.google.com/marketplace/product/bigquery-public-data/arco-era5 -# https://github.com/google-research/arco-era5/ -ds = xr.open_zarr( - "gs://gcp-public-data-arco-era5/ar/full_37-1h-0p25deg-chunk-1.zarr-v3", - chunks=None, - storage_options=dict(token="anon"), -) -ar_full_37_1h = ds.sel( - time=slice(ds.attrs["valid_time_start"], ds.attrs["valid_time_stop"]) -) -# This is multiple TiB. - -# Want climatology so last 30 years-ish. -# Chunk so doesn't read all data in 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}) -) -# This is ~1 TiB. - -annual_tas = ( - clipped_window.groupby("time.year") - .mean("time") - .to_dataset() - .chunk({"year": "auto", "latitude": -1, "longitude": -1}) - .rename_vars({"2m_temperature": "tas"}) -) - - -clipped_window_daily = clipped_window.resample(time="D").mean() -clipped_window_daily = clipped_window_daily.chunk( - {"time": "auto", "latitude": -1, "longitude": -1} -) - -# 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(TARGET_REGRID_URI).isel( - {"forecast_reference_time": 0, "forecastMonth": 0}, drop=True -) -regridder = xe.Regridder(clipped_window_daily, target, method="bilinear", periodic=True) -clipped_window_daily_regrid = regridder(clipped_window_daily) - -clipped_window_daily_regrid.name = "tas" -clipped_window_daily_regrid.attrs |= clipped_window_daily.attrs - -# Seems to be an xarray bug? This only runs if we first compute() like this: -clipped_window_daily_regrid.to_dataset().chunk( - {"time": "auto", "latitude": -1, "longitude": -1} -).compute().to_zarr(OUT_ZARR, consolidated=False) -print(f"Output written to {OUT_ZARR}") - - -cluster.scale(0) -cluster.shutdown() diff --git a/scripts/01-parse_era5.py b/scripts/01-parse_era5.py new file mode 100644 index 0000000..e45c30f --- /dev/null +++ b/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 +from dotenv import load_dotenv +import xarray as xr +import xesmf as xe + +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() diff --git a/scripts/01-parse_forecast.py b/scripts/01-parse_forecast.py new file mode 100644 index 0000000..d38a49c --- /dev/null +++ b/scripts/01-parse_forecast.py @@ -0,0 +1,87 @@ +# Parse ECMWF S51 forcast/hindcast files to create daily temperature data ready +# for the analysis workflow. + +import datetime +import os +import uuid + +from dotenv import load_dotenv +import xarray as xr + + +load_dotenv() + +OUT_ZARR = os.environ["POREALLAS_PARSED_FORECAST_URI"] +RAW_FORECAST_FILE_PATTERN = "./data/raw/s51_hist_tasmin_tasmax/{var}-{year}-05.nc" +START_YEAR = 1981 +STOP_YEAR = 2026 +UID = str(uuid.uuid4()) +START_TIME = datetime.datetime.now(datetime.UTC).isoformat() + +print( + f""" + {START_TIME=} + {UID=} + """ +) + +# # We are intentionally merging hindcast and forecast datasets even though they +# have different ensemble sizes (the "number" dim). The hindcast gets stuffed +# with NaNs. This is needed when training QDM on the hindcast data and applying +# to forecast data, otherwise the forecast ensemble size gets cut to the hindcast +# size. +target_tasmax_paths = [ + RAW_FORECAST_FILE_PATTERN.format(year=yr, var="tasmax") + for yr in range(START_YEAR, STOP_YEAR + 1) +] +target_tasmin_paths = [ + RAW_FORECAST_FILE_PATTERN.format(year=yr, var="tasmin") + for yr in range(START_YEAR, STOP_YEAR + 1) +] + +mfdataset_kwargs = dict( + join="outer", + coords="different", + compat="no_conflicts", +) +s51 = ( + xr.open_mfdataset(target_tasmax_paths, **mfdataset_kwargs)["mx2t24"] + + xr.open_mfdataset(target_tasmin_paths, **mfdataset_kwargs)["mn2t24"] +) / 2 +s51.name = "tas" +s51 = s51.to_dataset() + + +# Make "valid_time" the "time" dim and main time dim and reducing it to a single +# dimension so xsdba's QDM can run this against the `ref` data. Version of +# xsdba we're running with seems to require a "time" dim to train QDM. So, this +# *needs* to be named "time". +s51 = s51.set_coords("valid_time").rename({"valid_time": "time"}) +# Collapse multidimension time-like dims into to one dimension we can swap in "time" for. +s51 = ( + s51.stack({"_time_placeholder": ["forecast_reference_time", "forecast_period"]}) + .swap_dims({"_time_placeholder": "time"}) + .squeeze(drop=True) + .drop_vars("_time_placeholder") +) + +# Need matching calendars. Removing leap days makes QDM easier. +s51 = s51.convert_calendar("noleap", dim="time") + +# Add additional general metadata. +s51.attrs |= { + "poreallas_created_at": START_TIME, + "poreallas_uid": UID, + "poreallas_description": "Parsed ECMWF S51 ensemble fields", +} +s51["tas"].attrs |= { + "poreallas_created_at": START_TIME, + "poreallas_uid": UID, + "poreallas_description": "Parsed ECMWF S51 ensemble tas fields", +} + +# Rechunking because all of "time", or whatever we're grouping QDM on, needs to be in one chunk. +s51 = s51.chunk({"number": -1, "time": -1, "latitude": 30, "longitude": "auto"}) + +s51.to_zarr(OUT_ZARR, consolidated=True) +print(f"Output written to {OUT_ZARR}") diff --git a/scripts/01-parse_gmfd.py b/scripts/01-parse_gmfd.py new file mode 100644 index 0000000..35f5cf4 --- /dev/null +++ b/scripts/01-parse_gmfd.py @@ -0,0 +1,105 @@ +# Clean archived GMFD NetCDF files. +# +# Run on notebooks.cilresearch.org with container image pangeo/pangeo-notebook:2026.06.04. +# The raw GMFD files this script is parsing are available in our internal cloud storage. +# This is run on the cluster 1) in order to access the raw dataset; 2) Data regridding also +# uses a compiled library which can be difficult to install on some platforms, but is readily +# available on the cluster. + +import datetime +import os +import uuid + +import dask +from dask_gateway import GatewayCluster +from dotenv import load_dotenv +import xarray as xr +import xesmf as xe + +load_dotenv() + + +OUT_ZARR = os.environ["POREALLAS_PARSED_GMFD_URI"] +START_YEAR = 1981 +STOP_YEAR = 2010 +TARGET_REGRID_URI = "s51_hcm.nc" +GMFD_FILE_PATTERN = ( + "/gcs/impactlab-data/climate/source_data/GMFD/tas/tas_0p25_daily_{year}-{year}.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_gmfd(file_pattern: str, start_year: int, stop_year: int) -> xr.Dataset: + """Open the GMFD dataset for a range of years""" + # We have these stored for years 1950 - 2010. Pattern for file name: + target_paths = [ + file_pattern.format(year=yr) for yr in range(start_year, stop_year + 1) + ] + gmfd = xr.open_mfdataset(target_paths) + gmfd.attrs["source_uris"] = str(target_paths) + return gmfd + + +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) + +gmfd = open_gmfd( + file_pattern=GMFD_FILE_PATTERN, + start_year=START_YEAR, + stop_year=STOP_YEAR, +) +# Cannot have leap years in QDM bias adjustment. +gmfd = gmfd.convert_calendar("noleap", dim="time") + +regridder = xe.Regridder(gmfd, regrid_target, method="bilinear", periodic=True) +gmfd_regrid = regridder(gmfd) +gmfd_regrid.attrs |= gmfd.attrs + +# Metadata on units is required later in the workflow. +gmfd_regrid["tas"].attrs["units"] = "K" + +# Add additional general metadata. +gmfd_regrid.attrs |= { + "poreallas_created_at": START_TIME, + "poreallas_uid": UID, + "poreallas_description": "Parsed GMFD climate fields", +} +gmfd_regrid["tas"].attrs |= { + "poreallas_created_at": START_TIME, + "poreallas_uid": UID, + "poreallas_description": "Parsed GMFD tas field", +} + +# All of time needs to be in a single chunk for QDM bias adjustment. +gmfd_regrid = gmfd_regrid.chunk({"time": -1, "latitude": 30, "longitude": "auto"}) + +gmfd_regrid.to_zarr(OUT_ZARR, consolidated=True) +print(f"Output written to {OUT_ZARR}") + +cluster.scale(0) +cluster.shutdown() diff --git a/scripts/01-prepare_GMFD_ERA5_bias_adjustment_climos.py b/scripts/01-prepare_GMFD_ERA5_bias_adjustment_climos.py deleted file mode 100644 index f104524..0000000 --- a/scripts/01-prepare_GMFD_ERA5_bias_adjustment_climos.py +++ /dev/null @@ -1,188 +0,0 @@ -import marimo - -__generated_with = "0.23.8" -app = marimo.App() - - -@app.cell -def _(): - import marimo as mo - - return (mo,) - - -@app.cell(hide_code=True) -def _(mo): - mo.md(r""" - # Calculate `tas` climatologies for ERA5 and GMFD for use in bias adjusting the 7-month forecast for ENSO mortality projections - - This notebook needs to be run on notebooks.cilresearch.org for data access. - - Was run with notebook image pangeo/pangeo-notebook:2026.06.04. - - K. McCusker - - July 2026 - """) - return - - -@app.cell -def _(): - import xarray as xr - import numpy as np - - import matplotlib.pyplot as plt - import cartopy.crs as ccrs - import cartopy.feature as cfeature - import xesmf as xe - - return ccrs, cfeature, np, plt, xe, xr - - -@app.cell -def _(): - gpath = "/gcs/impactlab-data/climate/source_data/GMFD/tas/tas_0p25_daily_{year}-{year}.nc" # 1950-2010 - epathzarr = "gs://poreallas-public-20260605/v20260612/era5_daily_tas_regrid.zarr" # 1993-2025. - epath = "/gcs/impactlab-data/climate/source_data/ERA-5/tas/daily/netcdf/v1.1/tas_daily_{year}-{year}.nc" # 1979-2020 - return epath, epathzarr, gpath - - -@app.cell -def _(epathzarr, xr): - # This is the ERA5 data used in the projection step for the baseline - era5zarr = xr.open_dataset(epathzarr) - - # Use for regridding to 1x1deg. - target_grid = era5zarr.isel(time=0, drop=True) - return (target_grid,) - - -@app.cell -def _(epath, gpath, xr): - # these are the .25deg grids - gmfd_grid = ( - xr.open_dataset(gpath.format(year="1990")) - .isel(time=0) - .rename({"lat": "latitude", "lon": "longitude"}) - ) - era5_grid = xr.open_dataset(epath.format(year="1990")).isel(time=0) - era5_grid - return era5_grid, gmfd_grid - - -@app.cell -def _(target_grid): - target_grid - return - - -@app.cell -def _(era5_grid, gmfd_grid, target_grid, xe): - # set up the regridders to regrid GMFD and ERA5 .25deg to the ERA5 1deg grid - gregridder = xe.Regridder(gmfd_grid, target_grid, method="bilinear", periodic=True) - eregridder = xe.Regridder(era5_grid, target_grid, method="bilinear", periodic=True) - return eregridder, gregridder - - -@app.cell -def _(epath, eregridder, gpath, gregridder, np, xr): - # Regrid ERA5 and GMFD for each year - # This could be distributed instead of looped if want to speed it up - YEARS = np.arange(1980, 2011) - - gmfd_regrid_lst = [] - era5_regrid_lst = [] - for year in YEARS: - print(year) - - gmfd = xr.open_dataset(gpath.format(year=year)).rename( - {"lat": "latitude", "lon": "longitude"} - ) - era5 = xr.open_dataset(epath.format(year=year)) - - gmfd_regrid_lst.append(gregridder(gmfd)) - era5_regrid_lst.append(eregridder(era5)) - return era5_regrid_lst, gmfd_regrid_lst - - -@app.cell -def _(era5_regrid_lst, gmfd_regrid_lst, xr): - # combine each year into one dataset - era5_regrid = xr.concat(era5_regrid_lst, dim="time") - gmfd_regrid = xr.concat(gmfd_regrid_lst, dim="time") - return era5_regrid, gmfd_regrid - - -@app.cell -def _(era5_regrid, gmfd_regrid): - # calculate the monthly mean timeseries and then the climatological monthly mean - era5_regrid_climo = ( - era5_regrid.tas.resample(time="MS").mean().groupby("time.month").mean() - ) - gmfd_regrid_climo = ( - gmfd_regrid.tas.resample(time="MS").mean().groupby("time.month").mean() - ) - return era5_regrid_climo, gmfd_regrid_climo - - -@app.cell -def _(era5_regrid_climo, gmfd_regrid_climo): - # save climos to file - era5_regrid_climo.to_netcdf("tas_era5_1deg_monthly_climo_1980-2010.nc") - gmfd_regrid_climo.to_netcdf("tas_gmfd_1deg_monthly_climo_1980-2010.nc") - return - - -@app.cell -def _(xr): - # check the data: - - era5_climo = xr.open_dataset("tas_era5_1deg_monthly_climo_1980-2010.nc") - gmfd_climo = xr.open_dataset("tas_gmfd_1deg_monthly_climo_1980-2010.nc") - return era5_climo, gmfd_climo - - -@app.cell -def _(era5_climo, gmfd_climo): - # Look at the difference in climos - monthly_diff = gmfd_climo - era5_climo - return (monthly_diff,) - - -@app.cell -def _(ccrs, cfeature, monthly_diff, plt): - fg = ( - monthly_diff["tas"] - .isel(month=[5, 6, 7, 8, 9]) - .plot( - col="month", - transform=ccrs.PlateCarree(), - subplot_kws={ - "projection": ccrs.Robinson() - }, # central_longitude=-95, central_latitude=45)}, - cbar_kwargs={"orientation": "horizontal", "shrink": 0.8, "aspect": 40}, - robust=True, - vmin=-5, - vmax=5, - cmap="RdBu_r", - ) - ) - plt.suptitle("tas (degC), GMFD - ERA5 (1980 - 2010 mean)") - - # lets add a coastline to each axis - # great reason to use FacetGrid.map - fg.map(lambda: plt.gca().coastlines()) - fg.map( - lambda: plt.gca().add_feature(cfeature.OCEAN, facecolor="lightgray", zorder=1) - ) - - plt.savefig( - "tas_enso_monthly_tas_difference_gmfd_v_era5_1deg.png", - dpi=300, - ) - return - - -if __name__ == "__main__": - app.run() diff --git a/scripts/02-bias_adjust_era5.py b/scripts/02-bias_adjust_era5.py new file mode 100644 index 0000000..6ef6acf --- /dev/null +++ b/scripts/02-bias_adjust_era5.py @@ -0,0 +1,69 @@ +# Created QDM bias adjustment of the parsed ERA5 dataset using the parsed GMFD dataset. + +import datetime +import os +import uuid + +from dotenv import load_dotenv +import xarray as xr +from xsdba.adjustment import QuantileDeltaMapping + +load_dotenv() + +ERA5_URI = os.environ["POREALLAS_PARSED_ERA5_URI"] +GMFD_URI = os.environ["POREALLAS_PARSED_GMFD_URI"] +OUT_ZARR = os.environ["POREALLAS_ERA5_URI"] +HISTREF_START_YEAR = 1981 +HISTREF_STOP_YEAR = 1997 +SIM_START_YEAR = 1996 +SIM_STOP_YEAR = 2025 +QDM_N_QUANTILES = 10 +UID = str(uuid.uuid4()) +START_TIME = datetime.datetime.now(datetime.UTC).isoformat() + +gmfd = xr.open_zarr(GMFD_URI) +era5 = xr.open_zarr(ERA5_URI) + +ref = gmfd.sel(time=slice(str(HISTREF_START_YEAR), str(HISTREF_STOP_YEAR))) +hist = era5.sel(time=slice(str(HISTREF_START_YEAR), str(HISTREF_STOP_YEAR))) +sim = era5.sel(time=slice(str(SIM_START_YEAR), str(SIM_STOP_YEAR))) + +# # "time" dim cannot be chunked for QDM. +ref = ref.chunk({"time": -1}) +hist = hist.chunk({"time": -1}) +sim = sim.chunk({"time": -1}) + +qdm = QuantileDeltaMapping.train( + ref["tas"], hist["tas"], nquantiles=QDM_N_QUANTILES, kind="+", group="time.month" +) + +sim_adj = qdm.adjust(sim["tas"]) + +sim_adj.name = "tas" +sim_adj = sim_adj.to_dataset() + +# Add additional general metadata. +sim_adj.attrs |= { + "poreallas_created_at": START_TIME, + "poreallas_uid": UID, + "poreallas_description": "QDM bias-adjusted ERA5 climate fields", +} +sim_adj["tas"].attrs |= { + "poreallas_created_at": START_TIME, + "poreallas_uid": UID, + "poreallas_description": "QDM bias-adjusted ERA5 tas field", + "poreallas_adjustment_method": "QDM", + "poreallas_histref_start_year": HISTREF_START_YEAR, + "poreallas_histref_stop_year": HISTREF_STOP_YEAR, + "poreallas_sim_start_year": SIM_START_YEAR, + "poreallas_sim_stop_year": SIM_STOP_YEAR, + "poreallas_qdm_nquantiles": QDM_N_QUANTILES, + "poreallas_ref_uri": GMFD_URI, + "poreallas_hist_uri": ERA5_URI, + "poreallas_sim_uri": ERA5_URI, +} + +sim_adj = sim_adj.chunk("auto") + +sim_adj.to_zarr(OUT_ZARR, consolidated=True) +print(f"Output written to {OUT_ZARR}") diff --git a/scripts/02-bias_adjust_forecast.py b/scripts/02-bias_adjust_forecast.py new file mode 100644 index 0000000..f724972 --- /dev/null +++ b/scripts/02-bias_adjust_forecast.py @@ -0,0 +1,163 @@ +# Create QDM bias adjustment of the parsed forecast dataset using the parsed GMFD dataset. + +import datetime +import os +import uuid + +from dotenv import load_dotenv +import xarray as xr +from xsdba.adjustment import QuantileDeltaMapping, TrainAdjust +from xsdba.base import Grouper + +load_dotenv() + +FORECAST_URI = os.environ["POREALLAS_PARSED_FORECAST_URI"] +GMFD_URI = os.environ["POREALLAS_PARSED_GMFD_URI"] +OUT_ZARR = os.environ["POREALLAS_TAS_FORECAST_URI"] +HISTREF_START_YEAR = 1981 +HISTREF_STOP_YEAR = 1997 +SIM_START_YEAR = 2007 +SIM_STOP_YEAR = 2026 +QDM_N_QUANTILES = 10 +FORECAST_LENGTH = 215 # ECMWF S51 is 215 days. +UID = str(uuid.uuid4()) +START_TIME = datetime.datetime.now(datetime.UTC).isoformat() + + +def adjust_month( + *, + ref: xr.DataArray, + hist: xr.DataArray, + sim: xr.DataArray, + target_month: int, + nquantiles: int, +) -> tuple[TrainAdjust, xr.DataArray]: + """ + Train and apply QDM for a particular `time.month` + """ + ref = ref.where(ref["time.month"] == target_month, drop=True) + hist = hist.where(hist["time.month"] == target_month, drop=True) + sim = sim.where(sim["time.month"] == target_month, drop=True) + + qdm = QuantileDeltaMapping.train( + ref, + hist, + nquantiles=nquantiles, + kind="+", + group=Grouper("time", add_dims=["number"]), + ) + adj = qdm.adjust(sim) + return qdm, adj + + +def adjust_months( + *, + ref: xr.DataArray, + hist: xr.DataArray, + sim: xr.DataArray, + nquantiles: int, +) -> xr.DataArray: + """ + Train and apply quantile delta mapping (QDM) for all `time.month` in a simulation. + + We need a custom algorithm for this because our forecast ensembles run for < + 365 days yet this QDM implementation does not allow us to group by "time.month" + when it does not have all 12 months. So we train and apply QDM to each of the + months in the simulation dataset and then concatenate them back together along + the time dimension. The concatenated data is then sorted by the time dimension + to return the data to chronological order. + + Parameters + ---------- + ref : + Reference dataset to compare against a historical simulation to train a QDM. + hist : + Historical simulation dataset to be compared against ref when training the QDM. + sim : + Simulation, or forecast ensemble to be adjusted by the trained QDM. + nquantiles : + Number of quantiles to use in the quantile mapping. + + Returns + ------- + combined : + Simulated, bias-adjusted by a QDM trained on a historical and reference dataset. + """ + adjusted = [] + for m in set(sim["time.month"].data): + _, adj = adjust_month( + ref=ref, + hist=hist, + sim=sim, + target_month=m, + nquantiles=nquantiles, + ) + adjusted.append(adj) + + combined = xr.concat(adjusted, dim="time").sortby("time") + return combined + + +gmfd = xr.open_zarr(GMFD_URI) +forecast = xr.open_zarr(FORECAST_URI) + +# Outline the datasets we need for the adjustment, grabbing the windows in time needed. +ref = gmfd.sel(time=slice(str(HISTREF_START_YEAR), str(HISTREF_STOP_YEAR))) +hist = forecast.sel(time=slice(str(HISTREF_START_YEAR), str(HISTREF_STOP_YEAR))) +sim = forecast.sel(time=slice(str(SIM_START_YEAR), str(SIM_STOP_YEAR))) + +# # Subset reference to only daysofyear that are in our forecast ensemble. The +# forecast ensemble has incomplete years. Ref/hist/sim need to have matching +# ragged ends in their time series for QDM. +ref = ref.where(ref["time.dayofyear"].isin(sim["time.dayofyear"]), drop=True) + +# Rechunking because all of "time", or whatever we're grouping QDM on, needs to be in one chunk. +ref = ref.chunk({"time": -1}) +hist = hist.chunk({"number": -1, "time": -1, "latitude": "30", "longitude": "auto"}) +sim = sim.chunk({"number": -1, "time": -1, "latitude": "30", "longitude": "auto"}) + +# Train QDM and adjust the forecast ensemble, for the months in the forecast ensemble. +sim_adj = adjust_months( + ref=ref["tas"], + hist=hist["tas"], + sim=sim["tas"], + nquantiles=QDM_N_QUANTILES, +) + +# QDM requires forecasts over several years to estimate distribution functions +# but we're only interested in the most recent forecast (of `FORECAST_LENGTH` +# days; e.g. 215 days for ECMWF S51). So we pop off and use those last +# FORECAST_LENGTH days because that's the real period of interest for the +# mortality projection. +# NOTE: If you turn this off for debugging, the output data can be orders of +# magnitude larger! +sim_adj = sim_adj.isel(time=slice(-int(FORECAST_LENGTH), None)) + +sim_adj.name = "tas" +sim_adj = sim_adj.to_dataset() + +# Add additional general metadata. +sim_adj.attrs |= { + "poreallas_created_at": START_TIME, + "poreallas_uid": UID, + "poreallas_description": "QDM bias-adjusted forecast ensemble fields", +} +sim_adj["tas"].attrs |= { + "poreallas_created_at": START_TIME, + "poreallas_uid": UID, + "poreallas_description": "QDM bias-adjusted forecast ensemble tas fields", + "poreallas_adjustment_method": "QDM", + "poreallas_histref_start_year": HISTREF_START_YEAR, + "poreallas_histref_stop_year": HISTREF_STOP_YEAR, + "poreallas_sim_start_year": SIM_START_YEAR, + "poreallas_sim_stop_year": SIM_STOP_YEAR, + "poreallas_qdm_nquantiles": QDM_N_QUANTILES, + "poreallas_ref_uri": GMFD_URI, + "poreallas_hist_uri": FORECAST_URI, + "poreallas_sim_uri": FORECAST_URI, +} + +sim_adj = sim_adj.chunk("auto") + +sim_adj.to_zarr(OUT_ZARR, consolidated=True) +print(f"Output written to {OUT_ZARR}") diff --git a/scripts/02-clean_estimate_s51_tas.py b/scripts/02-clean_estimate_s51_tas.py deleted file mode 100644 index dd876dc..0000000 --- a/scripts/02-clean_estimate_s51_tas.py +++ /dev/null @@ -1,92 +0,0 @@ -""" -Estimate forecast ensemble daily-average air temperature from daily 24-hour -maximum and minimum air temperature and bias-adjust by comparing monthly -hindcast and reanalysis. - -Write output to POREALLAS_TAS_FORECAST_URI -""" - -import datetime -import os -import uuid - -from dotenv import load_dotenv -import xarray as xr - - -load_dotenv() - -OUT_ZARR = os.environ["POREALLAS_TAS_FORECAST_URI"] -ERA5_URI = os.environ["POREALLAS_ERA5_URI"] -S51_TASMAX_URI = "./data/raw/s51_tasmax.nc" -S51_TASMIN_URI = "./data/raw/s51_tasmin.nc" -HC_URI = "./data/raw/s51_hcm.nc" -# It would be great if we could grab this from CDS/ECMWF file metadata but we can't, so be careful. -HINDCAST_PERIOD = slice("1993", "2016") -UID = str(uuid.uuid4()) -DATETIME_NOW = datetime.datetime.now(datetime.timezone.utc).isoformat() - - -reanalysis = xr.open_dataset(ERA5_URI) - -hc = xr.open_dataset(HC_URI) - - -# Estimate forecast daily tas from daily tasmax and tasmin average. -ds_tas = xr.merge( - [ - xr.open_dataset(S51_TASMAX_URI), - xr.open_dataset(S51_TASMIN_URI), - ], - compat="no_conflicts", -) -# Estimate daily tas from daily tasmax and daily tasmin. -ds_tas["tas"] = (ds_tas["mx2t24"] + ds_tas["mn2t24"]) / 2 - - -# Add a new month variable to represent the calendar month for each valid period of the forecast. -# Squeeze it so it only has the needed "ForecastMonth" dim so we can swap them out in place. -# Need this "month" so can align with reanalysis to estimate forecast bias. -hc["month"] = ( - hc["forecast_reference_time"].dt.month + hc["forecastMonth"] -).squeeze() - 1 -hc = hc.swap_dims({"forecastMonth": "month"}) - - -# Estimate forecast monthly bias: difference between monthly hindcasts and -# average of reanalysis for all 12 months in forecast's handcast period. -months_climatology = reanalysis.sel(time=HINDCAST_PERIOD).groupby("time.month").mean() -bias = (hc["t2m"] - months_climatology["tas"]).squeeze() - - -# Make "valid_time" the "time" dim and main time dim rather than -# "forecast_period", dropping everything else we don't need. -ds_tas = ds_tas.set_coords("valid_time").rename({"valid_time": "time"}) -ds_tas = ds_tas.swap_dims({"forecast_period": "time"}).squeeze(drop=True) - -# Only use months from the forecast-ensemble that are bias-adjustable. If a month is not here, it's likely incomplete. -ds_tas = ds_tas.where(ds_tas["time.month"].isin(bias["month"]), drop=True) -# Apply monthly bias adjustment to forecast ensemble. -adjusted = ds_tas["tas"] - bias.sel(month=ds_tas["time.month"]) - -# Each month also uniquely corresponds to a forecast lead time, so the bias -# adjustment accounts for the combined bias from the forecast ensemble systematic bias in simulating -# 1) a particular month of a year -# 2) lead time from when the forecast was initialized -- or, how long the model has been running - -forecast_ensemble = xr.Dataset({"tas": adjusted}) - -forecast_ensemble.attrs |= { - "poreallas_created_at": DATETIME_NOW, - "poreallas_uid": UID, - "poreallas_description": "Bias-adjusted forecast ensemble daily air temperature", - "poreallas_reanalysis_uri": ERA5_URI, - "poreallas_forecast_tasmax_uri": S51_TASMAX_URI, - "poreallas_forecast_tasmin_uri": S51_TASMIN_URI, - "poreallas_hindcast_uri": HC_URI, -} - -forecast_ensemble = forecast_ensemble.chunk("auto") - -forecast_ensemble.to_zarr(OUT_ZARR, consolidated=False) -print(f"s51 data written to {OUT_ZARR}") diff --git a/scripts/05-project_effects.py b/scripts/05-project_effects.py index 067bcf5..f46bcc9 100644 --- a/scripts/05-project_effects.py +++ b/scripts/05-project_effects.py @@ -164,6 +164,12 @@ def main(): ) fixed_beta["beta_hotonly"].attrs["long_name"] = "Hot temperature mortality rate" + # Do beta, allowing only COLD deaths by 0-ing out everything on the warm side of the minimum-mortality temperature. + fixed_beta["beta_coldonly"] = fixed_beta["beta"].where( + fixed_beta["tas_bin"] < fixed_beta["mmt"], other=0 + ) + fixed_beta["beta_coldonly"].attrs["long_name"] = "Cold temperature mortality rate" + # Project mortality. # Start with forecast ensemble. forecast_input = ( @@ -221,6 +227,34 @@ def main(): "long_name": "Hot temperature mortality", } + # Now cold-only projection + forecast_input = ( + xr.Dataset( + { + "histogram_tas": histogram_forecast_tas["histogram_tas"], + "beta": fixed_beta["beta_coldonly"], + } + ) + .dropna(dim="region") + .chunk( + { + "region": "auto", # "auto" is a sensible default. + "time": -1, + "tas_bin": -1, + "age_cohort": 1, + "number": 1, + }, + ) + .unify_chunks() + ) + projected_forecast_coldonly = isku.project( + forecast_input, model=mortality_effect_model + ).compute() + projected_forecast_coldonly["effect"].attrs = { + "units": "deaths per 100,000 people", + "long_name": "Cold temperature mortality", + } + # Now do the baseline period. # Stick everything together and make sure it aligns and matches. Rechunk all together. Also drop any regions with NaNs. hist_input = ( @@ -273,12 +307,41 @@ def main(): "long_name": "Hot temperature mortality", } + # Now cold-only projection + hist_input = ( + xr.Dataset( + { + "histogram_tas": histogram_hist_tas["histogram_tas"], + "beta": fixed_beta["beta_coldonly"], + } + ) + .dropna(dim="region") + .chunk( + { + "region": "auto", # "auto" is a sensible default. + "time": -1, + "tas_bin": -1, + "age_cohort": 1, + }, + ) + .unify_chunks() + ) + projected_hist_coldonly = isku.project( + hist_input, model=mortality_effect_model + ).compute() + projected_hist_coldonly["effect"].attrs = { + "units": "deaths per 100,000 people", + "long_name": "Cold temperature mortality", + } + # Collect everything and write to storage. _out = { "forecast": projected_forecast, "baseline": projected_hist, "forecast_hotonly": projected_forecast_hotonly, "baseline_hotonly": projected_hist_hotonly, + "forecast_coldonly": projected_forecast_coldonly, + "baseline_coldonly": projected_hist_coldonly, } _out_dt = xr.DataTree.from_dict(_out) @@ -332,8 +395,28 @@ def main(): "poreallas_regions_uri": REGIONS_URI, } + _out_dt["forecast_coldonly"].attrs |= { + "poreallas_created_at": _datetime_now, + "poreallas_uid": _uid, + "poreallas_description": "Forecast ensemble projected cold temperature mortality effects", + "poreallas_temperature_uri": TAS_FORECAST_URI, + "poreallas_socioeconomics_uri": SOCIOECONOMICS_URI, + "poreallas_model_parameters_uri": GAMMA_URI, + "poreallas_regions_uri": REGIONS_URI, + } + + _out_dt["baseline_coldonly"].attrs |= { + "poreallas_created_at": _datetime_now, + "poreallas_uid": _uid, + "poreallas_description": "Baseline projected cold temperature mortality effects", + "poreallas_temperature_uri": ERA5_URI, + "poreallas_socioeconomics_uri": SOCIOECONOMICS_URI, + "poreallas_model_parameters_uri": GAMMA_URI, + "poreallas_regions_uri": REGIONS_URI, + } + if EFFECTS_URI is not None: - _out_dt.to_zarr(EFFECTS_URI, consolidated=False) + _out_dt.to_zarr(EFFECTS_URI, consolidated=True) print(f"Effects written to {EFFECTS_URI}") diff --git a/uv.lock b/uv.lock index 0991161..1b94cc3 100644 --- a/uv.lock +++ b/uv.lock @@ -87,6 +87,15 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/fb/76/641ae371508676492379f16e2fa48f4e2c11741bd63c48be4b12a6b09cba/aiosignal-1.4.0-py3-none-any.whl", hash = "sha256:053243f8b92b990551949e63930a839ff0cf0b0ebbe0597b0f3fb19e1a0fe82e", size = 7490, upload-time = "2025-07-03T22:54:42.156Z" }, ] +[[package]] +name = "annotated-doc" +version = "0.0.5" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/5a/8e/38aa427ed5402449e226975b649c5dc73ccadfefeb95e6aecb8f8ea4b6b6/annotated_doc-0.0.5.tar.gz", hash = "sha256:c7e58ce09192557605d8bbd92836d7e1d520ac9580096042c0bfd197efacf1bb", size = 10758, upload-time = "2026-07-28T13:50:58.129Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/3e/30/e900b21425a860e195f32e37657aa1f7c7f2b1bfb26f03ca209b90933c06/annotated_doc-0.0.5-py3-none-any.whl", hash = "sha256:117bac03a25ede5df5440e855b32d556049ca169ead221505badf432fed4b101", size = 5302, upload-time = "2026-07-28T13:50:57.239Z" }, +] + [[package]] name = "anyio" version = "4.13.0" @@ -1261,6 +1270,15 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/7b/91/984aca2ec129e2757d1e4e3c81c3fcda9d0f85b74670a094cc443d9ee949/joblib-1.5.3-py3-none-any.whl", hash = "sha256:5fc3c5039fc5ca8c0276333a188bbd59d6b7ab37fe6632daa76bc7f9ec18e713", size = 309071, upload-time = "2025-12-15T08:41:44.973Z" }, ] +[[package]] +name = "jsonpickle" +version = "4.1.2" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/8d/c0/dde9b4b42cc415b9579573f967f12efbb034e427a2a37e93ad5139891d87/jsonpickle-4.1.2.tar.gz", hash = "sha256:8afed18aa189fd81e2e833b426bb4af485594921f0b1d36c2001fc5637a2f210", size = 319120, upload-time = "2026-05-28T03:50:11.892Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/a1/7b/fd3c7a09649aea9da1d3587aea624d8f9b29963dfd84a1bdb2aa93b36dac/jsonpickle-4.1.2-py3-none-any.whl", hash = "sha256:7ffe34426bc797684dbf1dc84185558bd864cd25b1ff5fb01b7405e392d0a937", size = 47203, upload-time = "2026-05-28T03:50:10.605Z" }, +] + [[package]] name = "kiwisolver" version = "1.5.0" @@ -1500,6 +1518,18 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/de/1f/77fa3081e4f66ca3576c896ae5d31c3002ac6607f9747d2e3aa49227e464/markdown-3.10.2-py3-none-any.whl", hash = "sha256:e91464b71ae3ee7afd3017d9f358ef0baf158fd9a298db92f1d4761133824c36", size = 108180, upload-time = "2026-02-09T14:57:25.787Z" }, ] +[[package]] +name = "markdown-it-py" +version = "4.2.0" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "mdurl" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/06/ff/7841249c247aa650a76b9ee4bbaeae59370dc8bfd2f6c01f3630c35eb134/markdown_it_py-4.2.0.tar.gz", hash = "sha256:04a21681d6fbb623de53f6f364d352309d4094dd4194040a10fd51833e418d49", size = 82454, upload-time = "2026-05-07T12:08:28.36Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/b3/81/4da04ced5a082363ecfa159c010d200ecbd959ae410c10c0264a38cac0f5/markdown_it_py-4.2.0-py3-none-any.whl", hash = "sha256:9f7ebbcd14fe59494226453aed97c1070d83f8d24b6fc3a3bcf9a38092641c4a", size = 91687, upload-time = "2026-05-07T12:08:27.182Z" }, +] + [[package]] name = "markupsafe" version = "3.0.3" @@ -1563,6 +1593,15 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/78/23/92493c3e6e1b635ccfff146f7b99e674808787915420373ac399283764c2/matplotlib-3.10.9-cp314-cp314t-win_arm64.whl", hash = "sha256:a49f1eadc84ca85fd72fa4e89e70e61bf86452df6f971af04b12c60761a0772c", size = 8324785, upload-time = "2026-04-24T00:13:53.633Z" }, ] +[[package]] +name = "mdurl" +version = "0.1.2" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/d6/54/cfe61301667036ec958cb99bd3efefba235e65cdeb9c84d24a8293ba1d90/mdurl-0.1.2.tar.gz", hash = "sha256:bb413d29f5eea38f31dd4754dd7377d4465116fb207585f97bf925588687c1ba", size = 8729, upload-time = "2022-08-14T12:40:10.846Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/b3/38/89ba8ad64ae25be8de66a6d463314cf1eb366222074cfda9ee839c56a4b4/mdurl-0.1.2-py3-none-any.whl", hash = "sha256:84008a41e51615a49fc9966191ff91509e3c40b939176e643fd50a5c2196b8f8", size = 9979, upload-time = "2022-08-14T12:40:09.779Z" }, +] + [[package]] name = "metacsv" version = "0.1.1" @@ -1915,6 +1954,18 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/71/e7/40fb618334dcdf7c5a316c0e7343c5cd82d3d866edc100d98e29bc945ecd/partd-1.4.2-py3-none-any.whl", hash = "sha256:978e4ac767ec4ba5b86c6eaa52e5a2a3bc748a2ca839e8cc798f1cc6ce6efb0f", size = 18905, upload-time = "2024-05-06T19:51:39.271Z" }, ] +[[package]] +name = "patsy" +version = "1.0.2" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "numpy" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/be/44/ed13eccdd0519eff265f44b670d46fbb0ec813e2274932dc1c0e48520f7d/patsy-1.0.2.tar.gz", hash = "sha256:cdc995455f6233e90e22de72c37fcadb344e7586fb83f06696f54d92f8ce74c0", size = 399942, upload-time = "2025-10-20T16:17:37.535Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/f1/70/ba4b949bdc0490ab78d545459acd7702b211dfccf7eb89bbc1060f52818d/patsy-1.0.2-py2.py3-none-any.whl", hash = "sha256:37bfddbc58fcf0362febb5f54f10743f8b21dd2aa73dec7e7ef59d1b02ae668a", size = 233301, upload-time = "2025-10-20T16:17:36.563Z" }, +] + [[package]] name = "pillow" version = "12.2.0" @@ -2033,6 +2084,7 @@ dependencies = [ { name = "xarray", extra = ["complete"] }, { name = "xclim" }, { name = "xhistogram" }, + { name = "xsdba" }, ] [package.dev-dependencies] @@ -2060,6 +2112,7 @@ requires-dist = [ { name = "xarray", extras = ["complete"], specifier = ">=2026.4.0" }, { name = "xclim", specifier = ">=0.61.1" }, { name = "xhistogram", specifier = ">=0.3.2" }, + { name = "xsdba", specifier = ">=0.7.0" }, ] [package.metadata.requires-dev] @@ -2533,6 +2586,19 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/3b/5d/63d4ae3b9daea098d5d6f5da83984853c1bbacd5dc826764b249fe119d24/requests_oauthlib-2.0.0-py2.py3-none-any.whl", hash = "sha256:7dd8a5c40426b779b0868c404bdef9768deccf22749cde15852df527e6269b36", size = 24179, upload-time = "2024-03-22T20:32:28.055Z" }, ] +[[package]] +name = "rich" +version = "15.0.0" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "markdown-it-py" }, + { name = "pygments" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/c0/8f/0722ca900cc807c13a6a0c696dacf35430f72e0ec571c4275d2371fca3e9/rich-15.0.0.tar.gz", hash = "sha256:edd07a4824c6b40189fb7ac9bc4c52536e9780fbbfbddf6f1e2502c31b068c36", size = 230680, upload-time = "2026-04-12T08:24:00.75Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/82/3b/64d4899d73f91ba49a8c18a8ff3f0ea8f1c1d75481760df8c68ef5235bf5/rich-15.0.0-py3-none-any.whl", hash = "sha256:33bd4ef74232fb73fe9279a257718407f169c09b78a87ad3d296f548e27de0bb", size = 310654, upload-time = "2026-04-12T08:24:02.83Z" }, +] + [[package]] name = "ruff" version = "0.15.15" @@ -2657,6 +2723,15 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/9a/f6/f09272a71976dfc138129b8faf435d064a811ae2f708cb147dccdf7aacdb/shapely-2.1.2-cp314-cp314t-win_amd64.whl", hash = "sha256:0036ac886e0923417932c2e6369b6c52e38e0ff5d9120b90eef5cd9a5fc5cae9", size = 1796682, upload-time = "2025-09-24T13:51:39.233Z" }, ] +[[package]] +name = "shellingham" +version = "1.5.4" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/58/15/8b3609fd3830ef7b27b655beb4b4e9c62313a4e8da8c676e142cc210d58e/shellingham-1.5.4.tar.gz", hash = "sha256:8dbca0739d487e5bd35ab3ca4b36e11c4078f3a234bfce294b0a0291363404de", size = 10310, upload-time = "2023-10-24T04:13:40.426Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/e0/f9/0595336914c5619e5f28a1fb793285925a8cd4b432c9da0a987836c7f822/shellingham-1.5.4-py2.py3-none-any.whl", hash = "sha256:7ecfff8f2fd72616f7481040475a65b2bf8af90a56c89140852d1120324e8686", size = 9755, upload-time = "2023-10-24T04:13:38.866Z" }, +] + [[package]] name = "six" version = "1.17.0" @@ -2736,6 +2811,27 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/1c/54/196d0c1db10af76baa4f64894448505d60d3cdf70ef92cbb35f46a4e4c71/starlette-1.2.1-py3-none-any.whl", hash = "sha256:4de0082d08c8f6764a85a54cf1120d6939507a19905c7768acad2a9f875d2b89", size = 73350, upload-time = "2026-05-31T01:07:50.09Z" }, ] +[[package]] +name = "statsmodels" +version = "0.14.6" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "numpy" }, + { name = "packaging" }, + { name = "pandas" }, + { name = "patsy" }, + { name = "scipy" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/0d/81/e8d74b34f85285f7335d30c5e3c2d7c0346997af9f3debf9a0a9a63de184/statsmodels-0.14.6.tar.gz", hash = "sha256:4d17873d3e607d398b85126cd4ed7aad89e4e9d89fc744cdab1af3189a996c2a", size = 20689085, upload-time = "2025-12-05T23:08:39.522Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/71/de/09540e870318e0c7b58316561d417be45eff731263b4234fdd2eee3511a8/statsmodels-0.14.6-cp314-cp314-macosx_10_15_x86_64.whl", hash = "sha256:00781869991f8f02ad3610da6627fd26ebe262210287beb59761982a8fa88cae", size = 10069403, upload-time = "2025-12-05T23:12:48.424Z" }, + { url = "https://files.pythonhosted.org/packages/ab/f0/63c1bfda75dc53cee858006e1f46bd6d6f883853bea1b97949d0087766ca/statsmodels-0.14.6-cp314-cp314-macosx_11_0_arm64.whl", hash = "sha256:73f305fbf31607b35ce919fae636ab8b80d175328ed38fdc6f354e813b86ee37", size = 9989253, upload-time = "2025-12-05T23:13:05.274Z" }, + { url = "https://files.pythonhosted.org/packages/c1/98/b0dfb4f542b2033a3341aa5f1bdd97024230a4ad3670c5b0839d54e3dcab/statsmodels-0.14.6-cp314-cp314-manylinux2014_aarch64.manylinux_2_17_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:e443e7077a6e2d3faeea72f5a92c9f12c63722686eb80bb40a0f04e4a7e267ad", size = 10090802, upload-time = "2025-12-05T23:13:20.653Z" }, + { url = "https://files.pythonhosted.org/packages/34/0e/2408735aca9e764643196212f9069912100151414dd617d39ffc72d77eee/statsmodels-0.14.6-cp314-cp314-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:3414e40c073d725007a6603a18247ab7af3467e1af4a5e5a24e4c27bc26673b4", size = 10337587, upload-time = "2025-12-05T23:13:37.597Z" }, + { url = "https://files.pythonhosted.org/packages/0f/36/4d44f7035ab3c0b2b6a4c4ebb98dedf36246ccbc1b3e2f51ebcd7ac83abb/statsmodels-0.14.6-cp314-cp314-musllinux_1_2_x86_64.whl", hash = "sha256:a518d3f9889ef920116f9fa56d0338069e110f823926356946dae83bc9e33e19", size = 10363350, upload-time = "2025-12-05T23:13:53.08Z" }, + { url = "https://files.pythonhosted.org/packages/26/33/f1652d0c59fa51de18492ee2345b65372550501ad061daa38f950be390b6/statsmodels-0.14.6-cp314-cp314-win_amd64.whl", hash = "sha256:151b73e29f01fe619dbce7f66d61a356e9d1fe5e906529b78807df9189c37721", size = 9588010, upload-time = "2025-12-05T23:14:07.28Z" }, +] + [[package]] name = "tblib" version = "3.2.2" @@ -2826,6 +2922,21 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/04/87/369056ed46f1b235130ec0595393262f9cd2061ca3dab276d490980f9343/ty-0.0.40-py3-none-win_arm64.whl", hash = "sha256:07da2b09d9130e2c9a257d2a29beb53105835b0256ee5fdb288fe1aab83fee47", size = 11117369, upload-time = "2026-05-27T17:55:39.329Z" }, ] +[[package]] +name = "typer" +version = "0.27.0" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "annotated-doc" }, + { name = "colorama", marker = "sys_platform == 'win32'" }, + { name = "rich" }, + { name = "shellingham" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/37/78/fda3361b56efc27944f24225f6ecd13d96d6fcfe37bd0eb34e2f4c63f9fc/typer-0.27.0.tar.gz", hash = "sha256:629bd12ea5d13a17148125d9a264f949eb171fb3f120f9b04d85873cab054fa5", size = 203430, upload-time = "2026-07-15T19:21:07.007Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/40/03/26a383c9e58c213199d1aad1c3d353cfc22d4444ec6d2c0bf8ad02523843/typer-0.27.0-py3-none-any.whl", hash = "sha256:6f4b27631e47f077871b7dc30e933ec0131c1390fbe0e387ea5574b5bac9ccf1", size = 122716, upload-time = "2026-07-15T19:21:05.553Z" }, +] + [[package]] name = "typing-extensions" version = "4.15.0" @@ -2997,6 +3108,35 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/18/08/1432dd10193a5d45294bd42042a5631259ee5a12cd2e9075350546d07a03/xhistogram-0.3.2-py3-none-any.whl", hash = "sha256:ad55330d55296d273b3370678223fde0f50085e04cb744c7b3b0bb7702a2c6bf", size = 20227, upload-time = "2022-09-20T05:08:01.428Z" }, ] +[[package]] +name = "xsdba" +version = "0.7.0" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "boltons" }, + { name = "bottleneck" }, + { name = "cf-xarray" }, + { name = "cftime" }, + { name = "dask", extra = ["array"] }, + { name = "filelock" }, + { name = "jsonpickle" }, + { name = "matplotlib" }, + { name = "numba" }, + { name = "numpy" }, + { name = "packaging" }, + { name = "pandas" }, + { name = "pint" }, + { name = "rich" }, + { name = "scipy" }, + { name = "statsmodels" }, + { name = "typer" }, + { name = "xarray" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/3b/c0/e5ef3d69c24b19d958d48730326aaf651c5070844e34a8860083eb6dd4fb/xsdba-0.7.0.tar.gz", hash = "sha256:9fa47be804b4a478c761517504a99499c1d0fedfb231846d7843b24f04e28b5f", size = 240270, upload-time = "2026-06-16T17:32:18.663Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/3d/10/a1f1760dfc5e80bb68e0156c3a57ce1f07787ed660ed977cf5c0c48e3a14/xsdba-0.7.0-py3-none-any.whl", hash = "sha256:89a9aad5db4240346b412ee725a3ba96f0c856f7e8999cabfac9d32d8d76671f", size = 136932, upload-time = "2026-06-16T17:32:17.264Z" }, +] + [[package]] name = "xyzservices" version = "2026.3.0" From c7a32bd055d883d5869b784792f3614217f73c43 Mon Sep 17 00:00:00 2001 From: Brewster Malevich Date: Mon, 3 Aug 2026 20:21:56 -0700 Subject: [PATCH 2/5] Drop cruft notebook/script file coldonly projection was merged into the other project_effects script. Removing this to clean up the interface. --- scripts/05-project_effects_marimo.py | 678 --------------------------- 1 file changed, 678 deletions(-) delete mode 100644 scripts/05-project_effects_marimo.py diff --git a/scripts/05-project_effects_marimo.py b/scripts/05-project_effects_marimo.py deleted file mode 100644 index de0f6ec..0000000 --- a/scripts/05-project_effects_marimo.py +++ /dev/null @@ -1,678 +0,0 @@ -import marimo - -__generated_with = "0.23.8" -app = marimo.App(width="medium") - - -@app.cell(hide_code=True) -def _(): - import marimo as mo - - mo.md( - """ - # Project mortality effects - - Project mortality effects for forecast ensemble and baseline period. - """ - ) - return - - -@app.cell -def _(): - import datetime - import os - import uuid - - from dotenv import load_dotenv - import isku - import numpy as np - import xarray as xr - - from poreallas.extract import make_climtas, make_tas_monthly_histogram - from poreallas.project import calculate_beta, mortality_effect_model - - load_dotenv() - return ( - calculate_beta, - datetime, - isku, - make_climtas, - make_tas_monthly_histogram, - mortality_effect_model, - np, - os, - uuid, - xr, - ) - - -@app.cell -def _(os): - TAS_FORECAST_URI = os.environ["POREALLAS_TAS_FORECAST_URI"] - ERA5_URI = os.environ["POREALLAS_ERA5_URI"] - GAMMA_URI = os.environ["POREALLAS_GAMMA_URI"] - SOCIOECONOMICS_URI = os.environ["POREALLAS_SOCIOECONOMICS_URI"] - REGIONS_URI = os.environ["POREALLAS_REGIONS_URI"] - EFFECTS_URI = "gs://poreallas-public-20260605/v20260702/effects_net_hot_cold.zarr" # os.getenv("POREALLAS_EFFECTS_URI") - return ( - EFFECTS_URI, - ERA5_URI, - GAMMA_URI, - REGIONS_URI, - SOCIOECONOMICS_URI, - TAS_FORECAST_URI, - ) - - -@app.cell -def _(TAS_FORECAST_URI): - TAS_FORECAST_URI - return - - -@app.cell -def _(isku, np, xr): - # def make_daily_tas_avg(ds: xr.Dataset) -> xr.Dataset: - # _tas = xr.DataArray(units.convert_units_to(ds["tas"], "degC")) - - # _tas_daily_avg = _tas.groupby("time.day").mean(dim="time") - - # return _tas_daily_avg.to_dataset().astype("float32") - - def read_reanalysis(uri: str) -> xr.Dataset: - _ds = xr.load_dataset( - uri, - storage_options={"token": "anon"}, - ) - - # Clean up longitude. The data goes from longitude 0 to 360. It needs to go -180 to 180 in ascending order. - _ds["longitude"] = (_ds["longitude"] + 180) % 360 - 180 - _ds = _ds.sortby("longitude") - _ds = _ds.rename({"longitude": "lon", "latitude": "lat"}) - _ds = _ds.chunk("auto") - - return _ds - - def read_forecast_ensemble(uri: str) -> xr.Dataset: - _ds = xr.load_dataset( - uri, - storage_options={"token": "anon"}, - ) - - # Clean up longitude. The data goes from longitude 0 to 360. It needs to go -180 to 180 in ascending order. - _ds["longitude"] = (_ds["longitude"] + 180) % 360 - 180 - _ds = _ds.sortby("longitude") - _ds = _ds.rename({"latitude": "lat", "longitude": "lon"}) - _ds = _ds.chunk("auto") - - # TODO: We prob don't want this here. Should be in earlier cleaning. Here for backwards compatibility. - # Drop months without required number of obs. Forecast ensemble is for a fixed number of days so we expect to usually trim off the last month of the forecast if it is ragged and missing days beyond a threshold. - _dt_dim = "time" - _n_initial = _ds[_dt_dim].size - _number_obs = _ds[_dt_dim].resample(time="ME").count() - _days_in_month = _number_obs[_dt_dim].dt.days_in_month - required_percent = 0.9 - _min_req = np.round(_days_in_month * required_percent) - _qualifying_months = _number_obs.where(_number_obs >= _min_req, drop=True)[ - "time" - ].dt.month - _ds = _ds.where(_ds[_dt_dim].dt.month.isin(_qualifying_months), drop=True) - - _n_current = _ds[_dt_dim].size - _n_initial_months = _number_obs[_dt_dim].size - _n_qualifying_months = _qualifying_months["time"].size - - print( - f"continuing with {_n_qualifying_months} of {_n_initial_months} forecast months after removing incomplete months" - ) - print( - f"continuing with {_n_current} of {_n_initial} forecast periods after removing incomplete months" - ) - - assert (_n_qualifying_months - _n_initial_months) < 2, ( - "More than one incomplete month was removed from the forecast while checking for incomplete months. Something unexpected is happening." - ) - - return _ds - - def read_regions(uri: str) -> isku.GridWeightingRegions: - _region_weights = xr.load_dataset( - uri, - storage_options={"token": "anon"}, - )[["lat", "lon", "region", "weight"]] # Load only what we need. - # Apparently in this version of xarray the `.load()` method type-hints it'll return a DataArray instead of a Dataset. - # It is a Dataset (I checked). So telling ty to ignore it. - # # TODO: send bug upstream? - regions = isku.GridWeightingRegions(_region_weights) - return regions - - def read_gammas(uri: str) -> xr.Dataset: - return xr.load_dataset( - uri, - storage_options={"token": "anon"}, - ) - - def read_socioeconomics(uri: str) -> xr.Dataset: - return xr.load_dataset( - uri, - storage_options={"token": "anon"}, - ) - - return read_forecast_ensemble, read_reanalysis, read_regions - - -@app.cell -def _( - ERA5_URI, - REGIONS_URI, - TAS_FORECAST_URI, - read_forecast_ensemble, - read_reanalysis, - read_regions, -): - reanalysis = read_reanalysis(ERA5_URI) - forecast_ensemble = read_forecast_ensemble(TAS_FORECAST_URI) - regions = read_regions(REGIONS_URI) - # socioeconomics = read_socioeconomics(SOCIOECONOMICS_URI) - # gammas = read_gammas(GAMMA_URI) - return forecast_ensemble, reanalysis, regions - - -@app.cell -def _(ERA5_URI, TAS_FORECAST_URI, forecast_ensemble, reanalysis, xr): - bias_adjust = True # bias adjust to GMFD - save_adjusted = True # save the bias adjusted reanalysis and forecast data to gcp - if bias_adjust: - # here we want to delta shift both era5 and the forecast to GMFD I think - era5file = "../tas_era5_1deg_monthly_climo_1980-2010.nc" - gmfdfile = "../tas_gmfd_1deg_monthly_climo_1980-2010.nc" - era5 = xr.open_dataset(era5file) - gmfd = xr.open_dataset(gmfdfile) - # need to convert lons to -180 - 180 - era5["longitude"] = (era5["longitude"] + 180) % 360 - 180 - gmfd["longitude"] = (gmfd["longitude"] + 180) % 360 - 180 - - bias = gmfd - era5 - bias = bias.rename({"latitude": "lat", "longitude": "lon"}) - - # bias = bias.assign_coords({"month":reanalysis.time.dt.month}) - # want to match GMFD so if ERA5 is warmer, need to make it cooler. - # This means delta for gmfd-era5 is negative and can add it to era5 to cool era5 down. - # So bias = gmfd - era5, and then adjusted = forecast + bias. And era5 + bias - reanalysis_adjusted = reanalysis + bias.sel(month=reanalysis.time.dt.month) - - forecast_ensemble_adjusted = forecast_ensemble + bias.sel( - month=forecast_ensemble.time.dt.month - ) - - if save_adjusted: - reanalysis_adjusted.attrs.update( - { - "Description": "ERA5 adjusted to GMFD using a monthly climatology difference from ERA5-GMFD over 1980-2010.", - "reanalysis data": ERA5_URI, - "forecast_ensemble_data": TAS_FORECAST_URI, - } - ) - - reanalysis_adjusted.to_zarr("reanalysis_bias_adjusted_to_GMFD.zarr") - - forecast_ensemble_adjusted.attrs.update( - { - "Description": "Seasonal forecast ensemble adjusted to GMFD using a monthly climatology difference from ERA5-GMFD over 1980-2010.", - "reanalysis data": ERA5_URI, - "forecast_ensemble_data": TAS_FORECAST_URI, - "Bias adjustment": "TODO add this step to poreallas. Bias adjustment delta calculated on notebooks.cilresearch.org. https://notebooks.cilresearch.org/user/kemccusker/lab/tree/ClimateImpactLab/bias_adjust_ERA5_to_GMFD_monthly_for_ENSO_work.ipynb", - } - ) - - forecast_ensemble_adjusted.to_zarr( - "forecast_ensemble_bias_adjusted_to_GMFD.zarr" - ) - return bias_adjust, forecast_ensemble_adjusted, reanalysis_adjusted - - -@app.cell -def _(forecast_ensemble_adjusted): - forecast_ensemble_adjusted.tas.isnull().any().values - return - - -@app.cell -def _(xr): - testds = xr.open_zarr("forecast_ensemble_bias_adjusted_to_GMFD.zarr") - return (testds,) - - -@app.cell -def _(testds): - testds.tas.isnull().any().values - return - - -@app.cell -def _( - bias_adjust, - forecast_ensemble, - forecast_ensemble_adjusted, - isku, - make_climtas, - make_tas_monthly_histogram, - np, - reanalysis, - reanalysis_adjusted, - regions, - socioeconomics, - xr, -): - - if bias_adjust: - # Transform _adjusted_ gridded data, extracting regional data needed for projections. - histogram_hist_tas = isku.extract_regions( - reanalysis_adjusted, - template=make_tas_monthly_histogram, - regions=regions, - ) - histogram_forecast_tas = isku.extract_regions( - forecast_ensemble_adjusted, - template=make_tas_monthly_histogram, - regions=regions, - ) - else: - # Transform gridded data, extracting regional data needed for projections. - histogram_hist_tas = isku.extract_regions( - reanalysis, - template=make_tas_monthly_histogram, - regions=regions, - ) - histogram_forecast_tas = isku.extract_regions( - forecast_ensemble, - template=make_tas_monthly_histogram, - regions=regions, - ) - # Using the same static beta for forecast and reanalysis projection requires the histogram tas_bin for these data need to be equal, too. So we're calculating it here. - xr.testing.assert_allclose( - histogram_hist_tas["tas_bin"], - histogram_forecast_tas["tas_bin"], - ) - - # calculate covariates - # average recent climate - climtas = isku.extract_regions( - reanalysis, - template=make_climtas, - regions=regions, - ).sel(year=2025, drop=True) - # income - loggdppc = np.log(socioeconomics["gdppc"].sel(year=2023, drop=True)) - return climtas, histogram_forecast_tas, histogram_hist_tas, loggdppc - - -@app.cell -def _(calculate_beta, climtas, gammas, histogram_forecast_tas, loggdppc, xr): - # Calculate a fixed response function, i.e. beta. - # Single, static response function with no adaptation is used for both projections (does not vary in time). - # Stick everything together and make sure it aligns and matches. Rechunk all together. Also drop any regions with NaNs. - beta_input = ( - xr.Dataset( - { - "tas_bin": histogram_forecast_tas["tas_bin"], - "climtas": climtas["climtas"], - "loggdppc": loggdppc, - "gamma": gammas["gamma_mean"], - } - ) - .dropna(dim="region") # TODO: check these nan values. should they be nan? - .chunk( - { - "region": "auto", # "auto" is a sensible default. - "tas_bin": -1, # This also needs to be all in memory. - "age_cohort": 1, # We're doing all age_cohorts at once but could be done one-by-one. - "degree": -1, # For gammas and polynomial calculations. Should all be in memory. - }, - ) - .unify_chunks() - ) - fixed_beta = calculate_beta(beta_input).astype("float32").compute() - fixed_beta["beta"].attrs = { - "units": "deaths per 100,000 people", - "long_name": "Temperature mortality rate", - } - # Do beta, allowing only hot deaths by 0-ing out everything on the cold side of the minimum-mortality temperature. - fixed_beta["beta_hotonly"] = fixed_beta["beta"].where( - fixed_beta["tas_bin"] > fixed_beta["mmt"], other=0 - ) - fixed_beta["beta_hotonly"].attrs["long_name"] = "Hot temperature mortality rate" - - # Do beta, allowing only COLD deaths by 0-ing out everything on the warm side of the minimum-mortality temperature. - fixed_beta["beta_coldonly"] = fixed_beta["beta"].where( - fixed_beta["tas_bin"] < fixed_beta["mmt"], other=0 - ) - fixed_beta["beta_coldonly"].attrs["long_name"] = "Cold temperature mortality rate" - return (fixed_beta,) - - -@app.cell -def _(fixed_beta, histogram_forecast_tas, isku, mortality_effect_model, xr): - # Project mortality. - # Start with forecast ensemble. - _forecast_input = ( - xr.Dataset( - { - "histogram_tas": histogram_forecast_tas["histogram_tas"], - "beta": fixed_beta["beta"], - } - ) - .dropna(dim="region") # TODO check the nans are expected - .chunk( - { - "region": "auto", # "auto" is a sensible default. - "time": -1, - "tas_bin": -1, - "age_cohort": 1, - "number": 1, - }, - ) - .unify_chunks() - ) - projected_forecast = isku.project( - _forecast_input, model=mortality_effect_model - ).compute() - projected_forecast["effect"].attrs = { - "units": "deaths per 100,000 people", - "long_name": "Temperature mortality", - } - return (projected_forecast,) - - -@app.cell -def _(fixed_beta, histogram_forecast_tas, isku, mortality_effect_model, xr): - # Now hot-only projection - _forecast_input = ( - xr.Dataset( - { - "histogram_tas": histogram_forecast_tas["histogram_tas"], - "beta": fixed_beta["beta_hotonly"], - } - ) - .dropna(dim="region") - .chunk( - { - "region": "auto", # "auto" is a sensible default. - "time": -1, - "tas_bin": -1, - "age_cohort": 1, - "number": 1, - }, - ) - .unify_chunks() - ) - projected_forecast_hotonly = isku.project( - _forecast_input, model=mortality_effect_model - ).compute() - projected_forecast_hotonly["effect"].attrs = { - "units": "deaths per 100,000 people", - "long_name": "Hot temperature mortality", - } - return (projected_forecast_hotonly,) - - -@app.cell -def _(fixed_beta, histogram_forecast_tas, isku, mortality_effect_model, xr): - # Now cold-only projection - _forecast_input = ( - xr.Dataset( - { - "histogram_tas": histogram_forecast_tas["histogram_tas"], - "beta": fixed_beta["beta_coldonly"], - } - ) - .dropna(dim="region") - .chunk( - { - "region": "auto", # "auto" is a sensible default. - "time": -1, - "tas_bin": -1, - "age_cohort": 1, - "number": 1, - }, - ) - .unify_chunks() - ) - projected_forecast_coldonly = isku.project( - _forecast_input, model=mortality_effect_model - ).compute() - projected_forecast_coldonly["effect"].attrs = { - "units": "deaths per 100,000 people", - "long_name": "Cold temperature mortality", - } - return (projected_forecast_coldonly,) - - -@app.cell -def _(fixed_beta, histogram_hist_tas, isku, mortality_effect_model, xr): - # Now do the baseline period. - # Stick everything together and make sure it aligns and matches. Rechunk all together. Also drop any regions with NaNs. - _hist_input = ( - xr.Dataset( - { - "histogram_tas": histogram_hist_tas["histogram_tas"], - "beta": fixed_beta["beta"], - } - ) - .dropna(dim="region") - .chunk( - { - "region": "auto", # "auto" is a sensible default. - "time": -1, - "tas_bin": -1, - "age_cohort": 1, - }, - ) - .unify_chunks() - ) - projected_hist = isku.project(_hist_input, model=mortality_effect_model).compute() - projected_hist["effect"].attrs = { - "units": "deaths per 100,000 people", - "long_name": "Temperature mortality", - } - return (projected_hist,) - - -@app.cell -def _(fixed_beta, histogram_hist_tas, isku, mortality_effect_model, xr): - # Now hot-only projection - _hist_input = ( - xr.Dataset( - { - "histogram_tas": histogram_hist_tas["histogram_tas"], - "beta": fixed_beta["beta_hotonly"], - } - ) - .dropna(dim="region") - .chunk( - { - "region": "auto", # "auto" is a sensible default. - "time": -1, - "tas_bin": -1, - "age_cohort": 1, - }, - ) - .unify_chunks() - ) - projected_hist_hotonly = isku.project( - _hist_input, model=mortality_effect_model - ).compute() - projected_hist_hotonly["effect"].attrs = { - "units": "deaths per 100,000 people", - "long_name": "Hot temperature mortality", - } - return (projected_hist_hotonly,) - - -@app.cell -def _(fixed_beta, histogram_hist_tas, isku, mortality_effect_model, xr): - # Now cold-only projection - _hist_input = ( - xr.Dataset( - { - "histogram_tas": histogram_hist_tas["histogram_tas"], - "beta": fixed_beta["beta_coldonly"], - } - ) - .dropna(dim="region") - .chunk( - { - "region": "auto", # "auto" is a sensible default. - "time": -1, - "tas_bin": -1, - "age_cohort": 1, - }, - ) - .unify_chunks() - ) - projected_hist_coldonly = isku.project( - _hist_input, model=mortality_effect_model - ).compute() - projected_hist_coldonly["effect"].attrs = { - "units": "deaths per 100,000 people", - "long_name": "Cold temperature mortality", - } - return (projected_hist_coldonly,) - - -@app.cell -def _(EFFECTS_URI): - # EFFECTS_URI_LOCAL = "local_output/v20260702/effects_net_hot_cold.zarr" - EFFECTS_URI_LOCAL = ( - "local_output/v20260702/effects_net_hot_cold_gmfd_biasadjusted.zarr" - ) - - print(EFFECTS_URI) - print(EFFECTS_URI_LOCAL) - return (EFFECTS_URI_LOCAL,) - - -@app.cell -def _( - EFFECTS_URI, - EFFECTS_URI_LOCAL, - ERA5_URI, - GAMMA_URI, - REGIONS_URI, - SOCIOECONOMICS_URI, - TAS_FORECAST_URI, - datetime, - projected_forecast, - projected_forecast_coldonly, - projected_forecast_hotonly, - projected_hist, - projected_hist_coldonly, - projected_hist_hotonly, - uuid, - xr, -): - - # Collect everything and write to storage. - _out = { - "forecast": projected_forecast, - "baseline": projected_hist, - "forecast_hotonly": projected_forecast_hotonly, - "baseline_hotonly": projected_hist_hotonly, - "forecast_coldonly": projected_forecast_coldonly, - "baseline_coldonly": projected_hist_coldonly, - } - _out_dt = xr.DataTree.from_dict(_out) - - # Add metadata - _uid = str(uuid.uuid4()) - _datetime_now = datetime.datetime.now(datetime.timezone.utc).isoformat() - - _out_dt.attrs |= { - "poreallas_created_at": _datetime_now, - "poreallas_uid": _uid, - "poreallas_description": "Projected temperature mortality effects", - } - - _out_dt["forecast"].attrs |= { - "poreallas_created_at": _datetime_now, - "poreallas_uid": _uid, - "poreallas_description": "Forecast ensemble projected temperature mortality effects", - "poreallas_temperature_uri": TAS_FORECAST_URI, - "poreallas_socioeconomics_uri": SOCIOECONOMICS_URI, - "poreallas_model_parameters_uri": GAMMA_URI, - "poreallas_regions_uri": REGIONS_URI, - } - - _out_dt["baseline"].attrs |= { - "poreallas_created_at": _datetime_now, - "poreallas_uid": _uid, - "poreallas_description": "Baseline projected temperature mortality effects", - "poreallas_temperature_uri": ERA5_URI, - "poreallas_socioeconomics_uri": SOCIOECONOMICS_URI, - "poreallas_model_parameters_uri": GAMMA_URI, - "poreallas_regions_uri": REGIONS_URI, - } - - _out_dt["forecast_hotonly"].attrs |= { - "poreallas_created_at": _datetime_now, - "poreallas_uid": _uid, - "poreallas_description": "Forecast ensemble projected hot temperature mortality effects", - "poreallas_temperature_uri": TAS_FORECAST_URI, - "poreallas_socioeconomics_uri": SOCIOECONOMICS_URI, - "poreallas_model_parameters_uri": GAMMA_URI, - "poreallas_regions_uri": REGIONS_URI, - } - - _out_dt["baseline_hotonly"].attrs |= { - "poreallas_created_at": _datetime_now, - "poreallas_uid": _uid, - "poreallas_description": "Baseline projected hot temperature mortality effects", - "poreallas_temperature_uri": ERA5_URI, - "poreallas_socioeconomics_uri": SOCIOECONOMICS_URI, - "poreallas_model_parameters_uri": GAMMA_URI, - "poreallas_regions_uri": REGIONS_URI, - } - - _out_dt["forecast_coldonly"].attrs |= { - "poreallas_created_at": _datetime_now, - "poreallas_uid": _uid, - "poreallas_description": "Forecast ensemble projected cold temperature mortality effects", - "poreallas_temperature_uri": TAS_FORECAST_URI, - "poreallas_socioeconomics_uri": SOCIOECONOMICS_URI, - "poreallas_model_parameters_uri": GAMMA_URI, - "poreallas_regions_uri": REGIONS_URI, - } - - _out_dt["baseline_coldonly"].attrs |= { - "poreallas_created_at": _datetime_now, - "poreallas_uid": _uid, - "poreallas_description": "Baseline projected cold temperature mortality effects", - "poreallas_temperature_uri": ERA5_URI, - "poreallas_socioeconomics_uri": SOCIOECONOMICS_URI, - "poreallas_model_parameters_uri": GAMMA_URI, - "poreallas_regions_uri": REGIONS_URI, - } - if EFFECTS_URI_LOCAL is not None: - _out_dt.to_zarr(EFFECTS_URI_LOCAL, consolidated=False) - print(f"Effects written to {EFFECTS_URI_LOCAL}") - if EFFECTS_URI is not None: - try: - _out_dt.to_zarr(EFFECTS_URI, consolidated=False) - print(f"Effects written to {EFFECTS_URI}") - except Exception: - print( - "Caught Exception. You probably don't have permissions to write to the bucket" - ) - - _out_dt - return - - -@app.cell -def _(): - return - - -if __name__ == "__main__": - app.run() From c25016edb622a6b5d3f0996998366da7541aa944 Mon Sep 17 00:00:00 2001 From: Brewster Malevich Date: Mon, 3 Aug 2026 20:31:14 -0700 Subject: [PATCH 3/5] Tell CI to ingore typing warnings and import errors from cluster-run scripts --- scripts/01-parse_era5.py | 4 ++-- scripts/01-parse_forecast.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/scripts/01-parse_era5.py b/scripts/01-parse_era5.py index e45c30f..cff55ec 100644 --- a/scripts/01-parse_era5.py +++ b/scripts/01-parse_era5.py @@ -16,10 +16,10 @@ import uuid import dask -from dask_gateway import GatewayCluster +from dask_gateway import GatewayCluster # type: ignore[ty:unresolved-import] from dotenv import load_dotenv import xarray as xr -import xesmf as xe +import xesmf as xe # type: ignore[ty:unresolved-import] load_dotenv() diff --git a/scripts/01-parse_forecast.py b/scripts/01-parse_forecast.py index d38a49c..aa39e8a 100644 --- a/scripts/01-parse_forecast.py +++ b/scripts/01-parse_forecast.py @@ -45,8 +45,8 @@ compat="no_conflicts", ) s51 = ( - xr.open_mfdataset(target_tasmax_paths, **mfdataset_kwargs)["mx2t24"] - + xr.open_mfdataset(target_tasmin_paths, **mfdataset_kwargs)["mn2t24"] + xr.open_mfdataset(target_tasmax_paths, **mfdataset_kwargs)["mx2t24"] # type: ignore[ty:invalid-argument-type] + + xr.open_mfdataset(target_tasmin_paths, **mfdataset_kwargs)["mn2t24"] # type: ignore[ty:invalid-argument-type] ) / 2 s51.name = "tas" s51 = s51.to_dataset() From a5470ed055bf1de42e6826897de4b535b4eac245 Mon Sep 17 00:00:00 2001 From: Brewster Malevich Date: Mon, 3 Aug 2026 20:33:51 -0700 Subject: [PATCH 4/5] Tell CI to ignore typiing warnings and import errors from cluster-run scripts --- scripts/01-parse_gmfd.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/01-parse_gmfd.py b/scripts/01-parse_gmfd.py index 35f5cf4..9fd8511 100644 --- a/scripts/01-parse_gmfd.py +++ b/scripts/01-parse_gmfd.py @@ -11,10 +11,10 @@ import uuid import dask -from dask_gateway import GatewayCluster +from dask_gateway import GatewayCluster # type: ignore[ty:unresolved-import] from dotenv import load_dotenv import xarray as xr -import xesmf as xe +import xesmf as xe # type: ignore[ty:unresolved-import] load_dotenv() From e018cb986c2965cc7c14ed536833ff695a49c54b Mon Sep 17 00:00:00 2001 From: Brewster Malevich Date: Tue, 4 Aug 2026 10:54:11 -0700 Subject: [PATCH 5/5] Add CHANGELOG entry for changes --- CHANGELOG.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index fe027aa..53b8cfa 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Changes + +- 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.2.0] - 2026-07-27 ### Added