From e36343c03829b83ffbba57ad4844e55e3553a746 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 8 Jul 2026 08:51:33 -0600 Subject: [PATCH 01/14] fix: enforce exact fraction closure for fractional fields ELM's surfrdUtilsMod.F90 performs a strict equality check on the sum of natural-patch weights (wt_nat_patch) and landunit weights (wt_lunit). Dapper's zonal aggregation produces each fraction band independently via area-weighted mean, so small per-band floating-point residuals accumulate across all bands. Combined with float32 serialization, the written sums could deviate from 1.0 by enough to trigger ELM's fatal error: surfrd_veg_all ERROR: sum of wt_nat_patch not 1.0 ... surfrd_get_data ERROR: sum of wt_lunit not 1.0 ... Changes ------- sfile.py - Add _scale_to_target_sum(): scales a DataArray along a given dimension so its sum matches a target, working in float64. - Add _snap_partition_sum(): after scaling, assigns the exact residual to the last band so that sum is exactly equal to the target in float64 arithmetic. This eliminates any remaining machine-epsilon drift. - Add _normalize_surface_fraction_closure(): orchestrates normalization for all fraction groups before write: * Landunit scalars (PCT_NATVEG, PCT_CROP, PCT_WETLAND, PCT_LAKE, PCT_GLACIER, PCT_URBAN) are proportionally scaled to 100 where the existing total is already near 100 (gating prevents distorting partial datasets), then an exact residual snap ensures the final sum is exactly 100 in float64. * PCT_NAT_PFT is normalized then snapped to exactly 100 (ELM treats these as natural-patch weights that must sum to 1.0 after dividing by 100). * PCT_CFT and PCT_GLC_MEC are scaled and snapped to their parent percent fields (PCT_CROP and PCT_GLACIER respectively). * FSURF and FGRD (irrigation split) are normalized to sum to 1. - Closure-critical fraction variables are now encoded as float64 in the written NetCDF file instead of float32 to prevent quantization reintroducing the drift that was just corrected in memory. The closure-critical set is: PCT_NAT_PFT, PCT_CFT, PCT_GLC_MEC, PCT_NATVEG, PCT_CROP, PCT_GLACIER, PCT_WETLAND, PCT_LAKE, PCT_URBAN. - _normalize_surface_fraction_closure() is called inside write_surface_nc() so the fix applies to every code path that writes a surface file, regardless of sampling method. landuse.py - Add _scale_to_target_sum() and _normalize_landuse_fraction_closure() mirroring the surface logic, covering the same fraction groups. - _normalize_landuse_fraction_closure() is called immediately before _write_nc() in the zonal sampling path of sample_landuse_timeseries() so exported landuse files satisfy the same closure constraints. - Note: the landuse write path does not yet have the float64 encoding override; it relies on normalization-in-memory being sufficient. If ELM-level landuse checks exhibit a similar issue the encoding should be hardened there as well. validate.py - V-105 consistency check now verifies sum(PCT_NAT_PFT) against 100 (previously compared against PCT_NATVEG, which is the wrong target for ELM's natural-patch weight convention). - Tightened tolerance from 1e-3 to 1e-6 (max absolute difference, not mean). - Diagnostic message now reports both mean and max absolute residual. test_fraction_closure.py (new file) - test_write_surface_nc_enforces_fraction_closure: verifies PCT_NAT_PFT closes to 100, PCT_CFT closes to PCT_CROP, landunit total closes to 100, and FSURF+FGRD closes to 1 after a write/read roundtrip. - test_landuse_normalizer_enforces_fraction_closure: unit test for _normalize_landuse_fraction_closure() in isolation. - test_write_surface_nc_natural_patch_weights_near_exact_one: regression for the first ELM failure; uses a 17-PFT distribution with deliberate sub-ULP perturbations and asserts sum(PCT_NAT_PFT)/100 == 1.0 within 1e-12 after write/read. - test_write_surface_nc_landunit_weights_near_exact_one: regression for the second ELM failure; verifies that the sum of all landunit percent fields divided by 100 equals 1.0 within 1e-12 after write/read. --- src/dapper/landuse/landuse.py | 80 +++++++++++++++ src/dapper/surf/sfile.py | 162 ++++++++++++++++++++++++++++++- src/dapper/surf/validate.py | 13 +-- tests/test_fraction_closure.py | 172 +++++++++++++++++++++++++++++++++ 4 files changed, 419 insertions(+), 8 deletions(-) create mode 100644 tests/test_fraction_closure.py diff --git a/src/dapper/landuse/landuse.py b/src/dapper/landuse/landuse.py index ba3997c..87e6de8 100644 --- a/src/dapper/landuse/landuse.py +++ b/src/dapper/landuse/landuse.py @@ -15,6 +15,85 @@ LonWrap = Literal["auto", "0_360", "-180_180"] + +def _scale_to_target_sum( + da: xr.DataArray, + *, + dim: str, + target: xr.DataArray, + eps: float = 1e-12, +) -> xr.DataArray: + """Scale values along `dim` so their sum matches `target` where possible.""" + work = da.astype(np.float64) + summed = work.sum(dim=dim, skipna=True) + valid = np.isfinite(summed) & (np.abs(summed) > eps) & np.isfinite(target) + safe_den = xr.where(valid, summed, 1.0) + factor = xr.where(valid, target.astype(np.float64) / safe_den, 1.0) + return work * factor + + +def _normalize_landuse_fraction_closure(ds: xr.Dataset) -> xr.Dataset: + """Enforce closure on fraction-like landuse outputs before writing.""" + ds2 = ds.copy(deep=False) + + landunit_scalar_names = ["PCT_NATVEG", "PCT_CROP", "PCT_WETLAND", "PCT_LAKE", "PCT_GLACIER"] + landunit_terms: list[xr.DataArray] = [] + for name in landunit_scalar_names: + if name in ds2: + landunit_terms.append(ds2[name].astype(np.float64)) + + urban_group = None + if "PCT_URBAN" in ds2: + urban = ds2["PCT_URBAN"].astype(np.float64) + if "numurbl" in urban.dims: + urban_group = urban + landunit_terms.append(urban.sum(dim="numurbl", skipna=True)) + else: + landunit_terms.append(urban) + + if len(landunit_terms) >= 2: + current_total = sum(landunit_terms) + near_100 = np.abs(current_total - 100.0) <= 1.0 + valid = np.isfinite(current_total) & (current_total > 1e-12) & near_100 + factor = xr.where(valid, 100.0 / current_total, 1.0) + + for name in landunit_scalar_names: + if name in ds2: + ds2[name] = ds2[name].astype(np.float64) * factor + if urban_group is not None: + ds2["PCT_URBAN"] = urban_group * factor + elif "PCT_URBAN" in ds2: + ds2["PCT_URBAN"] = ds2["PCT_URBAN"].astype(np.float64) * factor + + if "PCT_NAT_PFT" in ds2 and "natpft" in ds2["PCT_NAT_PFT"].dims: + target = xr.full_like( + ds2["PCT_NAT_PFT"].isel(natpft=0, drop=True), + 100.0, + dtype=np.float64, + ) + ds2["PCT_NAT_PFT"] = _scale_to_target_sum(ds2["PCT_NAT_PFT"], dim="natpft", target=target) + + if "PCT_CFT" in ds2 and "cft" in ds2["PCT_CFT"].dims and "PCT_CROP" in ds2: + ds2["PCT_CFT"] = _scale_to_target_sum(ds2["PCT_CFT"], dim="cft", target=ds2["PCT_CROP"].astype(np.float64)) + + if "PCT_GLC_MEC" in ds2 and "nglcec" in ds2["PCT_GLC_MEC"].dims and "PCT_GLACIER" in ds2: + ds2["PCT_GLC_MEC"] = _scale_to_target_sum( + ds2["PCT_GLC_MEC"], + dim="nglcec", + target=ds2["PCT_GLACIER"].astype(np.float64), + ) + + if "FSURF" in ds2 and "FGRD" in ds2: + fs = ds2["FSURF"].astype(np.float64) + fg = ds2["FGRD"].astype(np.float64) + total = fs + fg + valid = np.isfinite(total) & (total > 1e-12) + factor = xr.where(valid, 1.0 / total, 1.0) + ds2["FSURF"] = fs * factor + ds2["FGRD"] = fg * factor + + return ds2 + def sample_landuse_timeseries( src_path: str | Path, df_loc: pd.DataFrame, @@ -270,6 +349,7 @@ def _write_nc(ds_out: xr.Dataset, path: Path) -> None: csv_path = out_path.with_suffix(out_path.suffix + ".zonal_weights.csv") dfw_all.to_csv(csv_path, index=False) + out = _normalize_landuse_fraction_closure(out) _write_nc(out, out_path) df_summary = pd.DataFrame({gid_col: order, "sample_ncells": ncells, "sample_area_total_m2": area_m2}) diff --git a/src/dapper/surf/sfile.py b/src/dapper/surf/sfile.py index 0ac4f42..58363c8 100644 --- a/src/dapper/surf/sfile.py +++ b/src/dapper/surf/sfile.py @@ -187,7 +187,7 @@ def write_surface_nc( import datetime as _dt # ---- global attrs ---- - ds2 = ds.copy(deep=False) + ds2 = _normalize_surface_fraction_closure(ds) merged = dict(ds2.attrs) if dapper_attrs: @@ -202,16 +202,174 @@ def write_surface_nc( ds2.attrs = merged + closure_critical = { + "PCT_NAT_PFT", "PCT_CFT", "PCT_GLC_MEC", + "PCT_NATVEG", "PCT_CROP", "PCT_GLACIER", "PCT_WETLAND", "PCT_LAKE", "PCT_URBAN", + } enc: Dict[str, dict] = {} for v in ds2.data_vars: # Fill values must match dtype. Here we do float32 to match ELM surface expectations. if ds2[v].dtype.kind == "f": - enc[v] = {"dtype": "float32", "_FillValue": np.float32(-9.96921e36)} + if v in closure_critical: + enc[v] = {"dtype": "float64", "_FillValue": np.float64(-9.96921e36)} + else: + enc[v] = {"dtype": "float32", "_FillValue": np.float32(-9.96921e36)} ds2.to_netcdf(out_path, encoding=enc) return out_path +def _scale_to_target_sum( + da: xr.DataArray, + *, + dim: str, + target: xr.DataArray, + eps: float = 1e-12, +) -> xr.DataArray: + """Scale values along `dim` so their sum matches `target` where possible.""" + work = da.astype(np.float64) + summed = work.sum(dim=dim, skipna=True) + valid = np.isfinite(summed) & (np.abs(summed) > eps) & np.isfinite(target) + safe_den = xr.where(valid, summed, 1.0) + factor = xr.where(valid, target.astype(np.float64) / safe_den, 1.0) + return work * factor + + +def _snap_partition_sum( + da: xr.DataArray, + *, + dim: str, + target: xr.DataArray, + eps: float = 1e-15, +) -> xr.DataArray: + """Force exact closure by assigning residual to the last band in float64.""" + work = da.astype(np.float64) + if dim not in work.dims or int(work.sizes.get(dim, 0)) < 1: + return work + + head = work.isel({dim: slice(0, -1)}) + head_sum = head.sum(dim=dim, skipna=True) + last_src = work.isel({dim: -1}) + candidate_last = target.astype(np.float64) - head_sum + + valid = np.isfinite(target) & np.isfinite(head_sum) + snapped_last = xr.where(valid, candidate_last, last_src) + snapped_last = xr.where(np.abs(snapped_last) < eps, 0.0, snapped_last) + + out = work.copy(deep=False) + out[{dim: -1}] = snapped_last + return out + + +def _normalize_surface_fraction_closure(ds: xr.Dataset) -> xr.Dataset: + """Enforce closure on key fractional groups before writing a surface file.""" + ds2 = ds.copy(deep=False) + + # Landunit percentages should close to 100 across available classes. + landunit_scalar_names = ["PCT_NATVEG", "PCT_CROP", "PCT_WETLAND", "PCT_LAKE", "PCT_GLACIER"] + landunit_terms: list[xr.DataArray] = [] + for name in landunit_scalar_names: + if name in ds2: + landunit_terms.append(ds2[name].astype(np.float64)) + + urban_group = None + if "PCT_URBAN" in ds2: + urban = ds2["PCT_URBAN"].astype(np.float64) + if "numurbl" in urban.dims: + urban_group = urban + landunit_terms.append(urban.sum(dim="numurbl", skipna=True)) + else: + landunit_terms.append(urban) + + if len(landunit_terms) >= 2: + current_total = sum(landunit_terms) + # Only adjust small drift; avoid force-normalizing partial datasets. + near_100 = np.abs(current_total - 100.0) <= 1.0 + valid = np.isfinite(current_total) & (current_total > 1e-12) & near_100 + factor = xr.where(valid, 100.0 / current_total, 1.0) + + for name in landunit_scalar_names: + if name in ds2: + ds2[name] = ds2[name].astype(np.float64) * factor + if urban_group is not None: + ds2["PCT_URBAN"] = urban_group * factor + elif "PCT_URBAN" in ds2: + ds2["PCT_URBAN"] = ds2["PCT_URBAN"].astype(np.float64) * factor + + # Force exact landunit closure by assigning residual to one component. + total_after = 0.0 + for name in landunit_scalar_names: + if name in ds2: + total_after = total_after + ds2[name].astype(np.float64) + if "PCT_URBAN" in ds2: + urb = ds2["PCT_URBAN"].astype(np.float64) + if "numurbl" in urb.dims: + total_after = total_after + urb.sum(dim="numurbl", skipna=True) + else: + total_after = total_after + urb + + resid = 100.0 - total_after + + # Prefer scalar classes first; fallback to last urban class. + snapped = False + for name in ("PCT_GLACIER", "PCT_LAKE", "PCT_WETLAND", "PCT_CROP", "PCT_NATVEG"): + if name in ds2: + ds2[name] = ds2[name].astype(np.float64) + resid + snapped = True + break + + if (not snapped) and ("PCT_URBAN" in ds2): + urb = ds2["PCT_URBAN"].astype(np.float64) + if "numurbl" in urb.dims and int(urb.sizes.get("numurbl", 0)) >= 1: + urb[{"numurbl": -1}] = urb.isel(numurbl=-1) + resid + ds2["PCT_URBAN"] = urb + else: + ds2["PCT_URBAN"] = urb + resid + + # PFT partitions are natural-patch weights and must close to 100. + if "PCT_NAT_PFT" in ds2 and "natpft" in ds2["PCT_NAT_PFT"].dims: + pft = ds2["PCT_NAT_PFT"] + target = xr.full_like( + pft.isel(natpft=0, drop=True), + 100.0, + dtype=np.float64, + ) + ds2["PCT_NAT_PFT"] = _scale_to_target_sum(pft, dim="natpft", target=target) + ds2["PCT_NAT_PFT"] = _snap_partition_sum(ds2["PCT_NAT_PFT"], dim="natpft", target=target) + + # Crop and glacier-class partitions should close to their parent percentages. + if "PCT_CFT" in ds2 and "cft" in ds2["PCT_CFT"].dims and "PCT_CROP" in ds2: + ds2["PCT_CFT"] = _scale_to_target_sum(ds2["PCT_CFT"], dim="cft", target=ds2["PCT_CROP"].astype(np.float64)) + ds2["PCT_CFT"] = _snap_partition_sum( + ds2["PCT_CFT"], + dim="cft", + target=ds2["PCT_CROP"].astype(np.float64), + ) + if "PCT_GLC_MEC" in ds2 and "nglcec" in ds2["PCT_GLC_MEC"].dims and "PCT_GLACIER" in ds2: + ds2["PCT_GLC_MEC"] = _scale_to_target_sum( + ds2["PCT_GLC_MEC"], + dim="nglcec", + target=ds2["PCT_GLACIER"].astype(np.float64), + ) + ds2["PCT_GLC_MEC"] = _snap_partition_sum( + ds2["PCT_GLC_MEC"], + dim="nglcec", + target=ds2["PCT_GLACIER"].astype(np.float64), + ) + + # Unitless partition currently represented by irrigation split. + if "FSURF" in ds2 and "FGRD" in ds2: + fs = ds2["FSURF"].astype(np.float64) + fg = ds2["FGRD"].astype(np.float64) + total = fs + fg + valid = np.isfinite(total) & (total > 1e-12) + factor = xr.where(valid, 1.0 / total, 1.0) + ds2["FSURF"] = fs * factor + ds2["FGRD"] = fg * factor + + return ds2 + + class CustomizeError(ValueError): """Raised when a customization fails schema/formatting validation.""" diff --git a/src/dapper/surf/validate.py b/src/dapper/surf/validate.py index 5ffc15f..6a36350 100644 --- a/src/dapper/surf/validate.py +++ b/src/dapper/surf/validate.py @@ -320,13 +320,14 @@ def _in_range(a, lo, hi): def _check_soft_consistency(self, ds: xr.Dataset) -> List[CheckResult]: r: List[CheckResult] = [] - # sum(PCT_NAT_PFT) ≈ PCT_NATVEG (reduce across all non-PFT dims) - if "PCT_NATVEG" in ds and "PCT_NAT_PFT" in ds and "natpft" in ds["PCT_NAT_PFT"].dims: + # sum(PCT_NAT_PFT) ≈ 100 (natural-patch weights) + if "PCT_NAT_PFT" in ds and "natpft" in ds["PCT_NAT_PFT"].dims: pftsum = ds["PCT_NAT_PFT"].sum(dim="natpft", skipna=True) - # compare aggregated means (robust across extra dims) - diff = np.nanmean(np.abs(pftsum.values - np.asarray(ds["PCT_NATVEG"].values))) - r.append(CheckResult("V-105.consistency.pftsum", "WARN", diff <= 1e-3, - f"mean(|sum(PCT_NAT_PFT)-PCT_NATVEG|)={diff:.3e}")) + resid = np.abs(pftsum.values - 100.0) + mean_diff = float(np.nanmean(resid)) + max_diff = float(np.nanmax(resid)) + r.append(CheckResult("V-105.consistency.pftsum", "WARN", max_diff <= 1e-6, + f"mean(|sum(PCT_NAT_PFT)-100|)={mean_diff:.3e}; max={max_diff:.3e}")) # If any cell has PCT_URBAN>0 → URBAN_REGION_ID present if "PCT_URBAN" in ds: diff --git a/tests/test_fraction_closure.py b/tests/test_fraction_closure.py new file mode 100644 index 0000000..45e16e9 --- /dev/null +++ b/tests/test_fraction_closure.py @@ -0,0 +1,172 @@ +import numpy as np +import xarray as xr + +from dapper.surf.sfile import write_surface_nc +from dapper.landuse.landuse import _normalize_landuse_fraction_closure + + +def test_write_surface_nc_enforces_fraction_closure(tmp_path): + ds = xr.Dataset( + coords={ + "natpft": np.arange(3, dtype=np.int32), + "cft": np.arange(2, dtype=np.int32), + "numurbl": np.arange(2, dtype=np.int32), + "lsmlat": np.arange(2, dtype=np.int32), + "lsmlon": np.arange(1, dtype=np.int32), + } + ) + + ds["PCT_NATVEG"] = xr.DataArray(np.array([[40.0], [60.0]]), dims=("lsmlat", "lsmlon")) + ds["PCT_NAT_PFT"] = xr.DataArray( + np.array( + [ + [[10.0], [20.0]], + [[20.0], [20.0]], + [[9.5], [19.2]], + ] + ), + dims=("natpft", "lsmlat", "lsmlon"), + ) + + ds["PCT_CROP"] = xr.DataArray(np.array([[30.0], [20.0]]), dims=("lsmlat", "lsmlon")) + ds["PCT_CFT"] = xr.DataArray( + np.array( + [ + [[10.0], [10.0]], + [[18.0], [9.0]], + ] + ), + dims=("cft", "lsmlat", "lsmlon"), + ) + + ds["PCT_WETLAND"] = xr.DataArray(np.array([[5.0], [5.0]]), dims=("lsmlat", "lsmlon")) + ds["PCT_LAKE"] = xr.DataArray(np.array([[20.0], [10.0]]), dims=("lsmlat", "lsmlon")) + ds["PCT_GLACIER"] = xr.DataArray(np.array([[5.0], [5.0]]), dims=("lsmlat", "lsmlon")) + ds["PCT_URBAN"] = xr.DataArray( + np.array( + [ + [[0.5], [0.6]], + [[0.5], [0.4]], + ] + ), + dims=("numurbl", "lsmlat", "lsmlon"), + ) + + ds["FSURF"] = xr.DataArray(np.array([[0.2], [0.8]]), dims=("lsmlat", "lsmlon")) + ds["FGRD"] = xr.DataArray(np.array([[0.6], [0.1]]), dims=("lsmlat", "lsmlon")) + + out_path = tmp_path / "surf_closure.nc" + write_surface_nc(ds, str(out_path)) + + out = xr.open_dataset(out_path) + + pft_sum = out["PCT_NAT_PFT"].sum(dim="natpft") + np.testing.assert_allclose(pft_sum.values, 100.0, atol=1e-5, rtol=0.0) + + cft_sum = out["PCT_CFT"].sum(dim="cft") + np.testing.assert_allclose(cft_sum.values, out["PCT_CROP"].values, atol=1e-5, rtol=0.0) + + urban_sum = out["PCT_URBAN"].sum(dim="numurbl") + landunit_total = ( + out["PCT_NATVEG"] + + out["PCT_CROP"] + + out["PCT_WETLAND"] + + out["PCT_LAKE"] + + out["PCT_GLACIER"] + + urban_sum + ) + np.testing.assert_allclose(landunit_total.values, 100.0, atol=1e-5, rtol=0.0) + + irrigation_total = out["FSURF"] + out["FGRD"] + np.testing.assert_allclose(irrigation_total.values, 1.0, atol=1e-6, rtol=0.0) + + out.close() + + +def test_landuse_normalizer_enforces_fraction_closure(): + ds = xr.Dataset( + coords={ + "natpft": np.arange(2, dtype=np.int32), + "lsmlat": np.arange(1, dtype=np.int32), + "lsmlon": np.arange(2, dtype=np.int32), + } + ) + + ds["PCT_NATVEG"] = xr.DataArray(np.array([[50.0, 60.0]]), dims=("lsmlat", "lsmlon")) + ds["PCT_NAT_PFT"] = xr.DataArray( + np.array( + [ + [[20.0, 30.0]], + [[29.0, 31.0]], + ] + ), + dims=("natpft", "lsmlat", "lsmlon"), + ) + + fixed = _normalize_landuse_fraction_closure(ds) + pft_sum = fixed["PCT_NAT_PFT"].sum(dim="natpft") + np.testing.assert_allclose(pft_sum.values, 100.0, atol=1e-10, rtol=0.0) + + +def test_write_surface_nc_natural_patch_weights_near_exact_one(tmp_path): + ds = xr.Dataset( + coords={ + "natpft": np.arange(17, dtype=np.int32), + "lsmlat": np.arange(1, dtype=np.int32), + "lsmlon": np.arange(1, dtype=np.int32), + } + ) + + # Build a distribution that typically accumulates roundoff when stored as float32. + vals = np.array([5.0] * 16 + [20.0], dtype=np.float64).reshape(17, 1, 1) + vals[3, 0, 0] = 4.9999997 + vals[7, 0, 0] = 5.0000003 + ds["PCT_NAT_PFT"] = xr.DataArray(vals, dims=("natpft", "lsmlat", "lsmlon")) + + out_path = tmp_path / "surf_nat_patch_precision.nc" + write_surface_nc(ds, str(out_path)) + out = xr.open_dataset(out_path) + + pft_sum_pct = float(out["PCT_NAT_PFT"].sum(dim="natpft").isel(lsmlat=0, lsmlon=0).item()) + wt_nat_patch_sum = pft_sum_pct / 100.0 + + assert abs(wt_nat_patch_sum - 1.0) <= 1e-12 + out.close() + + +def test_write_surface_nc_landunit_weights_near_exact_one(tmp_path): + ds = xr.Dataset( + coords={ + "numurbl": np.arange(2, dtype=np.int32), + "lsmlat": np.arange(1, dtype=np.int32), + "lsmlon": np.arange(1, dtype=np.int32), + } + ) + + ds["PCT_NATVEG"] = xr.DataArray(np.array([[65.432198765]], dtype=np.float64), dims=("lsmlat", "lsmlon")) + ds["PCT_CROP"] = xr.DataArray(np.array([[12.345678901]], dtype=np.float64), dims=("lsmlat", "lsmlon")) + ds["PCT_WETLAND"] = xr.DataArray(np.array([[4.210987654]], dtype=np.float64), dims=("lsmlat", "lsmlon")) + ds["PCT_LAKE"] = xr.DataArray(np.array([[8.765432109]], dtype=np.float64), dims=("lsmlat", "lsmlon")) + ds["PCT_GLACIER"] = xr.DataArray(np.array([[7.777777777]], dtype=np.float64), dims=("lsmlat", "lsmlon")) + ds["PCT_URBAN"] = xr.DataArray( + np.array([[[0.700000001]], [[0.768000002]]], dtype=np.float64), + dims=("numurbl", "lsmlat", "lsmlon"), + ) + + out_path = tmp_path / "surf_lunit_precision.nc" + write_surface_nc(ds, str(out_path)) + out = xr.open_dataset(out_path) + + urban_sum = out["PCT_URBAN"].sum(dim="numurbl") + pct_lunit = ( + out["PCT_NATVEG"] + + out["PCT_CROP"] + + out["PCT_WETLAND"] + + out["PCT_LAKE"] + + out["PCT_GLACIER"] + + urban_sum + ) + wt_lunit_sum = float(pct_lunit.isel(lsmlat=0, lsmlon=0).item()) / 100.0 + + assert abs(wt_lunit_sum - 1.0) <= 1e-12 + out.close() From 5451e57b5aee1b58bc05a47d350607883fcfa7f0 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 14 Jul 2026 16:13:44 -0600 Subject: [PATCH 02/14] Refactor previous commit and extend to topounits This commit refactors the previous commit by introducing a shared module, fraction_closure.py, that handles closure operations and is called from multiple places (surface and landuse variables). It also moves the closure requirements into the schema/specs, rather than having them ad hoc. It also extends the fractional closure calculations to topounits. Tests are added for the topounit closure. --- src/dapper/landuse/landuse.py | 82 +---------- src/dapper/surf/fraction_closure.py | 211 ++++++++++++++++++++++++++++ src/dapper/surf/schema.py | 2 + src/dapper/surf/sfile.py | 159 +-------------------- src/dapper/surf/validate.py | 30 ++++ tests/test_fraction_closure.py | 93 +++++++++++- 6 files changed, 339 insertions(+), 238 deletions(-) create mode 100644 src/dapper/surf/fraction_closure.py diff --git a/src/dapper/landuse/landuse.py b/src/dapper/landuse/landuse.py index 87e6de8..d31ad33 100644 --- a/src/dapper/landuse/landuse.py +++ b/src/dapper/landuse/landuse.py @@ -12,88 +12,10 @@ from dapper.geo import sampling from dapper.domains.domain import Domain +from dapper.surf.fraction_closure import normalize_fraction_closure LonWrap = Literal["auto", "0_360", "-180_180"] - -def _scale_to_target_sum( - da: xr.DataArray, - *, - dim: str, - target: xr.DataArray, - eps: float = 1e-12, -) -> xr.DataArray: - """Scale values along `dim` so their sum matches `target` where possible.""" - work = da.astype(np.float64) - summed = work.sum(dim=dim, skipna=True) - valid = np.isfinite(summed) & (np.abs(summed) > eps) & np.isfinite(target) - safe_den = xr.where(valid, summed, 1.0) - factor = xr.where(valid, target.astype(np.float64) / safe_den, 1.0) - return work * factor - - -def _normalize_landuse_fraction_closure(ds: xr.Dataset) -> xr.Dataset: - """Enforce closure on fraction-like landuse outputs before writing.""" - ds2 = ds.copy(deep=False) - - landunit_scalar_names = ["PCT_NATVEG", "PCT_CROP", "PCT_WETLAND", "PCT_LAKE", "PCT_GLACIER"] - landunit_terms: list[xr.DataArray] = [] - for name in landunit_scalar_names: - if name in ds2: - landunit_terms.append(ds2[name].astype(np.float64)) - - urban_group = None - if "PCT_URBAN" in ds2: - urban = ds2["PCT_URBAN"].astype(np.float64) - if "numurbl" in urban.dims: - urban_group = urban - landunit_terms.append(urban.sum(dim="numurbl", skipna=True)) - else: - landunit_terms.append(urban) - - if len(landunit_terms) >= 2: - current_total = sum(landunit_terms) - near_100 = np.abs(current_total - 100.0) <= 1.0 - valid = np.isfinite(current_total) & (current_total > 1e-12) & near_100 - factor = xr.where(valid, 100.0 / current_total, 1.0) - - for name in landunit_scalar_names: - if name in ds2: - ds2[name] = ds2[name].astype(np.float64) * factor - if urban_group is not None: - ds2["PCT_URBAN"] = urban_group * factor - elif "PCT_URBAN" in ds2: - ds2["PCT_URBAN"] = ds2["PCT_URBAN"].astype(np.float64) * factor - - if "PCT_NAT_PFT" in ds2 and "natpft" in ds2["PCT_NAT_PFT"].dims: - target = xr.full_like( - ds2["PCT_NAT_PFT"].isel(natpft=0, drop=True), - 100.0, - dtype=np.float64, - ) - ds2["PCT_NAT_PFT"] = _scale_to_target_sum(ds2["PCT_NAT_PFT"], dim="natpft", target=target) - - if "PCT_CFT" in ds2 and "cft" in ds2["PCT_CFT"].dims and "PCT_CROP" in ds2: - ds2["PCT_CFT"] = _scale_to_target_sum(ds2["PCT_CFT"], dim="cft", target=ds2["PCT_CROP"].astype(np.float64)) - - if "PCT_GLC_MEC" in ds2 and "nglcec" in ds2["PCT_GLC_MEC"].dims and "PCT_GLACIER" in ds2: - ds2["PCT_GLC_MEC"] = _scale_to_target_sum( - ds2["PCT_GLC_MEC"], - dim="nglcec", - target=ds2["PCT_GLACIER"].astype(np.float64), - ) - - if "FSURF" in ds2 and "FGRD" in ds2: - fs = ds2["FSURF"].astype(np.float64) - fg = ds2["FGRD"].astype(np.float64) - total = fs + fg - valid = np.isfinite(total) & (total > 1e-12) - factor = xr.where(valid, 1.0 / total, 1.0) - ds2["FSURF"] = fs * factor - ds2["FGRD"] = fg * factor - - return ds2 - def sample_landuse_timeseries( src_path: str | Path, df_loc: pd.DataFrame, @@ -349,7 +271,7 @@ def _write_nc(ds_out: xr.Dataset, path: Path) -> None: csv_path = out_path.with_suffix(out_path.suffix + ".zonal_weights.csv") dfw_all.to_csv(csv_path, index=False) - out = _normalize_landuse_fraction_closure(out) + out = normalize_fraction_closure(out) _write_nc(out, out_path) df_summary = pd.DataFrame({gid_col: order, "sample_ncells": ncells, "sample_area_total_m2": area_m2}) diff --git a/src/dapper/surf/fraction_closure.py b/src/dapper/surf/fraction_closure.py new file mode 100644 index 0000000..2511bfc --- /dev/null +++ b/src/dapper/surf/fraction_closure.py @@ -0,0 +1,211 @@ +"""Shared fraction-closure logic for surface and landuse datasets.""" + +from __future__ import annotations + +from typing import Iterable + +import numpy as np +import xarray as xr + +# Keep closures aligned with the canonical surface variable specs. +from dapper.surf.surface_var_specs import SURFACE_VAR_SPECS + + +def _scale_to_target_sum( + da: xr.DataArray, + *, + dim: str, + target: xr.DataArray, + eps: float = 1e-12, +) -> xr.DataArray: + """Scale values along `dim` so their sum matches `target` where possible.""" + work = da.astype(np.float64) + summed = work.sum(dim=dim, skipna=True) + valid = np.isfinite(summed) & (np.abs(summed) > eps) & np.isfinite(target) + safe_den = xr.where(valid, summed, 1.0) + factor = xr.where(valid, target.astype(np.float64) / safe_den, 1.0) + return work * factor + + +def _snap_partition_sum( + da: xr.DataArray, + *, + dim: str, + target: xr.DataArray, + eps: float = 1e-15, +) -> xr.DataArray: + """Force exact closure by assigning residual to the last band in float64.""" + work = da.astype(np.float64) + if dim not in work.dims or int(work.sizes.get(dim, 0)) < 1: + return work + + head = work.isel({dim: slice(0, -1)}) + head_sum = head.sum(dim=dim, skipna=True) + last_src = work.isel({dim: -1}) + candidate_last = target.astype(np.float64) - head_sum + + valid = np.isfinite(target) & np.isfinite(head_sum) + snapped_last = xr.where(valid, candidate_last, last_src) + snapped_last = xr.where(np.abs(snapped_last) < eps, 0.0, snapped_last) + + out = work.copy(deep=False) + out[{dim: -1}] = snapped_last + return out + + +def _close_partition( + ds: xr.Dataset, + *, + var_name: str, + dim: str, + target: xr.DataArray, +) -> None: + if var_name not in ds: + return + da = ds[var_name] + if dim not in da.dims: + return + ds[var_name] = _scale_to_target_sum(da, dim=dim, target=target) + ds[var_name] = _snap_partition_sum(ds[var_name], dim=dim, target=target) + + +def _close_unit_partition(ds: xr.Dataset, *, left: str, right: str) -> None: + if left not in ds or right not in ds: + return + lhs = ds[left].astype(np.float64) + rhs = ds[right].astype(np.float64) + total = lhs + rhs + valid = np.isfinite(total) & (total > 1e-12) + factor = xr.where(valid, 1.0 / total, 1.0) + ds[left] = lhs * factor + ds[right] = rhs * factor + + +def _landunit_scalar_names(ds: xr.Dataset) -> list[str]: + # Explicitly constrained to ELM landunit classes. We gate against canonical + # specs so stale names are ignored automatically. + candidates = ["PCT_NATVEG", "PCT_CROP", "PCT_WETLAND", "PCT_LAKE", "PCT_GLACIER"] + return [v for v in candidates if (v in ds and v in SURFACE_VAR_SPECS)] + + +def _apply_landunit_total_closure(ds: xr.Dataset) -> None: + """Enforce landunit totals to close to 100 (including urban aggregate).""" + scalar_names = _landunit_scalar_names(ds) + landunit_terms: list[xr.DataArray] = [ds[name].astype(np.float64) for name in scalar_names] + + urban_group = None + if "PCT_URBAN" in ds: + urban = ds["PCT_URBAN"].astype(np.float64) + if "numurbl" in urban.dims: + urban_group = urban + landunit_terms.append(urban.sum(dim="numurbl", skipna=True)) + else: + landunit_terms.append(urban) + + if len(landunit_terms) < 2: + return + + current_total = sum(landunit_terms) + near_100 = np.abs(current_total - 100.0) <= 1.0 + valid = np.isfinite(current_total) & (current_total > 1e-12) & near_100 + factor = xr.where(valid, 100.0 / current_total, 1.0) + + for name in scalar_names: + ds[name] = ds[name].astype(np.float64) * factor + if urban_group is not None: + ds["PCT_URBAN"] = urban_group * factor + elif "PCT_URBAN" in ds: + ds["PCT_URBAN"] = ds["PCT_URBAN"].astype(np.float64) * factor + + total_after = 0.0 + for name in scalar_names: + total_after = total_after + ds[name].astype(np.float64) + if "PCT_URBAN" in ds: + urb = ds["PCT_URBAN"].astype(np.float64) + if "numurbl" in urb.dims: + total_after = total_after + urb.sum(dim="numurbl", skipna=True) + else: + total_after = total_after + urb + + resid = 100.0 - total_after + + # Prefer scalar classes first; fallback to last urban class. + snapped = False + for name in ("PCT_GLACIER", "PCT_LAKE", "PCT_WETLAND", "PCT_CROP", "PCT_NATVEG"): + if name in ds: + ds[name] = ds[name].astype(np.float64) + resid + snapped = True + break + + if (not snapped) and ("PCT_URBAN" in ds): + urb = ds["PCT_URBAN"].astype(np.float64) + if "numurbl" in urb.dims and int(urb.sizes.get("numurbl", 0)) >= 1: + urb[{"numurbl": -1}] = urb.isel(numurbl=-1) + resid + ds["PCT_URBAN"] = urb + else: + ds["PCT_URBAN"] = urb + resid + + +def _full_like_from_partition(ds: xr.Dataset, *, var_name: str, dim: str, value: float) -> xr.DataArray | None: + if var_name not in ds or dim not in ds[var_name].dims: + return None + return xr.full_like(ds[var_name].isel({dim: 0}, drop=True), value, dtype=np.float64) + + +def normalize_fraction_closure(ds: xr.Dataset) -> xr.Dataset: + """Apply canonical fraction closure for surface/landuse datasets.""" + ds2 = ds.copy(deep=False) + + _apply_landunit_total_closure(ds2) + + # Natural-patch weights should sum to 100 for every non-natpft index tuple. + tgt_nat = _full_like_from_partition(ds2, var_name="PCT_NAT_PFT", dim="natpft", value=100.0) + if tgt_nat is not None: + _close_partition(ds2, var_name="PCT_NAT_PFT", dim="natpft", target=tgt_nat) + + if "PCT_CFT" in ds2 and "cft" in ds2["PCT_CFT"].dims and "PCT_CROP" in ds2: + _close_partition( + ds2, + var_name="PCT_CFT", + dim="cft", + target=ds2["PCT_CROP"].astype(np.float64), + ) + + if "PCT_GLC_MEC" in ds2 and "nglcec" in ds2["PCT_GLC_MEC"].dims and "PCT_GLACIER" in ds2: + _close_partition( + ds2, + var_name="PCT_GLC_MEC", + dim="nglcec", + target=ds2["PCT_GLACIER"].astype(np.float64), + ) + + # Topounit area weights can appear under either canonical name. + for top_var in ("PCT_TOPUNIT", "TopounitFracArea"): + tgt_top = _full_like_from_partition(ds2, var_name=top_var, dim="topounit", value=100.0) + if tgt_top is not None: + _close_partition(ds2, var_name=top_var, dim="topounit", target=tgt_top) + + _close_unit_partition(ds2, left="FSURF", right="FGRD") + + return ds2 + + +def closure_critical_variables(present_vars: Iterable[str]) -> set[str]: + """Return vars that should be written as float64 to preserve closure.""" + present = set(present_vars) + critical = { + "PCT_NAT_PFT", + "PCT_CFT", + "PCT_GLC_MEC", + "PCT_NATVEG", + "PCT_CROP", + "PCT_GLACIER", + "PCT_WETLAND", + "PCT_LAKE", + "PCT_URBAN", + "PCT_TOPUNIT", + "TopounitFracArea", + "FSURF", + "FGRD", + } + return {v for v in critical if v in present} \ No newline at end of file diff --git a/src/dapper/surf/schema.py b/src/dapper/surf/schema.py index 35dc92d..2fb0060 100644 --- a/src/dapper/surf/schema.py +++ b/src/dapper/surf/schema.py @@ -108,6 +108,7 @@ def pdef( units: str = "", doc: str = "", required_level: str = "", + contexts: Tuple[str, ...] = (), **attrs, ) -> ParDef: """ @@ -128,6 +129,7 @@ def pdef( doc=doc, required_level=required_level, attrs=attrs or {}, + contexts=tuple(contexts or ()), ) diff --git a/src/dapper/surf/sfile.py b/src/dapper/surf/sfile.py index 58363c8..0d1a0f6 100644 --- a/src/dapper/surf/sfile.py +++ b/src/dapper/surf/sfile.py @@ -12,6 +12,7 @@ from dapper.surf import sample as SP # for from_halfdegree_point from dapper.geo import sampling # shared gridded sampler from dapper.surf.surface_var_specs import SURFACE_VAR_SPECS +from dapper.surf.fraction_closure import normalize_fraction_closure, closure_critical_variables ArrayLike = Union[np.ndarray, "xr.DataArray", float, int] @@ -187,7 +188,7 @@ def write_surface_nc( import datetime as _dt # ---- global attrs ---- - ds2 = _normalize_surface_fraction_closure(ds) + ds2 = normalize_fraction_closure(ds) merged = dict(ds2.attrs) if dapper_attrs: @@ -202,10 +203,7 @@ def write_surface_nc( ds2.attrs = merged - closure_critical = { - "PCT_NAT_PFT", "PCT_CFT", "PCT_GLC_MEC", - "PCT_NATVEG", "PCT_CROP", "PCT_GLACIER", "PCT_WETLAND", "PCT_LAKE", "PCT_URBAN", - } + closure_critical = closure_critical_variables(ds2.data_vars) enc: Dict[str, dict] = {} for v in ds2.data_vars: # Fill values must match dtype. Here we do float32 to match ELM surface expectations. @@ -219,157 +217,6 @@ def write_surface_nc( return out_path -def _scale_to_target_sum( - da: xr.DataArray, - *, - dim: str, - target: xr.DataArray, - eps: float = 1e-12, -) -> xr.DataArray: - """Scale values along `dim` so their sum matches `target` where possible.""" - work = da.astype(np.float64) - summed = work.sum(dim=dim, skipna=True) - valid = np.isfinite(summed) & (np.abs(summed) > eps) & np.isfinite(target) - safe_den = xr.where(valid, summed, 1.0) - factor = xr.where(valid, target.astype(np.float64) / safe_den, 1.0) - return work * factor - - -def _snap_partition_sum( - da: xr.DataArray, - *, - dim: str, - target: xr.DataArray, - eps: float = 1e-15, -) -> xr.DataArray: - """Force exact closure by assigning residual to the last band in float64.""" - work = da.astype(np.float64) - if dim not in work.dims or int(work.sizes.get(dim, 0)) < 1: - return work - - head = work.isel({dim: slice(0, -1)}) - head_sum = head.sum(dim=dim, skipna=True) - last_src = work.isel({dim: -1}) - candidate_last = target.astype(np.float64) - head_sum - - valid = np.isfinite(target) & np.isfinite(head_sum) - snapped_last = xr.where(valid, candidate_last, last_src) - snapped_last = xr.where(np.abs(snapped_last) < eps, 0.0, snapped_last) - - out = work.copy(deep=False) - out[{dim: -1}] = snapped_last - return out - - -def _normalize_surface_fraction_closure(ds: xr.Dataset) -> xr.Dataset: - """Enforce closure on key fractional groups before writing a surface file.""" - ds2 = ds.copy(deep=False) - - # Landunit percentages should close to 100 across available classes. - landunit_scalar_names = ["PCT_NATVEG", "PCT_CROP", "PCT_WETLAND", "PCT_LAKE", "PCT_GLACIER"] - landunit_terms: list[xr.DataArray] = [] - for name in landunit_scalar_names: - if name in ds2: - landunit_terms.append(ds2[name].astype(np.float64)) - - urban_group = None - if "PCT_URBAN" in ds2: - urban = ds2["PCT_URBAN"].astype(np.float64) - if "numurbl" in urban.dims: - urban_group = urban - landunit_terms.append(urban.sum(dim="numurbl", skipna=True)) - else: - landunit_terms.append(urban) - - if len(landunit_terms) >= 2: - current_total = sum(landunit_terms) - # Only adjust small drift; avoid force-normalizing partial datasets. - near_100 = np.abs(current_total - 100.0) <= 1.0 - valid = np.isfinite(current_total) & (current_total > 1e-12) & near_100 - factor = xr.where(valid, 100.0 / current_total, 1.0) - - for name in landunit_scalar_names: - if name in ds2: - ds2[name] = ds2[name].astype(np.float64) * factor - if urban_group is not None: - ds2["PCT_URBAN"] = urban_group * factor - elif "PCT_URBAN" in ds2: - ds2["PCT_URBAN"] = ds2["PCT_URBAN"].astype(np.float64) * factor - - # Force exact landunit closure by assigning residual to one component. - total_after = 0.0 - for name in landunit_scalar_names: - if name in ds2: - total_after = total_after + ds2[name].astype(np.float64) - if "PCT_URBAN" in ds2: - urb = ds2["PCT_URBAN"].astype(np.float64) - if "numurbl" in urb.dims: - total_after = total_after + urb.sum(dim="numurbl", skipna=True) - else: - total_after = total_after + urb - - resid = 100.0 - total_after - - # Prefer scalar classes first; fallback to last urban class. - snapped = False - for name in ("PCT_GLACIER", "PCT_LAKE", "PCT_WETLAND", "PCT_CROP", "PCT_NATVEG"): - if name in ds2: - ds2[name] = ds2[name].astype(np.float64) + resid - snapped = True - break - - if (not snapped) and ("PCT_URBAN" in ds2): - urb = ds2["PCT_URBAN"].astype(np.float64) - if "numurbl" in urb.dims and int(urb.sizes.get("numurbl", 0)) >= 1: - urb[{"numurbl": -1}] = urb.isel(numurbl=-1) + resid - ds2["PCT_URBAN"] = urb - else: - ds2["PCT_URBAN"] = urb + resid - - # PFT partitions are natural-patch weights and must close to 100. - if "PCT_NAT_PFT" in ds2 and "natpft" in ds2["PCT_NAT_PFT"].dims: - pft = ds2["PCT_NAT_PFT"] - target = xr.full_like( - pft.isel(natpft=0, drop=True), - 100.0, - dtype=np.float64, - ) - ds2["PCT_NAT_PFT"] = _scale_to_target_sum(pft, dim="natpft", target=target) - ds2["PCT_NAT_PFT"] = _snap_partition_sum(ds2["PCT_NAT_PFT"], dim="natpft", target=target) - - # Crop and glacier-class partitions should close to their parent percentages. - if "PCT_CFT" in ds2 and "cft" in ds2["PCT_CFT"].dims and "PCT_CROP" in ds2: - ds2["PCT_CFT"] = _scale_to_target_sum(ds2["PCT_CFT"], dim="cft", target=ds2["PCT_CROP"].astype(np.float64)) - ds2["PCT_CFT"] = _snap_partition_sum( - ds2["PCT_CFT"], - dim="cft", - target=ds2["PCT_CROP"].astype(np.float64), - ) - if "PCT_GLC_MEC" in ds2 and "nglcec" in ds2["PCT_GLC_MEC"].dims and "PCT_GLACIER" in ds2: - ds2["PCT_GLC_MEC"] = _scale_to_target_sum( - ds2["PCT_GLC_MEC"], - dim="nglcec", - target=ds2["PCT_GLACIER"].astype(np.float64), - ) - ds2["PCT_GLC_MEC"] = _snap_partition_sum( - ds2["PCT_GLC_MEC"], - dim="nglcec", - target=ds2["PCT_GLACIER"].astype(np.float64), - ) - - # Unitless partition currently represented by irrigation split. - if "FSURF" in ds2 and "FGRD" in ds2: - fs = ds2["FSURF"].astype(np.float64) - fg = ds2["FGRD"].astype(np.float64) - total = fs + fg - valid = np.isfinite(total) & (total > 1e-12) - factor = xr.where(valid, 1.0 / total, 1.0) - ds2["FSURF"] = fs * factor - ds2["FGRD"] = fg * factor - - return ds2 - - class CustomizeError(ValueError): """Raised when a customization fails schema/formatting validation.""" diff --git a/src/dapper/surf/validate.py b/src/dapper/surf/validate.py index 6a36350..6ebe516 100644 --- a/src/dapper/surf/validate.py +++ b/src/dapper/surf/validate.py @@ -329,6 +329,36 @@ def _check_soft_consistency(self, ds: xr.Dataset) -> List[CheckResult]: r.append(CheckResult("V-105.consistency.pftsum", "WARN", max_diff <= 1e-6, f"mean(|sum(PCT_NAT_PFT)-100|)={mean_diff:.3e}; max={max_diff:.3e}")) + # landunit closure ≈ 100 across available classes (incl urban aggregate) + landunit_terms: List[xr.DataArray] = [] + for name in ("PCT_NATVEG", "PCT_CROP", "PCT_WETLAND", "PCT_LAKE", "PCT_GLACIER"): + if name in ds: + landunit_terms.append(ds[name].astype("float64")) + if "PCT_URBAN" in ds: + urb = ds["PCT_URBAN"].astype("float64") + if "numurbl" in urb.dims: + landunit_terms.append(urb.sum(dim="numurbl", skipna=True)) + else: + landunit_terms.append(urb) + if len(landunit_terms) >= 2: + lsum = sum(landunit_terms) + lresid = np.abs(np.asarray(lsum.values, dtype="float64") - 100.0) + mean_diff = float(np.nanmean(lresid)) + max_diff = float(np.nanmax(lresid)) + r.append(CheckResult("V-108.consistency.landunitsum", "WARN", max_diff <= 1e-6, + f"mean(|landunit_sum-100|)={mean_diff:.3e}; max={max_diff:.3e}")) + + # topounit area weights should close to 100 across topounit dim. + for top_var in ("PCT_TOPUNIT", "TopounitFracArea"): + if top_var in ds and "topounit" in ds[top_var].dims: + tsum = ds[top_var].sum(dim="topounit", skipna=True) + tresid = np.abs(np.asarray(tsum.values, dtype="float64") - 100.0) + mean_diff = float(np.nanmean(tresid)) + max_diff = float(np.nanmax(tresid)) + r.append(CheckResult("V-109.consistency.topounitsum", "WARN", max_diff <= 1e-6, + f"{top_var}: mean(|sum-100|)={mean_diff:.3e}; max={max_diff:.3e}")) + break + # If any cell has PCT_URBAN>0 → URBAN_REGION_ID present if "PCT_URBAN" in ds: urb_max = np.nanmax(np.asarray(ds["PCT_URBAN"].values)) diff --git a/tests/test_fraction_closure.py b/tests/test_fraction_closure.py index 45e16e9..8b73f09 100644 --- a/tests/test_fraction_closure.py +++ b/tests/test_fraction_closure.py @@ -2,7 +2,7 @@ import xarray as xr from dapper.surf.sfile import write_surface_nc -from dapper.landuse.landuse import _normalize_landuse_fraction_closure +from dapper.surf.fraction_closure import normalize_fraction_closure def test_write_surface_nc_enforces_fraction_closure(tmp_path): @@ -103,7 +103,7 @@ def test_landuse_normalizer_enforces_fraction_closure(): dims=("natpft", "lsmlat", "lsmlon"), ) - fixed = _normalize_landuse_fraction_closure(ds) + fixed = normalize_fraction_closure(ds) pft_sum = fixed["PCT_NAT_PFT"].sum(dim="natpft") np.testing.assert_allclose(pft_sum.values, 100.0, atol=1e-10, rtol=0.0) @@ -170,3 +170,92 @@ def test_write_surface_nc_landunit_weights_near_exact_one(tmp_path): assert abs(wt_lunit_sum - 1.0) <= 1e-12 out.close() + + +def test_write_surface_nc_topounit_pft_closure_exact(tmp_path): + ds = xr.Dataset( + coords={ + "topounit": np.arange(2, dtype=np.int32), + "natpft": np.arange(4, dtype=np.int32), + "lsmlat": np.arange(1, dtype=np.int32), + "lsmlon": np.arange(1, dtype=np.int32), + } + ) + + # Two topounits with tiny partition drift in opposite directions. + vals = np.array( + [ + [[[25.0]], [[25.0]]], + [[[25.0]], [[25.0]]], + [[[25.0]], [[25.0]]], + [[[24.9999998]], [[25.0000002]]], + ], + dtype=np.float64, + ) + ds["PCT_NAT_PFT"] = xr.DataArray(vals, dims=("natpft", "topounit", "lsmlat", "lsmlon")) + + out_path = tmp_path / "surf_topounit_natpft_precision.nc" + write_surface_nc(ds, str(out_path)) + + out = xr.open_dataset(out_path) + pft_sum = out["PCT_NAT_PFT"].sum(dim="natpft") + np.testing.assert_allclose(pft_sum.values, 100.0, atol=1e-12, rtol=0.0) + out.close() + + +def test_topounit_landunit_closure_exact(): + ds = xr.Dataset( + coords={ + "topounit": np.arange(2, dtype=np.int32), + "numurbl": np.arange(2, dtype=np.int32), + "lsmlat": np.arange(1, dtype=np.int32), + "lsmlon": np.arange(1, dtype=np.int32), + } + ) + + ds["PCT_NATVEG"] = xr.DataArray(np.array([[[50.0]], [[40.0]]]), dims=("topounit", "lsmlat", "lsmlon")) + ds["PCT_CROP"] = xr.DataArray(np.array([[[20.0]], [[30.0]]]), dims=("topounit", "lsmlat", "lsmlon")) + ds["PCT_WETLAND"] = xr.DataArray(np.array([[[10.0]], [[10.0]]]), dims=("topounit", "lsmlat", "lsmlon")) + ds["PCT_LAKE"] = xr.DataArray(np.array([[[9.0]], [[8.0]]]), dims=("topounit", "lsmlat", "lsmlon")) + ds["PCT_GLACIER"] = xr.DataArray(np.array([[[9.0]], [[10.0]]]), dims=("topounit", "lsmlat", "lsmlon")) + ds["PCT_URBAN"] = xr.DataArray( + np.array( + [ + [[[1.0]], [[1.1]]], + [[[1.0]], [[0.9]]], + ], + dtype=np.float64, + ), + dims=("numurbl", "topounit", "lsmlat", "lsmlon"), + ) + + fixed = normalize_fraction_closure(ds) + urban_sum = fixed["PCT_URBAN"].sum(dim="numurbl") + total = ( + fixed["PCT_NATVEG"] + + fixed["PCT_CROP"] + + fixed["PCT_WETLAND"] + + fixed["PCT_LAKE"] + + fixed["PCT_GLACIER"] + + urban_sum + ) + np.testing.assert_allclose(total.values, 100.0, atol=1e-12, rtol=0.0) + + +def test_topounit_weight_aliases_close_to_100(): + ds = xr.Dataset( + coords={ + "topounit": np.arange(3, dtype=np.int32), + "lsmlat": np.arange(1, dtype=np.int32), + "lsmlon": np.arange(1, dtype=np.int32), + } + ) + + ds["PCT_TOPUNIT"] = xr.DataArray( + np.array([[[40.0]], [[30.0]], [[29.999999]]], dtype=np.float64), + dims=("topounit", "lsmlat", "lsmlon"), + ) + + fixed = normalize_fraction_closure(ds) + tot = fixed["PCT_TOPUNIT"].sum(dim="topounit") + np.testing.assert_allclose(tot.values, 100.0, atol=1e-12, rtol=0.0) From 8aea1acaac07ad2faa1d7a3ca9c141f2dac9758d Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 15 Jul 2026 11:18:33 -0600 Subject: [PATCH 03/14] correct a comment --- src/dapper/surf/sfile.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/dapper/surf/sfile.py b/src/dapper/surf/sfile.py index 0d1a0f6..aaf1498 100644 --- a/src/dapper/surf/sfile.py +++ b/src/dapper/surf/sfile.py @@ -206,7 +206,8 @@ def write_surface_nc( closure_critical = closure_critical_variables(ds2.data_vars) enc: Dict[str, dict] = {} for v in ds2.data_vars: - # Fill values must match dtype. Here we do float32 to match ELM surface expectations. + # Fill values must match dtype. Most float vars are written as float32, + # while closure-critical fractions stay float64 to preserve exact sums. if ds2[v].dtype.kind == "f": if v in closure_critical: enc[v] = {"dtype": "float64", "_FillValue": np.float64(-9.96921e36)} From 120c54dc7e7bebc4091a155b481d52a7e81a0ebc Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 16 Jul 2026 13:06:26 -0600 Subject: [PATCH 04/14] fix: use TopounitFracArea instead of PCT_TOPUNIT and normalize to decimal fraction - Update fraction_closure to handle only TopounitFracArea with decimal target (1.0) - Remove PCT_TOPUNIT from closure_critical_variables set - Update validate.py to check TopounitFracArea closure (target=1.0, not 100) - Consolidate separate validation checks for TopounitFracArea - Update test to verify TopounitFracArea closes to 1.0 (not 100) This fixes the unit mismatch where ELM expects TopounitFracArea as a decimal fraction (0-1) but dapper was producing percent units (0-100). --- src/dapper/surf/fraction_closure.py | 10 ++++------ src/dapper/surf/validate.py | 18 ++++++++---------- 2 files changed, 12 insertions(+), 16 deletions(-) diff --git a/src/dapper/surf/fraction_closure.py b/src/dapper/surf/fraction_closure.py index 2511bfc..afcaf18 100644 --- a/src/dapper/surf/fraction_closure.py +++ b/src/dapper/surf/fraction_closure.py @@ -179,11 +179,10 @@ def normalize_fraction_closure(ds: xr.Dataset) -> xr.Dataset: target=ds2["PCT_GLACIER"].astype(np.float64), ) - # Topounit area weights can appear under either canonical name. - for top_var in ("PCT_TOPUNIT", "TopounitFracArea"): - tgt_top = _full_like_from_partition(ds2, var_name=top_var, dim="topounit", value=100.0) - if tgt_top is not None: - _close_partition(ds2, var_name=top_var, dim="topounit", target=tgt_top) + # TopounitFracArea is a decimal fraction (0-1) and should sum to 1.0 over topounit. + tgt_top = _full_like_from_partition(ds2, var_name="TopounitFracArea", dim="topounit", value=1.0) + if tgt_top is not None: + _close_partition(ds2, var_name="TopounitFracArea", dim="topounit", target=tgt_top) _close_unit_partition(ds2, left="FSURF", right="FGRD") @@ -203,7 +202,6 @@ def closure_critical_variables(present_vars: Iterable[str]) -> set[str]: "PCT_WETLAND", "PCT_LAKE", "PCT_URBAN", - "PCT_TOPUNIT", "TopounitFracArea", "FSURF", "FGRD", diff --git a/src/dapper/surf/validate.py b/src/dapper/surf/validate.py index 6ebe516..05cfdd2 100644 --- a/src/dapper/surf/validate.py +++ b/src/dapper/surf/validate.py @@ -348,16 +348,14 @@ def _check_soft_consistency(self, ds: xr.Dataset) -> List[CheckResult]: r.append(CheckResult("V-108.consistency.landunitsum", "WARN", max_diff <= 1e-6, f"mean(|landunit_sum-100|)={mean_diff:.3e}; max={max_diff:.3e}")) - # topounit area weights should close to 100 across topounit dim. - for top_var in ("PCT_TOPUNIT", "TopounitFracArea"): - if top_var in ds and "topounit" in ds[top_var].dims: - tsum = ds[top_var].sum(dim="topounit", skipna=True) - tresid = np.abs(np.asarray(tsum.values, dtype="float64") - 100.0) - mean_diff = float(np.nanmean(tresid)) - max_diff = float(np.nanmax(tresid)) - r.append(CheckResult("V-109.consistency.topounitsum", "WARN", max_diff <= 1e-6, - f"{top_var}: mean(|sum-100|)={mean_diff:.3e}; max={max_diff:.3e}")) - break + # TopounitFracArea is a decimal fraction; sum over topounit should close to 1. + if "TopounitFracArea" in ds and "topounit" in ds["TopounitFracArea"].dims: + tsum = ds["TopounitFracArea"].sum(dim="topounit", skipna=True) + tresid = np.abs(np.asarray(tsum.values, dtype="float64") - 1.0) + mean_diff = float(np.nanmean(tresid)) + max_diff = float(np.nanmax(tresid)) + r.append(CheckResult("V-109.consistency.topounitfracsum", "WARN", max_diff <= 1e-6, + f"TopounitFracArea: mean(|sum-1|)={mean_diff:.3e}; max={max_diff:.3e}")) # If any cell has PCT_URBAN>0 → URBAN_REGION_ID present if "PCT_URBAN" in ds: From da0332226fcd7ee748ed4672d0909d8b974b972f Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 16 Jul 2026 13:06:35 -0600 Subject: [PATCH 05/14] fix: expand topounit-indexed variables to include topounit dimension After add_topounits_from_domain(), all SURFACE_VAR_SPECS variables with 'topounit' in their dims are now expanded to include the topounit dimension if they lack it. This uses uniform broadcast (all topounits inherit the parent cell's distribution). Fixes the issue where PCT_NAT_PFT and other topounit-dependent variables were written without a topounit dimension, causing ELM to fail with dimension mismatch errors. Also updates add_topounits_from_domain() to use TopounitFracArea instead of PCT_TOPUNIT. Includes comprehensive test to verify: - PCT_NAT_PFT gains topounit dimension with correct ordering - PFT fractions still sum to 100 per topounit - all topounit-indexed vars are properly expanded --- src/dapper/surf/sfile.py | 34 ++++++++++++++++--- tests/test_fraction_closure.py | 61 +++++++++++++++++++++++++++++++--- 2 files changed, 85 insertions(+), 10 deletions(-) diff --git a/src/dapper/surf/sfile.py b/src/dapper/surf/sfile.py index aaf1498..00d7b02 100644 --- a/src/dapper/surf/sfile.py +++ b/src/dapper/surf/sfile.py @@ -892,7 +892,7 @@ def export( agg_policy=agg_policy, ) - # Attach topounit parameters (and PCT_TOPUNIT) exactly once + # Attach topounit parameters (and TopounitFracArea) exactly once if attach_topounits and getattr(run_dom, "topounits", None) is not None and run_dom.topounits is not None: sf.add_topounits_from_domain(run_dom) @@ -1027,7 +1027,7 @@ def add_topounits_from_domain( id_col: str = "topounit_id", pct_col: str = "TopounitPctOfCell", dim_name: str = "topounit", - pct_var_name: str = "PCT_TOPUNIT", + pct_var_name: str = "TopounitFracArea", ) -> None: """ Attach topounits + per-cell weights to the surface dataset. @@ -1097,8 +1097,8 @@ def add_topounits_from_domain( s = float(np.nansum(vals)) if not np.isfinite(s) or s <= 0: raise ValueError(f"Topounit pct weights for gid={gid} are invalid (sum={s}).") - # normalize to 100 just in case - vals = 100.0 * (vals / s) + # normalize to 1.0 (decimal fraction) just in case + vals = 1.0 * (vals / s) for tid, v in zip(grp[id_col].astype(str).tolist(), vals): k = id_to_k[tid] @@ -1114,7 +1114,31 @@ def add_topounits_from_domain( raise ValueError(f"Existing {dim_name} coord does not match topounit ids from domain.") ds[pct_var_name] = xr.DataArray(pct, dims=(dim_name, lat_dim, lon_dim)) - ds[pct_var_name].attrs.update({"long_name": "percent of gridcell in each topounit", "units": "percent"}) + ds[pct_var_name].attrs.update({"long_name": "fraction of gridcell area in each topounit", "units": "unitless"}) + + # Expand topounit-indexed variables that exist in ds but currently lack the + # topounit dimension. All topounits in a grid cell inherit the parent cell's + # distribution uniformly (per-topounit differentiation is a future extension). + top_coord = ds.coords[dim_name] + n_top = len(top_ids) + for _var_name, _spec in SURFACE_VAR_SPECS.items(): + if _var_name not in ds: + continue + _spec_dims = [d.strip() for d in _spec.get("dims", "").split(",")] + if dim_name not in _spec_dims: + continue + _da = ds[_var_name] + if dim_name in _da.dims: + continue # already has the topounit dim + # Repeat identical values across all topounits + _expanded = xr.concat([_da] * n_top, + dim=xr.DataArray(top_coord.values, dims=[dim_name], name=dim_name)) + _expanded[dim_name] = top_coord + # Reorder dims to match spec order (topounit first, spatial last) + _existing = set(_expanded.dims) + _ordered = [d for d in _spec_dims if d in _existing] + _extra = [d for d in _expanded.dims if d not in _ordered] + ds[_var_name] = _expanded.transpose(*_ordered, *_extra) self.ds = ds diff --git a/tests/test_fraction_closure.py b/tests/test_fraction_closure.py index 8b73f09..2a05c97 100644 --- a/tests/test_fraction_closure.py +++ b/tests/test_fraction_closure.py @@ -242,7 +242,8 @@ def test_topounit_landunit_closure_exact(): np.testing.assert_allclose(total.values, 100.0, atol=1e-12, rtol=0.0) -def test_topounit_weight_aliases_close_to_100(): +def test_topounit_fracarea_closure_exact(): + """TopounitFracArea should close to 1.0 (not 100) across topounit dimension.""" ds = xr.Dataset( coords={ "topounit": np.arange(3, dtype=np.int32), @@ -251,11 +252,61 @@ def test_topounit_weight_aliases_close_to_100(): } ) - ds["PCT_TOPUNIT"] = xr.DataArray( - np.array([[[40.0]], [[30.0]], [[29.999999]]], dtype=np.float64), + ds["TopounitFracArea"] = xr.DataArray( + np.array([[[0.40]], [[0.30]], [[0.29999999]]], dtype=np.float64), dims=("topounit", "lsmlat", "lsmlon"), ) fixed = normalize_fraction_closure(ds) - tot = fixed["PCT_TOPUNIT"].sum(dim="topounit") - np.testing.assert_allclose(tot.values, 100.0, atol=1e-12, rtol=0.0) + tot = fixed["TopounitFracArea"].sum(dim="topounit") + np.testing.assert_allclose(tot.values, 1.0, atol=1e-12, rtol=0.0) + + + + + +def test_add_topounits_expands_pct_nat_pft(): + """After add_topounits_from_domain, PCT_NAT_PFT must gain the topounit dimension.""" + import pandas as pd + from unittest.mock import MagicMock + from dapper.surf.sfile import SurfaceFile + + # 2 cells (lsmlat=2), 4 PFTs, 1 longitude — uniform 25% per PFT per cell + natpft_vals = np.ones((4, 2, 1), dtype=np.float64) * 25.0 + ds = xr.Dataset( + { + "PCT_NAT_PFT": xr.DataArray(natpft_vals, dims=("natpft", "lsmlat", "lsmlon")), + "PCT_NATVEG": xr.DataArray(np.ones((2, 1)) * 100.0, dims=("lsmlat", "lsmlon")), + }, + coords={ + "natpft": np.arange(4, dtype=np.int32), + "lsmlat": np.arange(2, dtype=np.int32), + "lsmlon": np.arange(1, dtype=np.int32), + }, + ) + + sf = SurfaceFile(ds) + + # Build a minimal mock domain with 2 cells and 2 topounits each + topounits_df = pd.DataFrame({ + "gid": ["0", "0", "1", "1"], + "topounit_id": ["tu0", "tu1", "tu2", "tu3"], + "TopounitPctOfCell": [60.0, 40.0, 50.0, 50.0], + }) + df_loc = pd.DataFrame({"gid": ["0", "1"]}) + + domain = MagicMock() + domain.topounits = topounits_df + domain.to_df_loc.return_value = df_loc + + sf.add_topounits_from_domain(domain) + + out = sf.ds + assert "topounit" in out["PCT_NAT_PFT"].dims, "PCT_NAT_PFT must have topounit dim after add_topounits_from_domain" + # Dim order: topounit before natpft, spatial last + dims = list(out["PCT_NAT_PFT"].dims) + assert dims.index("topounit") < dims.index("natpft"), "topounit must precede natpft" + assert dims[-2:] == ["lsmlat", "lsmlon"], "spatial dims must be last" + # Each topounit still has PFT fractions summing to 100 + pft_sum = out["PCT_NAT_PFT"].sum(dim="natpft") + np.testing.assert_allclose(pft_sum.values, 100.0, atol=1e-10, rtol=0.0) From 06666b14f7283b2ee27d76c7e8e0e1799abac099 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 16 Jul 2026 14:51:26 -0600 Subject: [PATCH 06/14] docs: update add_topounits_from_domain docstring and regenerate surface variables table --- src/dapper/surf/sfile.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/dapper/surf/sfile.py b/src/dapper/surf/sfile.py index 00d7b02..6f96e42 100644 --- a/src/dapper/surf/sfile.py +++ b/src/dapper/surf/sfile.py @@ -1032,6 +1032,14 @@ def add_topounits_from_domain( """ Attach topounits + per-cell weights to the surface dataset. + This method: + - Adds a topounit dimension and coordinate to the dataset + - Creates a topounit fraction variable (default: TopounitFracArea) with values + normalized to decimal fractions (0.0-1.0, summing to 1.0 per cell) + - Expands topounit-indexed variables in SURFACE_VAR_SPECS (e.g., PCT_NAT_PFT) + to include the topounit dimension via uniform broadcast, where all topounits + in a cell inherit the parent cell's distribution + Expects domain.topounits to exist and contain: - gid_col (links topounit -> cell gid) - id_col (unique id per topounit across the whole run) From 641aa7ea93d88620b257f8ba8388767f9414b77d Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 20 Jul 2026 11:28:44 -0600 Subject: [PATCH 07/14] Fix surface variable dimensions for topounit framework - Add topounit dimension to phosphorus pools (APATITE_P, LABILE_P, OCCLUDED_P, SECONDARY_P) with corrected units (gP/m2) from ncdump metadata - Add topounit dimension to soil properties (ORGANIC, PCT_SAND, PCT_CLAY, PCT_GRVL) with pattern: nlevsoi,topounit,lsmlat,lsmlon - Reorder class dimensions before topounit for consistency: PCT_URBAN, PCT_CFT, NFERT, PFERT, PCT_NAT_PFT, PCT_GLC_MEC, TOPO_GLC_MEC (class dimension now precedes topounit, not vice versa) - Add lsmpft and topounit dimensions to monthly vegetation variables (MONTHLY_LAI, MONTHLY_SAI, MONTHLY_HEIGHT_TOP, MONTHLY_HEIGHT_BOT) with dimension order: time,lsmpft,topounit,lsmlat,lsmlon - Update MaxTopounitElv and TOPO2 documentation to clarify topounit relationships - Add topounits context tag to all affected variables Ensures variables align with actual ELM surface file structure (topounit-based discretization). --- src/dapper/surf/surface_var_specs.py | 100 +++++++++++++-------------- 1 file changed, 50 insertions(+), 50 deletions(-) diff --git a/src/dapper/surf/surface_var_specs.py b/src/dapper/surf/surface_var_specs.py index cacf037..44567e3 100644 --- a/src/dapper/surf/surface_var_specs.py +++ b/src/dapper/surf/surface_var_specs.py @@ -217,7 +217,7 @@ "contexts": ["land_cover", "glaciers", "topounits"], }, "PCT_URBAN": { - "dims": "topounit,numurbl,lsmlat,lsmlon", + "dims": "numurbl,topounit,lsmlat,lsmlon", "doc": "Fraction of each grid cell that is urban, for each " "density class in the multi-density urban scheme.", "required_level": "conditional", @@ -239,7 +239,7 @@ "contexts": ["urban", "topounits"], }, "PCT_GLC_MEC": { - "dims": "topounit,nglcec,lsmlat,lsmlon", + "dims": "nglcec,topounit,lsmlat,lsmlon", "doc": "Percent of grid cell area assigned to mechanistic " "glacier classes (accumulation/ablation, etc.).", "required_level": "conditional", @@ -251,7 +251,7 @@ "contexts": ["land_cover", "glaciers", "topounits"], }, "TOPO_GLC_MEC": { - "dims": "topounit,nglcec,lsmlat,lsmlon", + "dims": "nglcec,topounit,lsmlat,lsmlon", "doc": "Elevation (m) of each mechanistic glacier class.", "required_level": "conditional", "attrs": { @@ -312,7 +312,7 @@ "contexts": ["land_cover", "polygonal_tundra", "topounits"], }, "PCT_CFT": { - "dims": "topounit,cft,lsmlat,lsmlon", + "dims": "cft,topounit,lsmlat,lsmlon", "doc": "Fraction of vegetated area allocated to each crop " "functional type; code aborts if missing when cft " "dimension exists.", @@ -325,7 +325,7 @@ "contexts": ["land_cover", "crops_irrigation", "topounits"], }, "NFERT": { - "dims": "topounit,cft,lsmlat,lsmlon", + "dims": "cft,topounit,lsmlat,lsmlon", "doc": "Nitrogen fertilizer application for each crop functional " "type; if absent, values default to zero.", "required_level": "optional", @@ -334,7 +334,7 @@ "contexts": ["crops_irrigation", "topounits"], }, "PFERT": { - "dims": "topounit,cft,lsmlat,lsmlon", + "dims": "cft,topounit,lsmlat,lsmlon", "doc": "Phosphorus fertilizer application for each crop functional " "type; treated like NFERT.", "required_level": "optional", @@ -342,7 +342,7 @@ "contexts": ["crops_irrigation", "topounits", "phosphorus_cycle"], }, "PCT_NAT_PFT": { - "dims": "topounit,natpft,lsmlat,lsmlon", + "dims": "natpft,topounit,lsmlat,lsmlon", "doc": "Fraction of vegetated area allocated to each natural " "plant functional type; code aborts if missing.", "required_level": "required", @@ -385,8 +385,8 @@ }, "MaxTopounitElv": { "dims": "lsmlat,lsmlon", - "doc": "Maximum elevation (m) among topounits for each " - "grid cell; read only if present.", + "doc": "Maximum topounits elevation in each grid cell; a summary " + "statistic reflecting the highest elevation among topounits.", "required_level": "optional", "attrs": { "requirement": "Optional (not required, but " @@ -443,10 +443,10 @@ }, "TOPO2": { "dims": "lsmlat,lsmlon", - "doc": "Second topography field used in the ELM topounit framework; " - "read only if present.", + "doc": "Weighted average of topounits elevation in each grid cell; " + "a summary statistic for topounit-based elevation characterization.", "required_level": "optional", - "units": "unknown", + "units": "m", "contexts": ["grid_topography", "topounits"], }, "AREA": { @@ -476,36 +476,36 @@ "contexts": ["grid_topography"], }, "PCT_SAND": { - "dims": "nlevsoi,lsmlat,lsmlon", + "dims": "nlevsoi,topounit,lsmlat,lsmlon", "doc": "Soil sand percentage by mass (0–100) in each soil layer; " "controls hydraulic and thermal properties.", "required_level": "required", "units": "percent", - "contexts": ["soil_properties"], + "contexts": ["soil_properties", "topounits"], }, "PCT_CLAY": { - "dims": "nlevsoi,lsmlat,lsmlon", + "dims": "nlevsoi,topounit,lsmlat,lsmlon", "doc": "Soil clay percentage by mass (0–100) in each soil layer; " "controls hydraulic and thermal properties.", "required_level": "required", "units": "percent", - "contexts": ["soil_properties"], + "contexts": ["soil_properties", "topounits"], }, "ORGANIC": { - "dims": "nlevsoi,lsmlat,lsmlon", - "doc": "Soil organic matter or organic carbon per layer; used in " + "dims": "nlevsoi,topounit,lsmlat,lsmlon", + "doc": "Organic matter density per soil layer; used in " "biogeochemical and thermal calculations.", "required_level": "optional", - "units": "unknown", - "contexts": ["soil_properties"], + "units": "kg/m3 (assumed carbon content 0.58 gC per gOM)", + "contexts": ["soil_properties", "topounits"], }, "PCT_GRVL": { - "dims": "nlevsoi,lsmlat,lsmlon", + "dims": "nlevsoi,topounit,lsmlat,lsmlon", "doc": "Percent gravel content (0–100) in each soil layer; " "affects soil water storage and hydraulic conductivity.", "required_level": "optional", "units": "percent", - "contexts": ["soil_properties"], + "contexts": ["soil_properties", "topounits"], }, "GLC_MEC": { "dims": "lsmlat,lsmlon", @@ -520,67 +520,67 @@ "contexts": ["glaciers"], }, "MONTHLY_LAI": { - "dims": "time,lsmlat,lsmlon", + "dims": "time,lsmpft,topounit,lsmlat,lsmlon", "doc": "Leaf area index (LAI; m2 leaf per m2 ground) monthly " "climatology; time is typically 12 months.", "required_level": "optional", "units": "unitless", - "contexts": ["vegetation_structure"], + "contexts": ["vegetation_structure", "topounits"], }, "MONTHLY_SAI": { - "dims": "time,lsmlat,lsmlon", + "dims": "time,lsmpft,topounit,lsmlat,lsmlon", "doc": "Stem area index (SAI) monthly climatology; time is " "typically 12 months.", "required_level": "optional", "units": "unitless", - "contexts": ["vegetation_structure"], + "contexts": ["vegetation_structure", "topounits"], }, "MONTHLY_HEIGHT_TOP": { - "dims": "time,lsmlat,lsmlon", + "dims": "time,lsmpft,topounit,lsmlat,lsmlon", "doc": "Monthly climatology of canopy top height (m) for vegetated landunits.", "required_level": "optional", "units": "m", - "contexts": ["vegetation_structure"], + "contexts": ["vegetation_structure", "topounits"], }, "MONTHLY_HEIGHT_BOT": { - "dims": "time,lsmlat,lsmlon", + "dims": "time,lsmpft,topounit,lsmlat,lsmlon", "doc": "Monthly climatology of canopy bottom height " "(m) for vegetated landunits.", "required_level": "optional", "units": "m", - "contexts": ["vegetation_structure"], + "contexts": ["vegetation_structure", "topounits"], }, "APATITE_P": { - "dims": "lsmlat,lsmlon", - "doc": "Soil phosphorus pool in apatite (primary mineral) form; " - "used by phosphorus biogeochemistry when enabled.", + "dims": "topounit,lsmlat,lsmlon", + "doc": "Apatite phosphorus; soil phosphorus pool in apatite " + "(primary mineral) form; used by phosphorus biogeochemistry when enabled.", "required_level": "optional", - "units": "unknown", - "contexts": ["phosphorus_cycle"], + "units": "gP/m2", + "contexts": ["phosphorus_cycle", "topounits"], }, "LABILE_P": { - "dims": "lsmlat,lsmlon", - "doc": "Soil labile (readily available) phosphorus pool; used by " - "P-cycle parameterizations.", + "dims": "topounit,lsmlat,lsmlon", + "doc": "Labile inorganic phosphorus; soil labile (readily available) " + "phosphorus pool; used by P-cycle parameterizations.", "required_level": "optional", - "units": "unknown", - "contexts": ["phosphorus_cycle"], + "units": "gP/m2", + "contexts": ["phosphorus_cycle", "topounits"], }, "OCCLUDED_P": { - "dims": "lsmlat,lsmlon", - "doc": "Soil occluded phosphorus pool (sorbed or otherwise " - "inaccessible); part of multi-pool P parameterization.", + "dims": "topounit,lsmlat,lsmlon", + "doc": "Occluded phosphorus; soil occluded phosphorus pool (sorbed or " + "otherwise inaccessible); part of multi-pool P parameterization.", "required_level": "optional", - "units": "unknown", - "contexts": ["phosphorus_cycle"], + "units": "gP/m2", + "contexts": ["phosphorus_cycle", "topounits"], }, "SECONDARY_P": { - "dims": "lsmlat,lsmlon", - "doc": "Soil secondary mineral phosphorus pool; intermediate " - "in reactivity between apatite and labile pools.", + "dims": "topounit,lsmlat,lsmlon", + "doc": "Secondary mineral phosphorus; soil secondary mineral phosphorus " + "pool; intermediate in reactivity between apatite and labile pools.", "required_level": "optional", - "units": "unknown", - "contexts": ["phosphorus_cycle"], + "units": "gP/m2", + "contexts": ["phosphorus_cycle", "topounits"], }, } From 800bdc00bb1054e5d49b2f347a3d9bafc437e691 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 20 Jul 2026 11:29:41 -0600 Subject: [PATCH 08/14] Add urban albedo variables for multi-spectral radiation Add 8 new urban albedo variables (diffuse and direct components) for impervious/pervious roads, roofs, and walls across spectral bands. These variables control solar radiation absorption in urban land model. Variables added: - ALB_IMPROAD_DIF, ALB_IMPROAD_DIR - ALB_PERROAD_DIF, ALB_PERROAD_DIR - ALB_ROOF_DIF, ALB_ROOF_DIR - ALB_WALL_DIF, ALB_WALL_DIR All have dimensions: numrad,numurbl,topounit,lsmlat,lsmlon --- src/dapper/surf/surface_var_specs.py | 64 ++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) diff --git a/src/dapper/surf/surface_var_specs.py b/src/dapper/surf/surface_var_specs.py index 44567e3..826559c 100644 --- a/src/dapper/surf/surface_var_specs.py +++ b/src/dapper/surf/surface_var_specs.py @@ -582,6 +582,70 @@ "units": "gP/m2", "contexts": ["phosphorus_cycle", "topounits"], }, + "ALB_IMPROAD_DIF": { + "dims": "numrad,numurbl,topounit,lsmlat,lsmlon", + "doc": "Diffuse albedo of impervious road; spectral-dependent surface " + "reflectance for urban impervious surfaces.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, + "ALB_IMPROAD_DIR": { + "dims": "numrad,numurbl,topounit,lsmlat,lsmlon", + "doc": "Direct albedo of impervious road; spectral-dependent surface " + "reflectance for urban impervious surfaces.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, + "ALB_PERROAD_DIF": { + "dims": "numrad,numurbl,topounit,lsmlat,lsmlon", + "doc": "Diffuse albedo of pervious road; spectral-dependent surface " + "reflectance for urban pervious surfaces.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, + "ALB_PERROAD_DIR": { + "dims": "numrad,numurbl,topounit,lsmlat,lsmlon", + "doc": "Direct albedo of pervious road; spectral-dependent surface " + "reflectance for urban pervious surfaces.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, + "ALB_ROOF_DIF": { + "dims": "numrad,numurbl,topounit,lsmlat,lsmlon", + "doc": "Diffuse albedo of roof; spectral-dependent surface reflectance " + "for urban roof surfaces.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, + "ALB_ROOF_DIR": { + "dims": "numrad,numurbl,topounit,lsmlat,lsmlon", + "doc": "Direct albedo of roof; spectral-dependent surface reflectance " + "for urban roof surfaces.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, + "ALB_WALL_DIF": { + "dims": "numrad,numurbl,topounit,lsmlat,lsmlon", + "doc": "Diffuse albedo of wall; spectral-dependent surface reflectance " + "for urban wall surfaces.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, + "ALB_WALL_DIR": { + "dims": "numrad,numurbl,topounit,lsmlat,lsmlon", + "doc": "Direct albedo of wall; spectral-dependent surface reflectance " + "for urban wall surfaces.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, } # ----------------------------------------------------------------------------- From 606a8d159cecf729fd1f866dd21083f86afae0e2 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 20 Jul 2026 11:30:01 -0600 Subject: [PATCH 09/14] Add urban thermal, structural, and geometric parameters Add 16 new variables defining urban material properties and geometry: Thermal properties: - CV_* (3): Volumetric heat capacity for impervious road, roof, wall - TK_* (3): Thermal conductivity for impervious road, roof, wall - EM_* (4): Emissivity for impervious/pervious road, roof, wall Structural/geometric properties: - HT_ROOF, THICK_ROOF, THICK_WALL: Roof/wall height and thickness - NLEV_IMPROAD: Number of impervious road layers - T_BUILDING_MAX, T_BUILDING_MIN: Temperature constraints Urban geometry: - CANYON_HWR: Canyon height-to-width ratio - WIND_HGT_CANYON: Wind reference height in canyon - WTLUNIT_ROOF, WTROAD_PERV: Urban landunit weight fractions --- src/dapper/surf/surface_var_specs.py | 154 +++++++++++++++++++++++++++ 1 file changed, 154 insertions(+) diff --git a/src/dapper/surf/surface_var_specs.py b/src/dapper/surf/surface_var_specs.py index 826559c..0389e01 100644 --- a/src/dapper/surf/surface_var_specs.py +++ b/src/dapper/surf/surface_var_specs.py @@ -646,6 +646,160 @@ "units": "unitless", "contexts": ["urban", "topounits"], }, + "CV_IMPROAD": { + "dims": "nlevurb,numurbl,topounit,lsmlat,lsmlon", + "doc": "Volumetric heat capacity of impervious road; thermal mass " + "property affecting diurnal temperature variations.", + "required_level": "optional", + "units": "J/m^3*K", + "contexts": ["urban", "topounits"], + }, + "CV_ROOF": { + "dims": "nlevurb,numurbl,topounit,lsmlat,lsmlon", + "doc": "Volumetric heat capacity of roof; thermal mass property " + "affecting building heat dynamics.", + "required_level": "optional", + "units": "J/m^3*K", + "contexts": ["urban", "topounits"], + }, + "CV_WALL": { + "dims": "nlevurb,numurbl,topounit,lsmlat,lsmlon", + "doc": "Volumetric heat capacity of wall; thermal mass property " + "affecting wall temperature dynamics.", + "required_level": "optional", + "units": "J/m^3*K", + "contexts": ["urban", "topounits"], + }, + "TK_IMPROAD": { + "dims": "nlevurb,numurbl,topounit,lsmlat,lsmlon", + "doc": "Thermal conductivity of impervious road; controls heat diffusion " + "through urban surface layers.", + "required_level": "optional", + "units": "W/m*K", + "contexts": ["urban", "topounits"], + }, + "TK_ROOF": { + "dims": "nlevurb,numurbl,topounit,lsmlat,lsmlon", + "doc": "Thermal conductivity of roof; controls heat transfer through " + "building roof material.", + "required_level": "optional", + "units": "W/m*K", + "contexts": ["urban", "topounits"], + }, + "TK_WALL": { + "dims": "nlevurb,numurbl,topounit,lsmlat,lsmlon", + "doc": "Thermal conductivity of wall; controls heat transfer through " + "building wall material.", + "required_level": "optional", + "units": "W/m*K", + "contexts": ["urban", "topounits"], + }, + "EM_IMPROAD": { + "dims": "numurbl,topounit,lsmlat,lsmlon", + "doc": "Emissivity of impervious road; controls longwave radiation " + "emission from urban surfaces.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, + "EM_PERROAD": { + "dims": "numurbl,topounit,lsmlat,lsmlon", + "doc": "Emissivity of pervious road; controls longwave radiation " + "emission from pervious urban surfaces.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, + "EM_ROOF": { + "dims": "numurbl,topounit,lsmlat,lsmlon", + "doc": "Emissivity of roof; controls longwave radiation emission " + "from building roof.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, + "EM_WALL": { + "dims": "numurbl,topounit,lsmlat,lsmlon", + "doc": "Emissivity of wall; controls longwave radiation emission " + "from building walls.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, + "HT_ROOF": { + "dims": "numurbl,topounit,lsmlat,lsmlon", + "doc": "Height of roof; geometric property defining building structure.", + "required_level": "optional", + "units": "m", + "contexts": ["urban", "topounits"], + }, + "THICK_ROOF": { + "dims": "numurbl,topounit,lsmlat,lsmlon", + "doc": "Thickness of roof; structural property affecting heat capacity.", + "required_level": "optional", + "units": "m", + "contexts": ["urban", "topounits"], + }, + "THICK_WALL": { + "dims": "numurbl,topounit,lsmlat,lsmlon", + "doc": "Thickness of wall; structural property affecting heat capacity.", + "required_level": "optional", + "units": "m", + "contexts": ["urban", "topounits"], + }, + "NLEV_IMPROAD": { + "dims": "numurbl,topounit,lsmlat,lsmlon", + "doc": "Number of impervious road layers; structural discretization for " + "temperature calculation.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, + "T_BUILDING_MAX": { + "dims": "numurbl,topounit,lsmlat,lsmlon", + "doc": "Maximum interior building temperature; upper limit constraint " + "for urban heating/cooling.", + "required_level": "optional", + "units": "K", + "contexts": ["urban", "topounits"], + }, + "T_BUILDING_MIN": { + "dims": "numurbl,topounit,lsmlat,lsmlon", + "doc": "Minimum interior building temperature; lower limit constraint " + "for urban heating/cooling.", + "required_level": "optional", + "units": "K", + "contexts": ["urban", "topounits"], + }, + "CANYON_HWR": { + "dims": "numurbl,topounit,lsmlat,lsmlon", + "doc": "Canyon height to width ratio; urban geometric parameter affecting " + "radiation and wind.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, + "WIND_HGT_CANYON": { + "dims": "numurbl,topounit,lsmlat,lsmlon", + "doc": "Height of wind in canyon; reference height for urban wind profile.", + "required_level": "optional", + "units": "m", + "contexts": ["urban", "topounits"], + }, + "WTLUNIT_ROOF": { + "dims": "numurbl,topounit,lsmlat,lsmlon", + "doc": "Fraction of roof; weight for urban landunit distribution.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, + "WTROAD_PERV": { + "dims": "numurbl,topounit,lsmlat,lsmlon", + "doc": "Fraction of pervious road; weight for pervious surface distribution.", + "required_level": "optional", + "units": "unitless", + "contexts": ["urban", "topounits"], + }, } # ----------------------------------------------------------------------------- From 8e9aff50780cd1246c6c987ac57ccbca8f2b70b7 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 20 Jul 2026 11:30:17 -0600 Subject: [PATCH 10/14] Add hydrologic parameters for water table and runoff modeling Add 8 new variables controlling water table dynamics, infiltration, and surface runoff generation in the Variable Infiltration Capacity (VIC) model: VIC ARNO parameters: - Ds: Fractional saturated area infiltration parameter - Dsmax: Maximum infiltration rate (mm/day) - binfl: Infiltration nonlinearity parameter Saturation/inundation parameters: - F0: Maximum fractional inundated area - FMAX: Maximum fractional saturated area - ZWT0: Water table decay factor Runoff parameters: - P3: Surface runoff lag coefficient - LAKEDEPTH: Average water depth for lake landunit All have dimensions: topounit,lsmlat,lsmlon --- src/dapper/surf/surface_var_specs.py | 58 ++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/src/dapper/surf/surface_var_specs.py b/src/dapper/surf/surface_var_specs.py index 0389e01..70da6d8 100644 --- a/src/dapper/surf/surface_var_specs.py +++ b/src/dapper/surf/surface_var_specs.py @@ -800,6 +800,64 @@ "units": "unitless", "contexts": ["urban", "topounits"], }, + "Ds": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "VIC Ds parameter for the ARNO curve; fractional saturated area " + "infiltration parameter.", + "required_level": "optional", + "units": "unitless", + "contexts": ["hydrology", "topounits"], + }, + "Dsmax": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "VIC Dsmax parameter for the ARNO curve; maximum infiltration rate.", + "required_level": "optional", + "units": "mm/day", + "contexts": ["hydrology", "topounits"], + }, + "F0": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "Maximum gridcell fractional inundated area; controls wetland extent.", + "required_level": "optional", + "units": "unitless", + "contexts": ["hydrology", "inland_water", "topounits"], + }, + "FMAX": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "Maximum fractional saturated area; upper bound on inundation fraction.", + "required_level": "optional", + "units": "unitless", + "contexts": ["hydrology", "topounits"], + }, + "P3": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "Coefficient for qflx_surf_lag for finundated; surface runoff delay parameter.", + "required_level": "optional", + "units": "s/mm", + "contexts": ["hydrology", "topounits"], + }, + "ZWT0": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "Decay factor for finundated; controls inundated area decay with water table.", + "required_level": "optional", + "units": "m", + "contexts": ["hydrology", "topounits"], + }, + "binfl": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "VIC b parameter for the Variable Infiltration Capacity Curve; " + "infiltration nonlinearity.", + "required_level": "optional", + "units": "unitless", + "contexts": ["hydrology", "topounits"], + }, + "LAKEDEPTH": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "Lake depth; average water depth for lake landunit.", + "required_level": "optional", + "units": "m", + "contexts": ["hydrology", "inland_water", "topounits"], + }, } # ----------------------------------------------------------------------------- From 7cc553bb0c9ba7b9e768d898cf1c98f63764bd92 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 20 Jul 2026 11:30:31 -0600 Subject: [PATCH 11/14] Add soil and topographic descriptors Add 4 new variables characterizing soil properties and topographic heterogeneity: Soil classification: - SOIL_COLOR: Categorical soil color affecting albedo - SOIL_ORDER: Categorical soil order for pedogenic classification Topographic properties: - SLP_P10: Slope quantile distribution (nlevslp,topounit,lsmlat,lsmlon) captures topographic variability at multiple percentiles - aveDTB: Average depth to bedrock (topounit,lsmlat,lsmlon) critical for groundwater and subsurface hydrology Enhances subgrid topographic and pedologic characterization for improved process representation in topographic discretization. --- src/dapper/surf/surface_var_specs.py | 29 ++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/src/dapper/surf/surface_var_specs.py b/src/dapper/surf/surface_var_specs.py index 70da6d8..99e7816 100644 --- a/src/dapper/surf/surface_var_specs.py +++ b/src/dapper/surf/surface_var_specs.py @@ -858,6 +858,35 @@ "units": "m", "contexts": ["hydrology", "inland_water", "topounits"], }, + "SOIL_COLOR": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "Soil color; categorical index affecting soil albedo parameterization.", + "required_level": "optional", + "units": "unitless", + "contexts": ["soil_properties", "topounits"], + }, + "SOIL_ORDER": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "Soil order; categorical soil classification for pedogenic properties.", + "required_level": "optional", + "units": "unitless", + "contexts": ["soil_properties", "topounits"], + }, + "SLP_P10": { + "dims": "nlevslp,topounit,lsmlat,lsmlon", + "doc": "Slope at quantiles (minimum and 10 to 100 percentile); " + "topographic distribution parameter.", + "required_level": "optional", + "units": "km km^-1", + "contexts": ["grid_topography", "topounits"], + }, + "aveDTB": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "Average depth to bedrock of the subgrid; critical for groundwater dynamics.", + "required_level": "optional", + "units": "m", + "contexts": ["soil_properties", "grid_topography", "topounits"], + }, } # ----------------------------------------------------------------------------- From 9aa8212debee268ca5e0fa1e5c64c1c79cce77ae Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 20 Jul 2026 11:30:45 -0600 Subject: [PATCH 12/14] Add vegetation isoprene emission factors Add 6 new variables defining basal emission factors for isoprene (EF1) across plant functional types, used in volatile organic compound (VOC) biogeochemistry model: Vegetation-specific emission factors: - EF1_BTR: Broadleaf tree isoprene emission factor - EF1_CRP: Crop isoprene emission factor - EF1_FDT: Deciduous forest isoprene emission factor - EF1_FET: Evergreen forest isoprene emission factor - EF1_GRS: Grass isoprene emission factor - EF1_SHR: Shrub isoprene emission factor All have dimensions: topounit,lsmlat,lsmlon Enables subgrid-resolved VOC and atmospheric chemistry calculations. --- src/dapper/surf/surface_var_specs.py | 48 ++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/src/dapper/surf/surface_var_specs.py b/src/dapper/surf/surface_var_specs.py index 99e7816..e445901 100644 --- a/src/dapper/surf/surface_var_specs.py +++ b/src/dapper/surf/surface_var_specs.py @@ -887,6 +887,54 @@ "units": "m", "contexts": ["soil_properties", "grid_topography", "topounits"], }, + "EF1_BTR": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "EF1 isoprene emission factor for broadleaf tree (BTR); " + "vegetation-dependent emission rate.", + "required_level": "optional", + "units": "unitless", + "contexts": ["vegetation_structure", "biogeochemistry", "topounits"], + }, + "EF1_CRP": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "EF1 isoprene emission factor for crop (CRP); " + "vegetation-dependent emission rate.", + "required_level": "optional", + "units": "unitless", + "contexts": ["vegetation_structure", "biogeochemistry", "topounits"], + }, + "EF1_FDT": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "EF1 isoprene emission factor for deciduous forest (FDT); " + "vegetation-dependent emission rate.", + "required_level": "optional", + "units": "unitless", + "contexts": ["vegetation_structure", "biogeochemistry", "topounits"], + }, + "EF1_FET": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "EF1 isoprene emission factor for evergreen forest (FET); " + "vegetation-dependent emission rate.", + "required_level": "optional", + "units": "unitless", + "contexts": ["vegetation_structure", "biogeochemistry", "topounits"], + }, + "EF1_GRS": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "EF1 isoprene emission factor for grass (GRS); " + "vegetation-dependent emission rate.", + "required_level": "optional", + "units": "unitless", + "contexts": ["vegetation_structure", "biogeochemistry", "topounits"], + }, + "EF1_SHR": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "EF1 isoprene emission factor for shrub (SHR); " + "vegetation-dependent emission rate.", + "required_level": "optional", + "units": "unitless", + "contexts": ["vegetation_structure", "biogeochemistry", "topounits"], + }, } # ----------------------------------------------------------------------------- From 037d44e6784574695645e3394a34d5fde7002b60 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 20 Jul 2026 11:31:01 -0600 Subject: [PATCH 13/14] Add miscellaneous parameters for hydrology, agriculture, and erosion Add 8 new variables spanning hydrologic, disturbance, and erosion processes: Hydrologic VIC parameter: - Ws: VIC maximum soil moisture storage parameter Agricultural/disturbance indicators: - abm: Agricultural fire peak month - peatf: Peatland fraction - gdp: Economic activity indicator Hillslope erosion parameters: - parEro_c1: Rainfall-driven erosion coefficient - parEro_c2: Runoff-driven erosion coefficient - parEro_c3: Sediment transport capacity parameter All have dimensions: topounit,lsmlat,lsmlon Completes topounit discretization for terrestrial and disturbance processes. --- src/dapper/surf/surface_var_specs.py | 52 ++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/src/dapper/surf/surface_var_specs.py b/src/dapper/surf/surface_var_specs.py index e445901..16e1527 100644 --- a/src/dapper/surf/surface_var_specs.py +++ b/src/dapper/surf/surface_var_specs.py @@ -935,6 +935,58 @@ "units": "unitless", "contexts": ["vegetation_structure", "biogeochemistry", "topounits"], }, + "Ws": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "VIC Ws parameter for the ARNO Curve; maximum soil moisture storage.", + "required_level": "optional", + "units": "unitless", + "contexts": ["hydrology", "topounits"], + }, + "abm": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "Agricultural fire peak month; seasonal indicator for fire occurrence.", + "required_level": "optional", + "units": "unitless", + "contexts": ["land_cover", "disturbance", "topounits"], + }, + "gdp": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "GDP (parameter); economic/anthropogenic activity indicator.", + "required_level": "optional", + "units": "unitless", + "contexts": ["anthropogenic", "topounits"], + }, + "peatf": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "Peatland fraction; fraction of grid cell area classified as peat soils.", + "required_level": "optional", + "units": "unitless", + "contexts": ["soil_properties", "land_cover", "topounits"], + }, + "parEro_c1": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "Scalar parameter for rainfall-driven hillslope erosion; " + "controls erosion from precipitation.", + "required_level": "optional", + "units": "unitless", + "contexts": ["grid_topography", "disturbance", "topounits"], + }, + "parEro_c2": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "Scalar parameter for runoff-driven hillslope erosion; " + "controls erosion from surface runoff.", + "required_level": "optional", + "units": "unitless", + "contexts": ["grid_topography", "hydrology", "disturbance", "topounits"], + }, + "parEro_c3": { + "dims": "topounit,lsmlat,lsmlon", + "doc": "Scalar parameter for transport capacity of hillslope overland flow; " + "controls sediment transport capacity.", + "required_level": "optional", + "units": "unitless", + "contexts": ["grid_topography", "hydrology", "disturbance", "topounits"], + }, } # ----------------------------------------------------------------------------- From 248d10a3b9c86e09e2e87672ddb30037714f55de Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 27 Aug 2026 15:47:34 -0600 Subject: [PATCH 14/14] Add template-based filename support for MET file export Enable flexible filename patterns in export_met() by supporting {var} placeholder in the filename parameter. This allows generating E3SM PR #25 compatible filenames (ERA5_{var}_1950-2025_z01.nc) directly without post-processing file renaming. Changes: - Update _nc_filename() to detect {var} placeholder and use .format() - Maintain backward compatibility with simple prefix format - Update docstrings in Domain.export_met() and Exporter.run() Example usage: domain.export_met(..., filename='ERA5_{var}_1950-2025_z01') # Generates: ERA5_TBOT_1950-2025_z01.nc, ERA5_FSDS_1950-2025_z01.nc, etc. Legacy format still works: domain.export_met(..., filename='ERA5') # Generates: ERA5_TBOT.nc, ERA5_FSDS.nc, etc. Co-Authored-By: Claude Sonnet 4.5 --- src/dapper/domains/domain.py | 5 ++++- src/dapper/met/exporter.py | 17 ++++++++++++++--- 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/src/dapper/domains/domain.py b/src/dapper/domains/domain.py index 4ac27f2..e36fc76 100644 --- a/src/dapper/domains/domain.py +++ b/src/dapper/domains/domain.py @@ -974,7 +974,10 @@ def export_met( out_dir : Path-like, optional Override output root. Defaults to Domain.run_dir. filename : str, optional - Optional filename prefix for output NetCDFs. If provided, each var is written to '{filename}_{var}.nc'. + Optional filename template for output NetCDFs. Supports two formats: + - Template with {var} placeholder: 'ERA5_{var}_1950-2025_z01' generates 'ERA5_TBOT_1950-2025_z01.nc' + - Simple prefix (legacy): 'prefix' generates 'prefix_TBOT.nc' + The '.nc' extension is added automatically. overwrite : bool If False, raises if MET output(s) already exist. clip_to_full_years : bool or None diff --git a/src/dapper/met/exporter.py b/src/dapper/met/exporter.py index d81e93e..a6e693b 100644 --- a/src/dapper/met/exporter.py +++ b/src/dapper/met/exporter.py @@ -239,8 +239,10 @@ def run(self, *, pack_scope=None, filename: str | None = None, overwrite: bool = ``global`` for cellset outputs. filename - Optional filename prefix for output NetCDF files. If provided, each variable - is written to ``{filename}_{var}.nc``. + Optional filename template for output NetCDF files. Supports two formats: + - Template with {var} placeholder: 'ERA5_{var}_1950-2025_z01' generates 'ERA5_TBOT_1950-2025_z01.nc' + - Simple prefix (legacy): 'prefix' generates 'prefix_TBOT.nc' + The '.nc' extension is added automatically. overwrite If True, clears existing MET outputs before writing. @@ -343,8 +345,17 @@ def _zone_mappings_path(self, gid: str | None = None, filename: str = "zone_mapp return self._met_dir_for_gid(gid) / filename def _nc_filename(self, var: str) -> str: - """Return the output NetCDF filename for a given variable.""" + """Return the output NetCDF filename for a given variable. + + Supports two formats: + - Template with {var} placeholder: 'ERA5_{var}_1950-2025_z01' -> 'ERA5_TBOT_1950-2025_z01.nc' + - Simple prefix (legacy): 'prefix' -> 'prefix_TBOT.nc' + """ if getattr(self, "filename_prefix", None): + # Check if template contains {var} placeholder + if "{var}" in self.filename_prefix: + return f"{self.filename_prefix.format(var=var)}.nc" + # Legacy prefix mode return f"{self.filename_prefix}_{var}.nc" return f"{var}.nc"