fix: warp rendered GeoTIFFs to Web Mercator so overlays land in place

server/mapcontrol_server/services/geotiff_service.py

@@ -26,8 +26,9 @@
 import httpx
 import numpy as np
 import rasterio
-from rasterio.enums import MaskFlags
-from rasterio.warp import transform_bounds
+from rasterio.crs import CRS
+from rasterio.enums import MaskFlags, Resampling
+from rasterio.warp import calculate_default_transform, reproject, transform_bounds
 from matplotlib import colormaps
 from PIL import Image
 
@@ -186,6 +187,54 @@ def _compute_bounds_4326(dataset) -> list[float]:
     return list(bounds)  # [minLon, minLat, maxLon, maxLat]
 
 
+WEB_MERCATOR = CRS.from_epsg(3857)
+# Web Mercator is undefined past this latitude.
+_MERCATOR_MAX_LAT = 85.05
+
+
+def _to_web_mercator(rgba: np.ndarray, dataset) -> tuple[np.ndarray, list[float]]:
+    """Warp a rendered RGBA image onto a Web Mercator grid.
+
+    The browser places an overlay by its four corners and stretches it
+    linearly in Web Mercator, so only an image on a Web Mercator grid lands
+    where its pixels belong. A UTM or equal-area raster placed by its
+    lon/lat bounding box is rotated and stretched instead: kilometres off
+    for a large raster at high latitude, and two rasters in different
+    projections stop overlapping each other.
+
+    Warps the finished RGBA (nearest neighbour, so mask edges and classes
+    stay crisp; pixels outside the source footprint get alpha 0) and returns
+    it with its exact [minLon, minLat, maxLon, maxLat], an axis-aligned box
+    in both CRSs. Masking already ran on the native grid, unchanged.
+    """
+    bounds_4326 = _compute_bounds_4326(dataset)  # also rejects unplaceable files
+    src_crs = dataset.crs or CRS.from_epsg(4326)
+    if src_crs == WEB_MERCATOR:
+        return rgba, bounds_4326
+    if max(abs(bounds_4326[1]), abs(bounds_4326[3])) > _MERCATOR_MAX_LAT:
+        logger.warning("GeoTIFF reaches past %s° latitude; placing it unwarped", _MERCATOR_MAX_LAT)
+        return rgba, bounds_4326
+
+    height, width = rgba.shape[:2]
+    dst_transform, dst_width, dst_height = calculate_default_transform(
+        src_crs, WEB_MERCATOR, width, height, *dataset.bounds
+    )
+    warped = np.zeros((4, dst_height, dst_width), dtype=np.uint8)
+    reproject(
+        source=np.moveaxis(rgba, -1, 0),
+        destination=warped,
+        src_transform=dataset.transform,
+        src_crs=src_crs,
+        dst_transform=dst_transform,
+        dst_crs=WEB_MERCATOR,
+        resampling=Resampling.nearest,
+    )
+    left, top = dst_transform * (0, 0)
+    right, bottom = dst_transform * (dst_width, dst_height)
+    bounds = list(transform_bounds(WEB_MERCATOR, "EPSG:4326", left, bottom, right, top))
+    return np.moveaxis(warped, 0, -1), bounds
+
+
 def _bounds_to_geojson_polygon(bounds: list[float]) -> dict:
     """Convert [minLon, minLat, maxLon, maxLat] to a GeoJSON Polygon."""
     min_lon, min_lat, max_lon, max_lat = bounds
@@ -445,7 +494,7 @@ async def process_geotiff_rgb(
         # a raster renders.
         def _render() -> GeoTIFFResult:
             rgba, ds = _process_rgb(tif_path, bands, alpha, nodata)
-            bounds = _compute_bounds_4326(ds)
+            rgba, bounds = _to_web_mercator(rgba, ds)
             band_count = ds.count
             crs_str = str(ds.crs) if ds.crs else "unknown"
             ds.close()
@@ -549,7 +598,7 @@ def _render() -> GeoTIFFResult:
                 tif_path, band, colormap, alpha,
                 vmin, vmax, percentile_min, percentile_max, nodata,
             )
-            bounds = _compute_bounds_4326(ds)
+            rgba, bounds = _to_web_mercator(rgba, ds)
             band_count = ds.count
             crs_str = str(ds.crs) if ds.crs else "unknown"
             ds.close()

server/tests/test_geotiff_placement.py

@@ -0,0 +1,130 @@
+"""Where GeoTIFF pixels land on the map: rendered images must be on a Web Mercator grid.
+
+The browser places an overlay by its four lon/lat corners and stretches it
+linearly in Web Mercator. A UTM raster rendered on its native grid and placed
+by its lon/lat bounding box lands rotated and stretched — kilometres off for a
+large raster at high latitude (eo-gpt's Enontekiö runs, 68-69°N).
+"""
+
+from __future__ import annotations
+
+import asyncio
+from types import SimpleNamespace
+
+import numpy as np
+import pytest
+import rasterio
+from PIL import Image
+from rasterio.transform import from_origin
+from rasterio.warp import transform as warp_points
+
+from mapcontrol_server.services import geotiff_service
+
+
+
+class _Points:
+    """lon/lat <-> projected point conversion (scalars or sequences) via rasterio."""
+
+    def __init__(self, src, dst):
+        self.src, self.dst = src, dst
+
+    def transform(self, x, y):
+        scalar = np.isscalar(x)
+        xs, ys = warp_points(self.src, self.dst, np.atleast_1d(x).tolist(), np.atleast_1d(y).tolist())
+        return (xs[0], ys[0]) if scalar else (xs, ys)
+
+
+TO_MERC = _Points("EPSG:4326", "EPSG:3857")
+FROM_MERC = _Points("EPSG:3857", "EPSG:4326")
+TO_UTM = _Points("EPSG:4326", "EPSG:32634")
+
+# Three markers spread over an 80 x 60 km UTM 34N raster at 100 m (Enontekiö).
+MARKERS = {0.4: (21.2, 68.95), 0.7: (22.0, 68.75), 1.0: (22.6, 68.6)}
+RES = 100.0
+
+
+def _utm_raster(path):
+    x0, y0 = TO_UTM.transform(21.0, 69.05)
+    width, height = 800, 600
+    data = np.zeros((1, height, width), dtype=np.float32)
+    for value, (lon, lat) in MARKERS.items():
+        x, y = TO_UTM.transform(lon, lat)
+        data[0, int((y0 - y) // RES), int((x - x0) // RES)] = value
+    with rasterio.open(path, "w", driver="GTiff", count=1, height=height, width=width,
+                       dtype="float32", crs="EPSG:32634", transform=from_origin(x0, y0, RES, RES),
+                       nodata=np.nan) as dst:
+        dst.write(data)
+    return str(path)
+
+
+def _config(tmp_path):
+    return SimpleNamespace(storage=SimpleNamespace(file_dir=str(tmp_path / "files")))
+
+
+def _displayed_lonlat(bounds, shape, row, col):
+    """Where the browser draws PNG pixel (row, col): linear in Web Mercator between the corners."""
+    (mx0, mx1), (my0, my1) = TO_MERC.transform([bounds[0], bounds[2]], [bounds[1], bounds[3]])
+    height, width = shape
+    mx = mx0 + (col + 0.5) / width * (mx1 - mx0)
+    my = my1 - (row + 0.5) / height * (my1 - my0)
+    return FROM_MERC.transform(mx, my)
+
+
+def _ground_m(lon1, lat1, lon2, lat2):
+    dx = (lon2 - lon1) * 111_320 * np.cos(np.radians((lat1 + lat2) / 2))
+    return float(np.hypot(dx, (lat2 - lat1) * 110_570))
+
+
+def test_utm_singleband_markers_land_where_they_belong(tmp_path):
+    result = asyncio.run(geotiff_service.process_geotiff_singleband(
+        _utm_raster(tmp_path / "utm.tif"), "asset-utm", band=1, colormap="gray",
+        vmin=0.0, vmax=1.0, config=_config(tmp_path)))
+    assert isinstance(result, geotiff_service.GeoTIFFResult), result
+    rgba = np.asarray(Image.open(result.png_path))
+    gray = rgba[:, :, 0].astype(int)
+    for value, (lon, lat) in MARKERS.items():
+        rows, cols = np.nonzero((np.abs(gray - round(value * 255)) <= 2) & (rgba[:, :, 3] == 255))
+        assert len(rows), f"marker {value} missing"
+        shown = _displayed_lonlat(result.bounds, gray.shape, rows.mean(), cols.mean())
+        # Within ~2 source pixels of the truth (was 1-4 km off on the native grid).
+        assert _ground_m(*shown, lon, lat) < 2.5 * RES, (value, shown, (lon, lat))
+
+
+def test_native_grid_placement_would_have_been_kilometres_off(tmp_path):
+    """The failure mode this guards against, measured on the same raster."""
+    path = _utm_raster(tmp_path / "utm.tif")
+    with rasterio.open(path) as ds:
+        bounds = geotiff_service._compute_bounds_4326(ds)
+        data = ds.read(1)
+    errors = []
+    for value, (lon, lat) in MARKERS.items():
+        rows, cols = np.nonzero(data == np.float32(value))
+        shown = _displayed_lonlat(bounds, data.shape, rows.mean(), cols.mean())
+        errors.append(_ground_m(*shown, lon, lat))
+    assert max(errors) > 1000, errors
+
+
+def test_rgb_footprint_outside_is_transparent_and_inside_opaque(tmp_path):
+    x0, y0 = TO_UTM.transform(21.0, 69.05)
+    path = tmp_path / "rgb.tif"
+    with rasterio.open(path, "w", driver="GTiff", count=3, height=300, width=400, dtype="uint8",
+                       crs="EPSG:32634", transform=from_origin(x0, y0, RES, RES)) as dst:
+        dst.write(np.full((3, 300, 400), 120, dtype=np.uint8))
+    result = asyncio.run(geotiff_service.process_geotiff_rgb(
+        str(path), "asset-rgb", config=_config(tmp_path)))
+    alpha = np.asarray(Image.open(result.png_path))[:, :, 3]
+    # The rotated UTM footprint leaves at least one corner of the Mercator box empty.
+    assert min(alpha[0, 0], alpha[0, -1], alpha[-1, 0], alpha[-1, -1]) == 0
+    assert alpha[alpha.shape[0] // 2, alpha.shape[1] // 2] == 255
+
+
+def test_web_mercator_input_is_not_resampled(tmp_path):
+    path = tmp_path / "merc.tif"
+    mx, my = TO_MERC.transform(10.0, 50.0)
+    with rasterio.open(path, "w", driver="GTiff", count=1, height=40, width=50, dtype="float32",
+                       crs="EPSG:3857", transform=from_origin(mx, my, 250, 250)) as dst:
+        dst.write(np.random.default_rng(0).random((1, 40, 50), dtype=np.float32))
+    result = asyncio.run(geotiff_service.process_geotiff_singleband(
+        str(path), "asset-merc", config=_config(tmp_path)))
+    assert (result.width, result.height) == (50, 40)
+    assert result.bounds[0] == pytest.approx(10.0) and result.bounds[3] == pytest.approx(50.0)