From 64d7b9c1609c48c892286c924cff4cca690c7d69 Mon Sep 17 00:00:00 2001 From: Arthur Jacquiau-Chamski Date: Wed, 1 Apr 2026 23:29:36 +0200 Subject: [PATCH 01/10] Add tests for TerracottaBackend to ensure correct coords projection when writing to Terracotta --- tests/backends/test_terracotta_alignment.py | 155 ++++++++++++++++++++ 1 file changed, 155 insertions(+) create mode 100644 tests/backends/test_terracotta_alignment.py diff --git a/tests/backends/test_terracotta_alignment.py b/tests/backends/test_terracotta_alignment.py new file mode 100644 index 0000000..da1789f --- /dev/null +++ b/tests/backends/test_terracotta_alignment.py @@ -0,0 +1,155 @@ +import os +from pathlib import Path + +import fsspec +import numpy as np +import pytest +import rasterio +import xarray as xr +from pyproj import Transformer +from rasterio.transform import rowcol, xy +from shapely.geometry import box, mapping +from terracotta import get_driver +from rasterio.enums import Resampling + +from nuthatch.backends.terracotta import TerracottaBackend + +GRID_RESOLUTIONS = [1.5, 0.25, 0.1] +WEB_MERCATOR = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True) + + +def _northup_source_dataset(resolution_degrees): + # Regular lat/lon grid with global longitude coverage and Mercator-safe latitudes. + half_res = resolution_degrees / 2.0 + lons = np.arange(-180.0 + half_res, 180.0, resolution_degrees) + lats = np.arange(84.0 - half_res, -84.0, -resolution_degrees) + values = np.arange(len(lats) * len(lons), dtype=float).reshape(len(lats), len(lons)) + return xr.Dataset( + {"precip": (("lat", "lon"), values)}, + coords={"lat": lats, "lon": lons}, + ) + +AFRICA_BOUNDS = (-19.5, -40.5, 55.5, 40.5) +def _clip_source_to_pseudo_africa(ds): + return ( + ds.rio.write_crs("EPSG:4326") + .rio.set_spatial_dims("lon", "lat") + .rio.clip( + [mapping(box(*AFRICA_BOUNDS))], + "EPSG:4326", + drop=True, + ) + ) + + +def _make_test_backend(tmp_path: Path, cache_key: str): + backend = TerracottaBackend.__new__(TerracottaBackend) + backend.lat_dim = "lat" + backend.lon_dim = "lon" + backend.time_dim = "time" + backend.resampling = Resampling.nearest + backend.cache_key = cache_key + backend.path = str(tmp_path / f"{cache_key}.terracotta") + backend.override_path = backend.path + backend.fs = fsspec.filesystem("file") + Path(backend.path).mkdir(parents=True, exist_ok=True) + backend.driver = get_driver(str(tmp_path / "terracotta_scope.sqlite")) + try: + backend.driver.get_keys() + except Exception: + backend.driver.create(["key"]) + return backend + + +def _reference_points_from_dataset(ds): + lon_indices = [ds.sizes["lon"] // 6, ds.sizes["lon"] // 2, (5 * ds.sizes["lon"]) // 6] + lat_indices = [ds.sizes["lat"] // 6, ds.sizes["lat"] // 2, (5 * ds.sizes["lat"]) // 6] + return [ + (float(ds.lon.values[lon_idx]), float(ds.lat.values[lat_idx])) + for lon_idx, lat_idx in zip(lon_indices, lat_indices, strict=False) + ] + + +@pytest.mark.parametrize("resolution_degrees", GRID_RESOLUTIONS) +def test_terracotta_write_preserves_latlon_resolution_for_crops(tmp_path, resolution_degrees): + source = _northup_source_dataset(resolution_degrees) + clipped = _clip_source_to_pseudo_africa(source) + + global_backend = _make_test_backend(tmp_path, "global") + africa_backend = _make_test_backend(tmp_path, "africa") + global_backend.write(source) + africa_backend.write(clipped) + + global_path = tmp_path / "global.terracotta" / "_.tif" + africa_path = tmp_path / "africa.terracotta" / "_.tif" + + with ( + rasterio.open(global_path) as global_tif, + rasterio.open(africa_path) as africa_tif, + ): + assert global_tif.crs == africa_tif.crs + assert global_tif.res == pytest.approx(africa_tif.res) + +@pytest.mark.parametrize("resolution_degrees", GRID_RESOLUTIONS) +def test_terracotta_returns_same_tile_values_for_global_and_africa_scopes(tmp_path, resolution_degrees): + source = _northup_source_dataset(resolution_degrees) + clipped = _clip_source_to_pseudo_africa(source) + + global_backend = _make_test_backend(tmp_path, "global") + africa_backend = _make_test_backend(tmp_path, "africa") + global_backend.write(source) + africa_backend.write(clipped) + driver = global_backend.driver + + tile_bounds = driver.get_metadata(("africa",))["bounds"] + global_tile = driver.get_raster_tile( + {"key": "global"}, + tile_bounds=tile_bounds, + tile_size=(256, 256), + preserve_values=True, + ) + africa_tile = driver.get_raster_tile( + {"key": "africa"}, + tile_bounds=tile_bounds, + tile_size=(256, 256), + preserve_values=True, + ) + + assert np.array_equal(global_tile.filled(-9999), africa_tile.filled(-9999)) + + +@pytest.mark.parametrize("resolution_degrees", GRID_RESOLUTIONS) +def test_terracotta_write_projects_reference_points_consistently(tmp_path, resolution_degrees): + source = _northup_source_dataset(resolution_degrees) + clipped = _clip_source_to_pseudo_africa(source) + + global_backend = _make_test_backend(tmp_path, "global") + africa_backend = _make_test_backend(tmp_path, "africa") + global_backend.write(source) + africa_backend.write(clipped) + + reference_points = _reference_points_from_dataset(clipped) + + global_path = tmp_path / "global.terracotta" / "_.tif" + africa_path = tmp_path / "africa.terracotta" / "_.tif" + + with ( + rasterio.open(global_path) as global_tif, + rasterio.open(africa_path) as africa_tif, + ): + for lon, lat in reference_points: + expected_x, expected_y = WEB_MERCATOR.transform(lon, lat) + observed_centers = [] + + for tif in (global_tif, africa_tif): + row, col = rowcol(tif.transform, expected_x, expected_y) + assert 0 <= row < tif.height + assert 0 <= col < tif.width + + observed_x, observed_y = xy(tif.transform, row, col, offset="center") + observed_centers.append((observed_x, observed_y)) + assert observed_x == pytest.approx(expected_x, abs=tif.res[0] / 2.0) + assert observed_y == pytest.approx(expected_y, abs=tif.res[1] / 2.0) + + assert observed_centers[0][0] == pytest.approx(observed_centers[1][0]) + assert observed_centers[0][1] == pytest.approx(observed_centers[1][1]) From 8c0232add94c1bb82777f55ade052ca42abf6950 Mon Sep 17 00:00:00 2001 From: Arthur Jacquiau-Chamski Date: Wed, 1 Apr 2026 23:38:28 +0200 Subject: [PATCH 02/10] fix underspecified Mercator projection in TerracottaBackend --- src/nuthatch/backends/terracotta.py | 102 +++++++++++++++++++++++++++- 1 file changed, 99 insertions(+), 3 deletions(-) diff --git a/src/nuthatch/backends/terracotta.py b/src/nuthatch/backends/terracotta.py index 0db6ed0..6af066e 100644 --- a/src/nuthatch/backends/terracotta.py +++ b/src/nuthatch/backends/terracotta.py @@ -5,8 +5,11 @@ import xarray as xr import rioxarray # Must import for .rio to work # noqa: F401 import numpy as np +from pyproj import CRS from rasterio.io import MemoryFile from rasterio.enums import Resampling +from rasterio.transform import Affine +from rasterio.warp import aligned_target, calculate_default_transform, transform from nuthatch.backend import DatabaseBackend, FileBackend, register_backend import logging @@ -93,6 +96,92 @@ def lon_base_change(ds, to_base="base180", lon_dim='lon'): return ds +def _get_regular_resolution(coords, dim_name): + """Infer the resolution of a regularly spaced coordinate axis.""" + values = np.asarray(coords.values, dtype=float) + if values.size < 2: + raise ValueError(f"Need at least two coordinates along {dim_name} to infer resolution.") + + diffs = np.diff(values) + resolution = float(np.abs(diffs[0])) + if resolution == 0.0: + raise ValueError(f"Coordinate spacing along {dim_name} must be non-zero.") + if not np.allclose(np.abs(diffs), resolution): + raise ValueError(f"Coordinates along {dim_name} must be regularly spaced.") + return resolution + + +def _get_aligned_mercator_target(ds): + """Build a scope-independent Web Mercator target for this raster resolution.""" + # Convert the source cell spacing to an approximate target spacing in Web Mercator meters. + projected_x, projected_y = transform( + "EPSG:4326", + "EPSG:3857", + [0.0, _get_regular_resolution(ds.x, "x")], + [0.0, _get_regular_resolution(ds.y, "y")], + ) + target_resolution = ( + abs(float(projected_x[1] - projected_x[0])), + abs(float(projected_y[1] - projected_y[0])), + ) + + # Ask rasterio for a projected target grid at that resolution, then snap it to aligned pixel edges. + target_transform, width, height = calculate_default_transform( + "EPSG:4326", + "EPSG:3857", + ds.sizes["x"], + ds.sizes["y"], + *ds.rio.bounds(), + resolution=target_resolution, + ) + target_transform, width, height = aligned_target( + target_transform, + width, + height, + target_resolution, + ) + # Derive the valid projected extent directly from the CRS area of use. + area_of_use = CRS.from_epsg(3857).area_of_use + valid_x, valid_y = transform( + "EPSG:4326", + "EPSG:3857", + [area_of_use.west, area_of_use.east], + [area_of_use.south, area_of_use.north], + ) + x_resolution, y_resolution = target_resolution + left = target_transform.c + top = target_transform.f + valid_left, valid_right = min(valid_x), max(valid_x) + valid_bottom, valid_top = min(valid_y), max(valid_y) + right = left + width * x_resolution + bottom = top - height * y_resolution + + # aligned_target can expand the raster past Web Mercator's valid extent, so trim whole rows/columns back. + if left < valid_left: + trim = int(np.ceil((valid_left - left) / x_resolution)) + left += trim * x_resolution + width -= trim + + if right > valid_right: + trim = int(np.ceil((right - valid_right) / x_resolution)) + width -= trim + + if top > valid_top: + trim = int(np.ceil((top - valid_top) / y_resolution)) + top -= trim * y_resolution + height -= trim + + if bottom < valid_bottom: + trim = int(np.ceil((valid_bottom - bottom) / y_resolution)) + height -= trim + + if width <= 0 or height <= 0: + raise ValueError("Aligned raster window fell outside the valid extent for EPSG:3857.") + + # Return a rasterio affine and shape that can be passed straight into rio.reproject(...). + return Affine(x_resolution, 0.0, left, 0.0, -y_resolution, top), width, height + + @register_backend class TerracottaBackend(DatabaseBackend, FileBackend): """ @@ -188,9 +277,17 @@ def write(self, ds, upsert=False, primary_keys=None): else: ds = ds.transpose('y', 'x') - # Adapt the CRS + # Adapt the CRS to Web Mercator. ds.rio.write_crs("epsg:4326", inplace=True) - ds = ds.rio.reproject('EPSG:3857', resampling=self.resampling, nodata=np.nan) + ds.rio.set_spatial_dims("x", "y", inplace=True) + mercator_transform, width, height = _get_aligned_mercator_target(ds) + ds = ds.rio.reproject( + "EPSG:3857", + transform=mercator_transform, + shape=(height, width), + resampling=self.resampling, + nodata=np.nan, + ) ds.rio.write_crs("epsg:3857", inplace=True) # Insert the parameters. @@ -268,4 +365,3 @@ def delete(self): for dataset in datasets: logger.info(f"Deleting datasets {datasets} from terracotta.") self.driver.delete({'key': dataset}) - From 281f333c1e06f0aef66c2fbea4f23e1a7847f400 Mon Sep 17 00:00:00 2001 From: Arthur Jacquiau-Chamski Date: Wed, 1 Apr 2026 23:51:25 +0200 Subject: [PATCH 03/10] fix/test: update SQLite URL format for Terracotta backend driver initialization in tests --- tests/backends/test_terracotta_alignment.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/backends/test_terracotta_alignment.py b/tests/backends/test_terracotta_alignment.py index da1789f..5aa7d40 100644 --- a/tests/backends/test_terracotta_alignment.py +++ b/tests/backends/test_terracotta_alignment.py @@ -53,7 +53,8 @@ def _make_test_backend(tmp_path: Path, cache_key: str): backend.override_path = backend.path backend.fs = fsspec.filesystem("file") Path(backend.path).mkdir(parents=True, exist_ok=True) - backend.driver = get_driver(str(tmp_path / "terracotta_scope.sqlite")) + sqlite_url = f"sqlite:///{(tmp_path / 'terracotta_scope.sqlite').as_posix()}" + backend.driver = get_driver(sqlite_url) try: backend.driver.get_keys() except Exception: From 2923326951cd042958e20fe2e1312eb69b33fc3a Mon Sep 17 00:00:00 2001 From: Arthur Jacquiau-Chamski Date: Thu, 2 Apr 2026 00:00:36 +0200 Subject: [PATCH 04/10] fix/test: update SQLite URL format for Terracotta backend driver initialization in tests --- tests/backends/test_terracotta_alignment.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/backends/test_terracotta_alignment.py b/tests/backends/test_terracotta_alignment.py index 5aa7d40..da69aa6 100644 --- a/tests/backends/test_terracotta_alignment.py +++ b/tests/backends/test_terracotta_alignment.py @@ -53,8 +53,7 @@ def _make_test_backend(tmp_path: Path, cache_key: str): backend.override_path = backend.path backend.fs = fsspec.filesystem("file") Path(backend.path).mkdir(parents=True, exist_ok=True) - sqlite_url = f"sqlite:///{(tmp_path / 'terracotta_scope.sqlite').as_posix()}" - backend.driver = get_driver(sqlite_url) + backend.driver = get_driver(tmp_path / "terracotta_scope.sqlite", provider="sqlite") try: backend.driver.get_keys() except Exception: From 8c370f14e2567e90dd2934bba6bc09474288c9d0 Mon Sep 17 00:00:00 2001 From: Arthur Jacquiau-Chamski Date: Thu, 2 Apr 2026 00:12:32 +0200 Subject: [PATCH 05/10] fix/test: update SQLite URL format for Terracotta backend driver initialization in tests --- tests/backends/test_terracotta_alignment.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/tests/backends/test_terracotta_alignment.py b/tests/backends/test_terracotta_alignment.py index da69aa6..e5767b0 100644 --- a/tests/backends/test_terracotta_alignment.py +++ b/tests/backends/test_terracotta_alignment.py @@ -42,6 +42,10 @@ def _clip_source_to_pseudo_africa(ds): ) +def _sqlite_database_path(path: Path) -> str: + return str(path) + + def _make_test_backend(tmp_path: Path, cache_key: str): backend = TerracottaBackend.__new__(TerracottaBackend) backend.lat_dim = "lat" @@ -53,7 +57,10 @@ def _make_test_backend(tmp_path: Path, cache_key: str): backend.override_path = backend.path backend.fs = fsspec.filesystem("file") Path(backend.path).mkdir(parents=True, exist_ok=True) - backend.driver = get_driver(tmp_path / "terracotta_scope.sqlite", provider="sqlite") + backend.driver = get_driver( + _sqlite_database_path(tmp_path / "terracotta_scope.sqlite"), + provider="sqlite", + ) try: backend.driver.get_keys() except Exception: From 9251f75a8ca7a398c3239dd5a7ef6b178eca755b Mon Sep 17 00:00:00 2001 From: Arthur Jacquiau-Chamski Date: Thu, 2 Apr 2026 00:20:30 +0200 Subject: [PATCH 06/10] fix/test: update SQLite URL format for Terracotta backend driver initialization in tests --- tests/backends/test_terracotta_alignment.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/backends/test_terracotta_alignment.py b/tests/backends/test_terracotta_alignment.py index e5767b0..27ffe01 100644 --- a/tests/backends/test_terracotta_alignment.py +++ b/tests/backends/test_terracotta_alignment.py @@ -43,7 +43,7 @@ def _clip_source_to_pseudo_africa(ds): def _sqlite_database_path(path: Path) -> str: - return str(path) + return os.path.relpath(path, Path.cwd()) def _make_test_backend(tmp_path: Path, cache_key: str): From 6634f0a523577ddbe69d7f2e01f6e6693a770771 Mon Sep 17 00:00:00 2001 From: Arthur Jacquiau-Chamski Date: Thu, 2 Apr 2026 11:30:23 +0200 Subject: [PATCH 07/10] refactor/tests: replace resolution-based dataset generation with grid specifications for Terracotta alignment tests --- tests/backends/test_terracotta_alignment.py | 86 ++++++++++----------- 1 file changed, 42 insertions(+), 44 deletions(-) diff --git a/tests/backends/test_terracotta_alignment.py b/tests/backends/test_terracotta_alignment.py index 27ffe01..5b0e949 100644 --- a/tests/backends/test_terracotta_alignment.py +++ b/tests/backends/test_terracotta_alignment.py @@ -14,15 +14,28 @@ from nuthatch.backends.terracotta import TerracottaBackend -GRID_RESOLUTIONS = [1.5, 0.25, 0.1] WEB_MERCATOR = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True) - -def _northup_source_dataset(resolution_degrees): - # Regular lat/lon grid with global longitude coverage and Mercator-safe latitudes. - half_res = resolution_degrees / 2.0 - lons = np.arange(-180.0 + half_res, 180.0, resolution_degrees) - lats = np.arange(84.0 - half_res, -84.0, -resolution_degrees) +# Test grids are regular north-up reference grids. +# Lats are offset by half a grid cell to avoid the poles, which create invalid coordinates in web mercator. +# The 0.1-degree case is kept away from poles to avoid the extreme output-height growth caused by near-pole +# Web Mercator Y values, so that tests are kept lightweight. +GRID_SPECS = { + "global1_5": {"grid_size": 1.5}, + "global0_25": {"grid_size": 0.25}, + "global0_1_avoidpoles": {"grid_size": 0.1, "lat_offset": 5.0}, +} +def _reference_grid_coordinates(grid_name): + grid_spec = GRID_SPECS[grid_name] + grid_size = grid_spec["grid_size"] + lat_offset = grid_spec.get("lat_offset", 0.0) + grid_size/2.0 + lons = np.arange(-180.0, 180.0, grid_size) + lats = np.arange(-90.0 + lat_offset, 90.0 - lat_offset, grid_size) + return lons, lats + + +def _reference_grid_dataset(grid_name): + lons, lats = _reference_grid_coordinates(grid_name) values = np.arange(len(lats) * len(lons), dtype=float).reshape(len(lats), len(lons)) return xr.Dataset( {"precip": (("lat", "lon"), values)}, @@ -77,9 +90,9 @@ def _reference_points_from_dataset(ds): ] -@pytest.mark.parametrize("resolution_degrees", GRID_RESOLUTIONS) -def test_terracotta_write_preserves_latlon_resolution_for_crops(tmp_path, resolution_degrees): - source = _northup_source_dataset(resolution_degrees) +@pytest.mark.parametrize("grid_name", GRID_SPECS) +def test_terracotta_write_preserves_latlon_resolution_for_crops(tmp_path, grid_name): + source = _reference_grid_dataset(grid_name) clipped = _clip_source_to_pseudo_africa(source) global_backend = _make_test_backend(tmp_path, "global") @@ -97,9 +110,9 @@ def test_terracotta_write_preserves_latlon_resolution_for_crops(tmp_path, resolu assert global_tif.crs == africa_tif.crs assert global_tif.res == pytest.approx(africa_tif.res) -@pytest.mark.parametrize("resolution_degrees", GRID_RESOLUTIONS) -def test_terracotta_returns_same_tile_values_for_global_and_africa_scopes(tmp_path, resolution_degrees): - source = _northup_source_dataset(resolution_degrees) +@pytest.mark.parametrize("grid_name", GRID_SPECS) +def test_terracotta_returns_same_tile_values_for_global_and_africa_scopes(tmp_path, grid_name): + source = _reference_grid_dataset(grid_name) clipped = _clip_source_to_pseudo_africa(source) global_backend = _make_test_backend(tmp_path, "global") @@ -125,38 +138,23 @@ def test_terracotta_returns_same_tile_values_for_global_and_africa_scopes(tmp_pa assert np.array_equal(global_tile.filled(-9999), africa_tile.filled(-9999)) -@pytest.mark.parametrize("resolution_degrees", GRID_RESOLUTIONS) -def test_terracotta_write_projects_reference_points_consistently(tmp_path, resolution_degrees): - source = _northup_source_dataset(resolution_degrees) - clipped = _clip_source_to_pseudo_africa(source) +@pytest.mark.parametrize("grid_name", GRID_SPECS) +@pytest.mark.parametrize("scope", ["global", "africa"]) +def test_terracotta_write_projects_reference_points_consistently(tmp_path, grid_name, scope): + source = _reference_grid_dataset(grid_name) + ds_to_write = source if scope == "global" else _clip_source_to_pseudo_africa(source) + backend = _make_test_backend(tmp_path, scope) + backend.write(ds_to_write) + reference_points = _reference_points_from_dataset(ds_to_write) + tif_path = tmp_path / f"{scope}.terracotta" / "_.tif" - global_backend = _make_test_backend(tmp_path, "global") - africa_backend = _make_test_backend(tmp_path, "africa") - global_backend.write(source) - africa_backend.write(clipped) - - reference_points = _reference_points_from_dataset(clipped) - - global_path = tmp_path / "global.terracotta" / "_.tif" - africa_path = tmp_path / "africa.terracotta" / "_.tif" - - with ( - rasterio.open(global_path) as global_tif, - rasterio.open(africa_path) as africa_tif, - ): + with rasterio.open(tif_path) as tif: for lon, lat in reference_points: expected_x, expected_y = WEB_MERCATOR.transform(lon, lat) - observed_centers = [] - - for tif in (global_tif, africa_tif): - row, col = rowcol(tif.transform, expected_x, expected_y) - assert 0 <= row < tif.height - assert 0 <= col < tif.width - - observed_x, observed_y = xy(tif.transform, row, col, offset="center") - observed_centers.append((observed_x, observed_y)) - assert observed_x == pytest.approx(expected_x, abs=tif.res[0] / 2.0) - assert observed_y == pytest.approx(expected_y, abs=tif.res[1] / 2.0) + row, col = rowcol(tif.transform, expected_x, expected_y) + assert 0 <= row < tif.height + assert 0 <= col < tif.width - assert observed_centers[0][0] == pytest.approx(observed_centers[1][0]) - assert observed_centers[0][1] == pytest.approx(observed_centers[1][1]) + observed_x, observed_y = xy(tif.transform, row, col, offset="center") + assert observed_x == pytest.approx(expected_x, abs=tif.res[0] / 2.0) + assert observed_y == pytest.approx(expected_y, abs=tif.res[1] / 2.0) From 0ee413db763d548d879cb014583a1dbc8f865c55 Mon Sep 17 00:00:00 2001 From: Arthur Jacquiau-Chamski Date: Thu, 2 Apr 2026 11:31:21 +0200 Subject: [PATCH 08/10] fix: simplify _get_aligned_mercator_target function by removing clipping to valid mercator --- src/nuthatch/backends/terracotta.py | 41 +---------------------------- 1 file changed, 1 insertion(+), 40 deletions(-) diff --git a/src/nuthatch/backends/terracotta.py b/src/nuthatch/backends/terracotta.py index 6af066e..939c29a 100644 --- a/src/nuthatch/backends/terracotta.py +++ b/src/nuthatch/backends/terracotta.py @@ -140,46 +140,7 @@ def _get_aligned_mercator_target(ds): height, target_resolution, ) - # Derive the valid projected extent directly from the CRS area of use. - area_of_use = CRS.from_epsg(3857).area_of_use - valid_x, valid_y = transform( - "EPSG:4326", - "EPSG:3857", - [area_of_use.west, area_of_use.east], - [area_of_use.south, area_of_use.north], - ) - x_resolution, y_resolution = target_resolution - left = target_transform.c - top = target_transform.f - valid_left, valid_right = min(valid_x), max(valid_x) - valid_bottom, valid_top = min(valid_y), max(valid_y) - right = left + width * x_resolution - bottom = top - height * y_resolution - - # aligned_target can expand the raster past Web Mercator's valid extent, so trim whole rows/columns back. - if left < valid_left: - trim = int(np.ceil((valid_left - left) / x_resolution)) - left += trim * x_resolution - width -= trim - - if right > valid_right: - trim = int(np.ceil((right - valid_right) / x_resolution)) - width -= trim - - if top > valid_top: - trim = int(np.ceil((top - valid_top) / y_resolution)) - top -= trim * y_resolution - height -= trim - - if bottom < valid_bottom: - trim = int(np.ceil((valid_bottom - bottom) / y_resolution)) - height -= trim - - if width <= 0 or height <= 0: - raise ValueError("Aligned raster window fell outside the valid extent for EPSG:3857.") - - # Return a rasterio affine and shape that can be passed straight into rio.reproject(...). - return Affine(x_resolution, 0.0, left, 0.0, -y_resolution, top), width, height + return target_transform, width, height @register_backend From 26edc915a4d5e60a4e8ce6f1a615169883c3a16d Mon Sep 17 00:00:00 2001 From: Arthur Jacquiau-Chamski Date: Sat, 11 Apr 2026 00:15:51 +0200 Subject: [PATCH 09/10] fix: remove pre-projection to Web Mercator and let terracotta project to web mercator on the fly for much better accuracy --- src/nuthatch/backends/terracotta.py | 62 +-------------------- tests/backends/test_terracotta_alignment.py | 21 +++---- 2 files changed, 10 insertions(+), 73 deletions(-) diff --git a/src/nuthatch/backends/terracotta.py b/src/nuthatch/backends/terracotta.py index 939c29a..e443313 100644 --- a/src/nuthatch/backends/terracotta.py +++ b/src/nuthatch/backends/terracotta.py @@ -5,11 +5,8 @@ import xarray as xr import rioxarray # Must import for .rio to work # noqa: F401 import numpy as np -from pyproj import CRS from rasterio.io import MemoryFile from rasterio.enums import Resampling -from rasterio.transform import Affine -from rasterio.warp import aligned_target, calculate_default_transform, transform from nuthatch.backend import DatabaseBackend, FileBackend, register_backend import logging @@ -95,54 +92,6 @@ def lon_base_change(ds, to_base="base180", lon_dim='lon'): ds = ds.sortby('lon') return ds - -def _get_regular_resolution(coords, dim_name): - """Infer the resolution of a regularly spaced coordinate axis.""" - values = np.asarray(coords.values, dtype=float) - if values.size < 2: - raise ValueError(f"Need at least two coordinates along {dim_name} to infer resolution.") - - diffs = np.diff(values) - resolution = float(np.abs(diffs[0])) - if resolution == 0.0: - raise ValueError(f"Coordinate spacing along {dim_name} must be non-zero.") - if not np.allclose(np.abs(diffs), resolution): - raise ValueError(f"Coordinates along {dim_name} must be regularly spaced.") - return resolution - - -def _get_aligned_mercator_target(ds): - """Build a scope-independent Web Mercator target for this raster resolution.""" - # Convert the source cell spacing to an approximate target spacing in Web Mercator meters. - projected_x, projected_y = transform( - "EPSG:4326", - "EPSG:3857", - [0.0, _get_regular_resolution(ds.x, "x")], - [0.0, _get_regular_resolution(ds.y, "y")], - ) - target_resolution = ( - abs(float(projected_x[1] - projected_x[0])), - abs(float(projected_y[1] - projected_y[0])), - ) - - # Ask rasterio for a projected target grid at that resolution, then snap it to aligned pixel edges. - target_transform, width, height = calculate_default_transform( - "EPSG:4326", - "EPSG:3857", - ds.sizes["x"], - ds.sizes["y"], - *ds.rio.bounds(), - resolution=target_resolution, - ) - target_transform, width, height = aligned_target( - target_transform, - width, - height, - target_resolution, - ) - return target_transform, width, height - - @register_backend class TerracottaBackend(DatabaseBackend, FileBackend): """ @@ -240,16 +189,7 @@ def write(self, ds, upsert=False, primary_keys=None): # Adapt the CRS to Web Mercator. ds.rio.write_crs("epsg:4326", inplace=True) - ds.rio.set_spatial_dims("x", "y", inplace=True) - mercator_transform, width, height = _get_aligned_mercator_target(ds) - ds = ds.rio.reproject( - "EPSG:3857", - transform=mercator_transform, - shape=(height, width), - resampling=self.resampling, - nodata=np.nan, - ) - ds.rio.write_crs("epsg:3857", inplace=True) + ds = ds.sortby("y", ascending=False) # Insert the parameters. with self.driver.connect(verify=False): diff --git a/tests/backends/test_terracotta_alignment.py b/tests/backends/test_terracotta_alignment.py index 5b0e949..d5b7978 100644 --- a/tests/backends/test_terracotta_alignment.py +++ b/tests/backends/test_terracotta_alignment.py @@ -6,7 +6,6 @@ import pytest import rasterio import xarray as xr -from pyproj import Transformer from rasterio.transform import rowcol, xy from shapely.geometry import box, mapping from terracotta import get_driver @@ -14,12 +13,7 @@ from nuthatch.backends.terracotta import TerracottaBackend -WEB_MERCATOR = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True) - -# Test grids are regular north-up reference grids. -# Lats are offset by half a grid cell to avoid the poles, which create invalid coordinates in web mercator. -# The 0.1-degree case is kept away from poles to avoid the extreme output-height growth caused by near-pole -# Web Mercator Y values, so that tests are kept lightweight. +# Test grids are regular north-up reference grids. GRID_SPECS = { "global1_5": {"grid_size": 1.5}, "global0_25": {"grid_size": 0.25}, @@ -140,7 +134,9 @@ def test_terracotta_returns_same_tile_values_for_global_and_africa_scopes(tmp_pa @pytest.mark.parametrize("grid_name", GRID_SPECS) @pytest.mark.parametrize("scope", ["global", "africa"]) -def test_terracotta_write_projects_reference_points_consistently(tmp_path, grid_name, scope): +def test_terracotta_write_preserves_reference_points_in_native_crs( + tmp_path, grid_name, scope +): source = _reference_grid_dataset(grid_name) ds_to_write = source if scope == "global" else _clip_source_to_pseudo_africa(source) backend = _make_test_backend(tmp_path, scope) @@ -149,12 +145,13 @@ def test_terracotta_write_projects_reference_points_consistently(tmp_path, grid_ tif_path = tmp_path / f"{scope}.terracotta" / "_.tif" with rasterio.open(tif_path) as tif: + assert tif.crs.to_epsg() == 4326 + for lon, lat in reference_points: - expected_x, expected_y = WEB_MERCATOR.transform(lon, lat) - row, col = rowcol(tif.transform, expected_x, expected_y) + row, col = rowcol(tif.transform, lon, lat) assert 0 <= row < tif.height assert 0 <= col < tif.width observed_x, observed_y = xy(tif.transform, row, col, offset="center") - assert observed_x == pytest.approx(expected_x, abs=tif.res[0] / 2.0) - assert observed_y == pytest.approx(expected_y, abs=tif.res[1] / 2.0) + assert observed_x == pytest.approx(lon, abs=tif.res[0] / 2.0) + assert observed_y == pytest.approx(lat, abs=abs(tif.res[1]) / 2.0) From 55c693a994b30228b3dfed0b8837df7d752b4b02 Mon Sep 17 00:00:00 2001 From: Arthur Jacquiau-Chamski Date: Mon, 20 Apr 2026 12:05:53 +0200 Subject: [PATCH 10/10] comment: fix comments on terracotta write --- src/nuthatch/backends/terracotta.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/nuthatch/backends/terracotta.py b/src/nuthatch/backends/terracotta.py index e443313..95ecb22 100644 --- a/src/nuthatch/backends/terracotta.py +++ b/src/nuthatch/backends/terracotta.py @@ -187,8 +187,10 @@ def write(self, ds, upsert=False, primary_keys=None): else: ds = ds.transpose('y', 'x') - # Adapt the CRS to Web Mercator. + # Label the dataset as WGS84. + # TODO: If the dataset already has a different CRS, we should reproject it to WGS84. ds.rio.write_crs("epsg:4326", inplace=True) + # Write rows north-to-south, as expected by terracotta ds = ds.sortby("y", ascending=False) # Insert the parameters.