Preserve exact native grids when combining soil auxiliaries

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

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)

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,22 +126,25 @@ 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:
                         if dim in result.variables and dim not in {'time','layer'}:
                             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:
         temporary.unlink(missing_ok=True)
     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
 

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:

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']

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)