From 020a2e0d446a874158cb4c71a03a89d7112b0de6 Mon Sep 17 00:00:00 2001 From: Yi Yu <72683152+yuyi13@users.noreply.github.com> Date: Sun, 18 Oct 2026 05:58:58 +1100 Subject: [PATCH] Preserve native auxiliary affine metadata and strict final grid checks --- docs/output-netcdf-layout.md | 9 ++++ pysweb/io/products.py | 54 ++++++++++++++++++-- pysweb/swb/paired.py | 12 +++-- pysweb/swb/preprocess.py | 1 + tests/io/test_native_grids.py | 67 ++++++++++++++++++++++++- tests/workflows/test_native_workflow.py | 2 + 6 files changed, 135 insertions(+), 10 deletions(-) diff --git a/docs/output-netcdf-layout.md b/docs/output-netcdf-layout.md index 9d2a5fe..3c8f3e6 100644 --- a/docs/output-netcdf-layout.md +++ b/docs/output-netcdf-layout.md @@ -17,6 +17,15 @@ Scene retrieval auxiliaries use `scene_time`; dT climatology uses `dayofyear`. T Intermediate single-scenario SWB files retain the component API's `rzsm_layer_*` variables under `products/intermediate/`. They support diagnosis and are not additional consolidated final products. Native calibration inputs remain separate under `calibration_inputs/`; they intentionally have a different spatial grid. See [the full native-grid contract](native-grids.md). +Final grid checks compare the declared affine metadata and coordinate arrays exactly. +Auxiliaries must have exactly equal pixel centres and CRS before copying. Affine +differences confined to floating-point coordinate roundoff are recorded and +normalised to the output grid mapping, without resampling or changing values. +Larger affine differences, shifted coordinates and different CRS are rejected. +The NetCDF attribute `auxiliary_grid_metadata_json` and its provenance manifest +retain the source and output transforms and any adjustment. Newly prepared soil +files preserve the selected template affine directly. + ## Streaming large records Streamed daily and scene-dependent arrays use temporal chunks of at most diff --git a/pysweb/io/products.py b/pysweb/io/products.py index 4028939..618d90e 100644 --- a/pysweb/io/products.py +++ b/pysweb/io/products.py @@ -1,36 +1,80 @@ """ Script: products.py -Objective: Publish exact spatial subsets and verify the shared final product grid. +Objective: Publish exact spatial subsets and verify product and auxiliary grid contracts. Author: Yi Yu (with assistance from Codex) Created: 2026-10-17 Last updated: 2026-10-17 Inputs: NetCDF source products and an exact Landsat template window. -Outputs: Streaming NetCDF subsets and strict grid-contract validation. +Outputs: Streaming subsets, strict final grid checks and auxiliary metadata provenance. Usage: from pysweb.io.products import crop_product, assert_identical_grids -Dependencies: netCDF4, numpy, pysweb +Dependencies: netCDF4, numpy, rasterio, xarray, pysweb """ from pathlib import Path import os import netCDF4 as nc import numpy as np +from rasterio import Affine +import xarray as xr from pysweb.io.grids import open_template, spatial from pysweb.io.provenance import input_manifest, save_output_manifest from pysweb.io.netcdf_stream import stream_chunks, sequential_read_cache +def _declared_transform(path, template): + """Read the stored affine too; coordinate inference can hide a bad origin.""" + if Path(path).suffix.lower() not in {'.tif', '.tiff'}: + with xr.open_dataset(path, decode_coords='all') as dataset: + mapping = dataset.rio.grid_mapping + if mapping in dataset and 'GeoTransform' in dataset[mapping].attrs: + values = [float(value) for value in dataset[mapping].attrs['GeoTransform'].split()] + if len(values) != 6 or not np.isfinite(values).all(): + raise ValueError('Invalid declared affine transform') + return Affine.from_gdal(*values) + return template.rio.transform() + + def assert_identical_grids(first_path, second_path): first, second = open_template(first_path), open_template(second_path) if first.rio.crs != second.rio.crs: raise ValueError('Final product CRS differs') - if not np.array_equal(tuple(first.rio.transform()), tuple(second.rio.transform())): + transform = _declared_transform(first_path, first) + if not np.array_equal(tuple(transform), tuple(_declared_transform(second_path, second))): raise ValueError('Final product affine transform differs') for dim in ('x', 'y'): if not np.array_equal(first[dim].values, second[dim].values): raise ValueError(f'Final product {dim} coordinates differ') - return {'crs': str(first.rio.crs), 'transform': list(first.rio.transform())[:6], + return {'crs': str(first.rio.crs), 'transform': list(transform)[:6], 'shape': [first.sizes['y'], first.sizes['x']], 'identical': True} +def auxiliary_grid_contract(output_path, auxiliary_path): + """Require exact centres and CRS before copying onto the output metadata. + + Reconstructing an affine from floating-point centres can change its last + bits. Only coordinate-scale roundoff is admissible for an auxiliary whose + centres match exactly. Final product comparisons remain bitwise strict. + """ + output, auxiliary = open_template(output_path), open_template(auxiliary_path) + if output.rio.crs != auxiliary.rio.crs: + raise ValueError('Auxiliary CRS differs') + for dim in ('x', 'y'): + if not np.array_equal(output[dim].values, auxiliary[dim].values): + raise ValueError(f'Auxiliary {dim} coordinates differ') + target = np.asarray(_declared_transform(output_path, output))[:6] + source = np.asarray(_declared_transform(auxiliary_path, auxiliary))[:6] + scale = max(1., *(float(np.max(np.abs(output[dim].values))) for dim in ('x', 'y'))) + tolerance = 8*np.finfo('float64').eps*scale + if (target[[1, 3]] != 0).any() or (source[[1, 3]] != 0).any(): + raise ValueError('Auxiliary rotated affine is unsupported') + if not np.isfinite(source).all() or not np.allclose(source, target, rtol=0, atol=tolerance): + raise ValueError('Auxiliary affine differs beyond coordinate roundoff') + return dict(identical_coordinates=True, crs=str(output.rio.crs), + source_transform=source.tolist(), output_transform=target.tolist(), + affine_metadata_normalized=not np.array_equal(source, target), + maximum_affine_adjustment=float(np.max(np.abs(source-target))), + coordinate_roundoff_tolerance=tolerance) + + def crop_product(source_path, output_path, template): """Copy every auxiliary in bounded slabs; crop only exact x/y source cells.""" template = spatial(template) diff --git a/pysweb/swb/paired.py b/pysweb/swb/paired.py index 0d654b5..aed236d 100644 --- a/pysweb/swb/paired.py +++ b/pysweb/swb/paired.py @@ -17,7 +17,7 @@ from pysweb.io.provenance import input_manifest, save_output_manifest from pysweb.swb.run import run_swb_workflow from pysweb.swb.core import SOIL_FILE_STEMS -from pysweb.io.products import assert_identical_grids +from pysweb.io.products import assert_identical_grids, auxiliary_grid_contract from pysweb.io.netcdf_stream import stream_chunks, sequential_read_cache @@ -73,6 +73,7 @@ def combine_scenarios(scenarios, output, auxiliaries=None, config=None): output.parent.mkdir(parents=True, exist_ok=True) temporary = output.with_name(f'.{output.name}.{os.getpid()}.tmp') auxiliaries = auxiliaries or {} + auxiliary_grids = {} first_path = next(iter(scenarios.values())) for path in scenarios.values(): assert_identical_grids(first_path,path) @@ -95,7 +96,7 @@ def combine_scenarios(scenarios, output, auxiliaries=None, config=None): depth[:] = [first[name].depth_bottom_mm for name in layers] depth.units = 'mm' result.setncatts({k:first.getncattr(k) for k in first.ncattrs()}) - result.title = 'SWB default and calibrated soil moisture on the Landsat grid' + result.title = 'SWB default and calibrated soil moisture' result.scenarios = ','.join(scenarios) result.configuration_json = json.dumps(config or {}, default=str, sort_keys=True) for name in ('precipitation', 'effective_precipitation'): @@ -125,7 +126,7 @@ def combine_scenarios(scenarios, output, auxiliaries=None, config=None): if key in source.ncattrs(): result.setncattr(f'{key}_{scenario}',source.getncattr(key)) for name,(path,variable_name) in auxiliaries.items(): - assert_identical_grids(first_path,path) + auxiliary_grids[name] = auxiliary_grid_contract(first_path,path) with nc.Dataset(path) as source: variable = source[variable_name] for dim in variable.dimensions: @@ -133,7 +134,9 @@ def combine_scenarios(scenarios, output, auxiliaries=None, config=None): if not np.array_equal(source[dim][:], result[dim][:]): raise ValueError(f'Auxiliary {name} is not on the output grid') indices = _time_indices(source,result) if 'time' in variable.dimensions else None - _copy_variable(result,variable,name,indices) + copied = _copy_variable(result,variable,name,indices) + copied.setncattr('grid_mapping', 'spatial_ref') + result.auxiliary_grid_metadata_json = json.dumps(auxiliary_grids, sort_keys=True) result.Conventions = 'CF-1.8' temporary.replace(output) finally: @@ -141,6 +144,7 @@ def combine_scenarios(scenarios, output, auxiliaries=None, config=None): paths = list(scenarios.values())+[path for path,_ in auxiliaries.values()] manifest = input_manifest(config or {}, paths) manifest['scenarios'] = {key:str(value) for key,value in scenarios.items()} + manifest['auxiliary_grids'] = auxiliary_grids save_output_manifest(output,manifest) return output diff --git a/pysweb/swb/preprocess.py b/pysweb/swb/preprocess.py index ded0420..f291913 100644 --- a/pysweb/swb/preprocess.py +++ b/pysweb/swb/preprocess.py @@ -1285,6 +1285,7 @@ def preprocess_inputs(**kwargs): reproject_to_template = _reproject_to_template, ) soil_arrays = {name:da.rio.write_crs(grid.crs).rio.set_spatial_dims(x_dim=grid.lon_dim,y_dim=grid.lat_dim) + .rio.write_transform(grid.template.rio.transform()) for name,da in soil_outputs.arrays.items()} aoi = None if args.template_extent is not None: diff --git a/tests/io/test_native_grids.py b/tests/io/test_native_grids.py index e7d1b9a..2515329 100644 --- a/tests/io/test_native_grids.py +++ b/tests/io/test_native_grids.py @@ -18,7 +18,7 @@ from rasterio.crs import CRS from pysweb.io.native_grid import verify_native_grid, native_projection, export_parameters from pysweb.io.grids import align_to_grid, subset_template -from pysweb.io.products import crop_product, assert_identical_grids +from pysweb.io.products import crop_product, assert_identical_grids, auxiliary_grid_contract def raster(values, transform=Affine(30,0,500000,0,-30,6100000), crs='EPSG:32632'): @@ -99,3 +99,68 @@ def test_global_origin_download_uses_bounded_native_pixel_window(): assert 'region' not in request and 'scale' not in request grid = SimpleNamespace(crs=CRS.from_epsg(4326),transform=Affine(*request['crs_transform']),width=73,height=73) assert verify_native_grid(grid,projection)['native_grid_verified'] + + +def test_auxiliary_affine_roundoff_keeps_final_contract_strict(tmp_path): + transform = Affine(.008983152841195215,0,146.25471140749926, + 0,-.008983152841195215,-34.96243085793178) + source = raster(np.ones((6,8)),transform,'EPSG:4326').rename('field') + derived = Affine(float((source.x[-1]-source.x[0])/7),0,transform.c, + 0,float((source.y[-1]-source.y[0])/5),transform.f) + source.to_netcdf(tmp_path/'canonical.nc') + source.rio.write_transform(derived).to_netcdf(tmp_path/'derived.nc') + contract = auxiliary_grid_contract(tmp_path/'canonical.nc',tmp_path/'derived.nc') + assert contract['identical_coordinates'] and contract['affine_metadata_normalized'] + assert 0 < contract['maximum_affine_adjustment'] < 1e-13 + with pytest.raises(ValueError,match='affine'): + assert_identical_grids(tmp_path/'canonical.nc',tmp_path/'derived.nc') + + +@pytest.mark.parametrize('change', ['origin', 'scale', 'coordinate', 'crs']) +def test_auxiliary_grid_rejects_mismatch_despite_metadata_normalization(tmp_path,change): + source = raster(np.ones((3,4))).rename('field') + transform = source.rio.transform() + source.to_netcdf(tmp_path/'canonical.nc') + if change == 'origin': + wrong = source.rio.write_transform(Affine(transform.a,0,transform.c+.01, + 0,transform.e,transform.f)) + elif change == 'scale': + wrong = source.rio.write_transform(Affine(transform.a+.01,0,transform.c, + 0,transform.e,transform.f)) + elif change == 'coordinate': + wrong = source.assign_coords(x=np.nextafter(source.x.values,np.inf)) + else: + wrong = source.rio.write_crs('EPSG:32633') + wrong.to_netcdf(tmp_path/'wrong.nc') + with pytest.raises(ValueError): + auxiliary_grid_contract(tmp_path/'canonical.nc',tmp_path/'wrong.nc') + with pytest.raises(ValueError): + assert_identical_grids(tmp_path/'canonical.nc',tmp_path/'wrong.nc') + + +def test_combined_auxiliary_uses_exact_output_metadata_and_preserves_values(tmp_path): + import json + from pysweb.swb.paired import combine_scenarios + transform = Affine(.008983152841195215,0,146.25471140749926, + 0,-.008983152841195215,-34.96243085793178) + grid = raster(np.arange(48).reshape(6,8),transform,'EPSG:4326') + days = np.array(['2014-01-01','2014-01-02'],dtype='datetime64[D]') + daily = grid.expand_dims(time=days) + scenario = xr.Dataset(dict(rzsm_layer_1=daily/200, + precipitation=daily*0,effective_precipitation=daily*0)) + scenario.rzsm_layer_1.attrs['depth_bottom_mm'] = 50. + scenario.attrs['state_timing'] = 'start' + scenario.rio.write_crs(grid.rio.crs).rio.write_transform(transform).to_netcdf(tmp_path/'scenario.nc') + derived = Affine(float((grid.x[-1]-grid.x[0])/7),0,transform.c, + 0,float((grid.y[-1]-grid.y[0])/5),transform.f) + auxiliary = grid.rename('porosity').rio.write_transform(derived) + auxiliary.to_netcdf(tmp_path/'soil.nc') + result = combine_scenarios({'default':tmp_path/'scenario.nc'},tmp_path/'combined.nc', + {'soil_porosity':(tmp_path/'soil.nc','porosity')}) + assert assert_identical_grids(result,tmp_path/'scenario.nc')['identical'] + with xr.open_dataset(result) as actual: + np.testing.assert_array_equal(actual.soil_porosity.values,auxiliary.values) + assert actual.soil_porosity.attrs['grid_mapping'] == 'spatial_ref' + assert json.loads(actual.auxiliary_grid_metadata_json)['soil_porosity']['affine_metadata_normalized'] + manifest = json.loads((tmp_path/'combined.nc.manifest.json').read_text()) + assert manifest['auxiliary_grids']['soil_porosity']['identical_coordinates'] diff --git a/tests/workflows/test_native_workflow.py b/tests/workflows/test_native_workflow.py index bdf4cbc..4f72088 100644 --- a/tests/workflows/test_native_workflow.py +++ b/tests/workflows/test_native_workflow.py @@ -91,6 +91,8 @@ def test_native_full_offline_workflow(tmp_path): with xr.open_dataset(base/'calibration_inputs/reference_ssm_daily_20200101_20200104.nc') as coarse: assert coarse.sizes['lat']*coarse.sizes['lon'] < 10 np.testing.assert_allclose(coarse.reference_ssm,.25) + assert assert_identical_grids(base/'calibration_inputs/effective_precip_daily_20200101_20200104.nc', + base/'calibration_inputs/soil_porosity.nc')['identical'] parameters = pd.read_csv(base/'calibration.csv') assert len(parameters) == 1 assert set(['diff_factor','sm_max_factor','sm_min_factor','root_beta']) <= set(parameters)