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/landuse/landuse.py b/src/dapper/landuse/landuse.py index ba3997c..d31ad33 100644 --- a/src/dapper/landuse/landuse.py +++ b/src/dapper/landuse/landuse.py @@ -12,6 +12,7 @@ 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"] @@ -270,6 +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_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/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" diff --git a/src/dapper/surf/fraction_closure.py b/src/dapper/surf/fraction_closure.py new file mode 100644 index 0000000..afcaf18 --- /dev/null +++ b/src/dapper/surf/fraction_closure.py @@ -0,0 +1,209 @@ +"""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), + ) + + # 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") + + 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", + "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 0ac4f42..6f96e42 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 = ds.copy(deep=False) + ds2 = normalize_fraction_closure(ds) merged = dict(ds2.attrs) if dapper_attrs: @@ -202,11 +203,16 @@ def write_surface_nc( ds2.attrs = merged + 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": - 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 @@ -886,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) @@ -1021,11 +1027,19 @@ 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. + 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) @@ -1091,8 +1105,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] @@ -1108,7 +1122,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/src/dapper/surf/surface_var_specs.py b/src/dapper/surf/surface_var_specs.py index cacf037..16e1527 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,472 @@ "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"], + }, + "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"], + }, + "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"], + }, + "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"], + }, + "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"], + }, + "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"], + }, + "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"], }, } diff --git a/src/dapper/surf/validate.py b/src/dapper/surf/validate.py index 5ffc15f..05cfdd2 100644 --- a/src/dapper/surf/validate.py +++ b/src/dapper/surf/validate.py @@ -320,13 +320,42 @@ 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}")) + + # 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}")) + + # 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: diff --git a/tests/test_fraction_closure.py b/tests/test_fraction_closure.py new file mode 100644 index 0000000..2a05c97 --- /dev/null +++ b/tests/test_fraction_closure.py @@ -0,0 +1,312 @@ +import numpy as np +import xarray as xr + +from dapper.surf.sfile import write_surface_nc +from dapper.surf.fraction_closure import normalize_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_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() + + +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_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), + "lsmlat": np.arange(1, dtype=np.int32), + "lsmlon": np.arange(1, dtype=np.int32), + } + ) + + 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["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)