#!/usr/bin/env python3 """Map decadal linear trends in July FATES aboveground carbon stocks. Run: python plot_fates_decadal_maps.py Use --start-year 1985 to fit only data from 1985 onwards. Uses the dependencies and HydroBASINS cache of plot_fates_basin_trends.py. Reads one monthly h0a file per year, once for both basins. Ordinary least squares slopes through July stocks are converted from kg C m-2 yr-1 to Mg C ha-1 yr-1 by multiplying by 10. This is an aboveground carbon-stock change proxy, NOT net ecosystem uptake or net primary production. Each basin figure has one panel per calendar decade. Partial decades are included and labelled with their actual years. All panels in both figures share one symmetric colour scale, defaulting to -1 to +1 Mg C ha-1 yr-1. Set --color-limit to change its positive magnitude. Values outside the limits saturate at the end colours; colour-bar arrows indicate saturation. NetCDF values remain unmodified. Positive slopes mean carbon gains. Basin membership uses land grid-cell centres inside HydroBASINS v1c level-3 watersheds (Lehner & Grill, 2013, https://doi.org/10.1002/hyp.9740). Cells missing any July stock within a decade are excluded from that decade. Slopes and basin masks are also saved in a NetCDF file. """ import argparse import calendar from pathlib import Path import warnings import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.colors import Normalize from matplotlib.ticker import FuncFormatter import numpy as np import shapely from shapely.plotting import plot_polygon import xarray as xr from plot_fates_basin_trends import ( BASINS, DATA_DIR, RUN_LABEL, VARIABLE, find_files, load_basin, ) def linear_trend(years: np.ndarray, stocks: np.ndarray) -> np.ndarray: """OLS slope in Mg C ha-1 yr-1, requiring a complete time series.""" if len(years) < 2 or stocks.shape[0] != len(years): raise ValueError("A trend needs at least two matching years and stock fields") centred_years = years.astype(float) - np.mean(years) denominator = np.sum(centred_years ** 2) if denominator == 0: raise ValueError("Trend years must span a nonzero interval") return np.einsum("t,tij->ij", centred_years, stocks) / denominator * 10 def read_fields(files: dict[int, Path], basins: dict): fields = [] masks = {} reference_lat = reference_lon = None order = None for year, path in files.items(): with xr.open_dataset( path, engine="netcdf4", decode_times=False, cache=False ) as dataset: variable = dataset[VARIABLE] if set(variable.dims) != {"time", "lat", "lon"}: raise ValueError(f"Unexpected dimensions in {path}: {variable.dims}") if variable.sizes["time"] != 1: raise ValueError(f"Expected one monthly record in {path}") if variable.attrs.get("units") != "kg m-2": raise ValueError(f"Unexpected carbon units in {path}: {variable.attrs}") lat = dataset["lat"].values lon = dataset["lon"].values if reference_lat is None: reference_lat, reference_lon = lat.copy(), lon.copy() signed_lon = (lon + 180) % 360 - 180 order = np.argsort(signed_lon) map_lon = signed_lon[order] lon2d, lat2d = np.meshgrid(map_lon, lat) points = shapely.points(lon2d, lat2d) landfrac = dataset["landfrac"].transpose("lat", "lon").values[:, order] land = np.isfinite(landfrac) & (landfrac > 0) & (landfrac <= 1) masks = { name: shapely.covers(geometry, points) & land for name, geometry in basins.items() } for name, mask in masks.items(): if not np.any(mask): raise ValueError(f"No land grid-cell centres inside {name}") elif not ( np.array_equal(lat, reference_lat) and np.array_equal(lon, reference_lon) ): raise ValueError(f"Grid changed in {path}") values = variable.isel(time=0).transpose("lat", "lon").values[:, order] fields.append(values.astype(float)) print(f"Read {path.name}", flush=True) return reference_lat, map_lon, np.stack(fields), masks def main() -> None: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--data-dir", type=Path, default=DATA_DIR) parser.add_argument("--month", type=int, choices=range(1, 13), default=7) parser.add_argument("--start-year", type=int, help="First stock year to include") parser.add_argument("--output-dir", type=Path, default=Path.cwd()) parser.add_argument( "--color-limit", type=float, default=1.0, help="Positive magnitude of shared symmetric colour limits (Mg C ha-1 yr-1)", ) parser.add_argument( "--cache-dir", type=Path, default=Path(__file__).resolve().parent / "hydrobasins_cache", ) args = parser.parse_args() if not np.isfinite(args.color_limit) or args.color_limit <= 0: parser.error("--color-limit must be finite and greater than zero") files = find_files(args.data_dir, args.month) if args.start_year is not None: files = {year: path for year, path in files.items() if year >= args.start_year} if len(files) < 2: parser.error("--start-year must leave at least two available stock years") years = np.array(list(files)) args.cache_dir.mkdir(parents=True, exist_ok=True) basins = { name: load_basin(args.cache_dir, continent, basin_id) for name, (continent, basin_id, _) in BASINS.items() } lat, lon, stocks, masks = read_fields(files, basins) union = np.logical_or.reduce(list(masks.values())) trends, starts, ends, decades = [], [], [], [] for decade in sorted(set(years // 10 * 10)): selected = (years >= decade) & (years <= decade + 9) if np.count_nonzero(selected) < 2: raise ValueError(f"Insufficient samples to fit decade {decade}") decade_stocks = stocks[selected] valid = np.all(np.isfinite(decade_stocks), axis=0) for name, mask in masks.items(): if not np.any(valid & mask): raise ValueError(f"No complete stock series for {name}, {decade}s") missing = np.count_nonzero(mask & ~valid) if missing: warnings.warn( f"{name}, {decade}s: excluding {missing} cells with missing stocks", stacklevel=2, ) finite_stocks = np.where(np.isfinite(decade_stocks), decade_stocks, np.nan) slope = linear_trend(years[selected], finite_stocks) trends.append(np.where(valid & union, slope, np.nan)) starts.append(int(years[selected][0])) ends.append(int(years[selected][-1])) decades.append(int(decade)) trends = np.stack(trends) limit = args.color_limit norm = Normalize(vmin=-limit, vmax=limit) below = np.any(trends < -limit) above = np.any(trends > limit) extend = "both" if below and above else "min" if below else "max" if above else "neither" args.output_dir.mkdir(parents=True, exist_ok=True) stem = f"fates_decadal_stock_change_{args.month:02d}_{years[0]}_{years[-1]}" output = xr.Dataset( { "aboveground_carbon_trend": ( ("decade", "lat", "lon"), trends, { "units": "Mg C ha-1 yr-1", "long_name": f"OLS trend in {calendar.month_name[args.month]} aboveground carbon stocks", }, ), "first_year": ("decade", starts), "last_year": ("decade", ends), **{ f"{name.lower()}_mask": (("lat", "lon"), mask.astype("int8")) for name, mask in masks.items() }, }, coords={ "decade": decades, "lat": ("lat", lat, {"units": "degrees_north"}), "lon": ("lon", lon, {"units": "degrees_east"}), }, attrs={ "method": f"OLS trend through {calendar.month_name[args.month]} monthly mean stocks; kg/m2 to Mg/ha factor 10", "interpretation": "Aboveground carbon-stock change, not net ecosystem uptake", "basin_source": "HydroBASINS v1c level 3; Lehner & Grill (2013), doi:10.1002/hyp.9740", "source_run": ( "ihist.f09_g17.ctsm5.4.042_noresm_v3_Trendy26." "T26_S2_pt2.2026-08-22" ), }, ) nc_path = args.output_dir / f"{stem}.nc" output.to_netcdf(nc_path) print(f"Saved {nc_path}") for name, geometry in basins.items(): ncols = 4 nrows = (len(decades) + ncols - 1) // ncols fig, axes = plt.subplots( nrows, ncols, figsize=(4 * ncols, 4.3 * nrows), squeeze=False, layout="constrained", ) xmin, ymin, xmax, ymax = geometry.bounds for index, ax in enumerate(axes.flat): if index >= len(decades): ax.set_visible(False) continue image = ax.pcolormesh( lon, lat, np.ma.masked_invalid(np.where(masks[name], trends[index], np.nan)), cmap="BrBG", norm=norm, shading="nearest", rasterized=True, ) plot_polygon( geometry, ax=ax, add_points=False, facecolor="none", edgecolor="black", linewidth=0.7, ) partial = starts[index] != decades[index] or ends[index] != decades[index] + 9 ax.set( xlim=(xmin - 1, xmax + 1), ylim=(ymin - 1, ymax + 1), title=f"{starts[index]}-{ends[index]}" + (" (partial)" if partial else ""), xlabel="Longitude", ylabel="Latitude", ) ax.set_aspect(1 / np.cos(np.deg2rad((ymin + ymax) / 2))) ax.xaxis.set_major_formatter(FuncFormatter( lambda value, _: f"{abs(value):g}" + r"$^\circ$" + ("W" if value < 0 else "E") )) ax.yaxis.set_major_formatter(FuncFormatter( lambda value, _: f"{abs(value):g}" + r"$^\circ$" + ("S" if value < 0 else "N") )) ax.grid(alpha=0.2) fig.colorbar( image, ax=list(axes.flat), shrink=0.8, extend=extend, label=r"Aboveground carbon-stock trend (Mg C ha$^{-1}$ yr$^{-1}$)", ) fig.suptitle( f"{name}: {RUN_LABEL} decadal trends in " f"{calendar.month_name[args.month]} aboveground carbon stocks\n" "Positive = gain; negative = loss | HydroBASINS v1c (Lehner & Grill, 2013)", fontsize=14, ) png_path = args.output_dir / f"{stem}_{name.lower()}.png" fig.savefig(png_path, dpi=180) plt.close(fig) print(f"Saved {png_path}") if __name__ == "__main__": main()