to_geotiff: round-trip descending x and ascending y (#1716)

The writer stored abs(pixel_width)/abs(pixel_height) in ModelPixelScaleTag,
and the reader hard-coded a north-up reconstruction. DataArrays with
descending x or ascending y silently round-tripped with the wrong sign.

Emit ModelTransformationTag (34264) for non-standard orientations: that
tag carries the full signed 4x4 affine, which ModelPixelScale cannot
because the spec requires positive scales. The reader already preferred
ModelTransformationTag when present; only the writer side needed work.

Also fix write_vrt's dst_y_off / dst_x_off math to use the top-left
corner of each source's geographic extent rather than blindly assuming
origin_y is the top. Without this, mosaics built from sources with
ascending y placed tiles outside the VRT extent.

Author: brendancol

xrspatial/geotiff/_geotags.py

@@ -888,16 +888,34 @@ def build_geo_tags(transform: GeoTransform, crs_epsg: int | None = None,
     """
     tags = {}
 
-    # ModelPixelScaleTag (33550): (ScaleX, ScaleY, ScaleZ)
-    sx = abs(transform.pixel_width)
-    sy = abs(transform.pixel_height)
-    tags[TAG_MODEL_PIXEL_SCALE] = (sx, sy, 0.0)
-
-    # ModelTiepointTag (33922): (I, J, K, X, Y, Z)
-    tags[TAG_MODEL_TIEPOINT] = (
-        0.0, 0.0, 0.0,
-        transform.origin_x, transform.origin_y, 0.0,
-    )
+    # Standard north-up rasters have pixel_width > 0 and pixel_height < 0.
+    # Anything else (descending x, ascending y) cannot be expressed via
+    # ModelPixelScale + ModelTiepoint because the spec requires the scales
+    # to be positive.
+    north_up = transform.pixel_width > 0 and transform.pixel_height < 0
+
+    if north_up:
+        # ModelPixelScaleTag (33550): (ScaleX, ScaleY, ScaleZ)
+        sx = transform.pixel_width
+        sy = -transform.pixel_height
+        tags[TAG_MODEL_PIXEL_SCALE] = (sx, sy, 0.0)
+
+        # ModelTiepointTag (33922): (I, J, K, X, Y, Z)
+        tags[TAG_MODEL_TIEPOINT] = (
+            0.0, 0.0, 0.0,
+            transform.origin_x, transform.origin_y, 0.0,
+        )
+    else:
+        # Why: GeoTIFF spec requires ModelPixelScale entries to be positive,
+        # so non-standard orientations (descending x, ascending y) must use
+        # ModelTransformationTag (34264) instead.  The 4x4 row-major matrix
+        # maps (col, row, 0, 1) -> (X, Y, 0, 1).
+        tags[TAG_MODEL_TRANSFORMATION] = (
+            transform.pixel_width, 0.0, 0.0, transform.origin_x,
+            0.0, transform.pixel_height, 0.0, transform.origin_y,
+            0.0, 0.0, 0.0, 0.0,
+            0.0, 0.0, 0.0, 1.0,
+        )
 
     # GeoKeyDirectoryTag (34735)
     geokeys = []

xrspatial/geotiff/_vrt.py

@@ -908,9 +908,21 @@ def write_vrt(vrt_path: str, source_files: list[str], *,
 
         for m in sources_meta:
             t = m['transform']
+            # Top edge of this source in geo coords -- origin_y is the
+            # *top* only for north-up rasters (pixel_height < 0).  Sources
+            # with ascending y (pixel_height > 0) place origin_y at the
+            # bottom, so the top is origin_y + height * pixel_height.
+            src_top = max(
+                t.origin_y,
+                t.origin_y + m['height'] * t.pixel_height,
+            )
+            src_left = min(
+                t.origin_x,
+                t.origin_x + m['width'] * t.pixel_width,
+            )
             # Pixel offset in the virtual raster
-            dst_x_off = int(round((t.origin_x - mosaic_x0) / abs(res_x)))
-            dst_y_off = int(round((mosaic_y_top - t.origin_y) / abs(res_y)))
+            dst_x_off = int(round((src_left - mosaic_x0) / abs(res_x)))
+            dst_y_off = int(round((mosaic_y_top - src_top) / abs(res_y)))
 
             fname = m['path']
             rel_attr = '0'

xrspatial/geotiff/_writer.py

@@ -40,6 +40,7 @@
     TAG_GDAL_NODATA,
     TAG_MODEL_PIXEL_SCALE,
     TAG_MODEL_TIEPOINT,
+    TAG_MODEL_TRANSFORMATION,
 )
 from ._header import (
     TAG_NEW_SUBFILE_TYPE,
@@ -873,6 +874,8 @@ def _assemble_tiff(width: int, height: int, dtype: np.dtype,
                     tags.append((gtag, DOUBLE, 3, list(gval)))
                 elif gtag == TAG_MODEL_TIEPOINT:
                     tags.append((gtag, DOUBLE, 6, list(gval)))
+                elif gtag == TAG_MODEL_TRANSFORMATION:
+                    tags.append((gtag, DOUBLE, 16, list(gval)))
                 elif gtag == TAG_GEO_KEY_DIRECTORY:
                     tags.append((gtag, SHORT, len(gval), list(gval)))
                 elif gtag == TAG_GEO_ASCII_PARAMS:
@@ -1516,6 +1519,8 @@ def write_streaming(dask_data, path: str, *,
             tags.append((gtag, DOUBLE, 3, list(gval)))
         elif gtag == TAG_MODEL_TIEPOINT:
             tags.append((gtag, DOUBLE, 6, list(gval)))
+        elif gtag == TAG_MODEL_TRANSFORMATION:
+            tags.append((gtag, DOUBLE, 16, list(gval)))
         elif gtag == TAG_GEO_KEY_DIRECTORY:
             tags.append((gtag, SHORT, len(gval), list(gval)))
         elif gtag == TAG_GEO_ASCII_PARAMS:

xrspatial/geotiff/tests/test_descending_coords_1716.py

@@ -0,0 +1,133 @@
+"""Regression tests for issue #1716.
+
+``to_geotiff`` previously stored ``abs(pixel_width)`` / ``abs(pixel_height)``
+in ModelPixelScaleTag and the reader hard-coded a north-up reconstruction.
+DataArrays with descending x or ascending y silently round-tripped with the
+wrong georeference.  The writer now emits ModelTransformationTag (34264)
+for non-standard orientations so the sign survives the round trip.
+"""
+from __future__ import annotations
+
+import os
+
+import numpy as np
+import pytest
+import xarray as xr
+
+from xrspatial.geotiff import open_geotiff, to_geotiff
+from xrspatial.geotiff._geotags import (
+    TAG_MODEL_PIXEL_SCALE,
+    TAG_MODEL_TIEPOINT,
+    TAG_MODEL_TRANSFORMATION,
+)
+from xrspatial.geotiff._header import parse_all_ifds, parse_header
+
+
+def _ifd_tag_ids(path: str) -> set[int]:
+    with open(path, 'rb') as fh:
+        data = fh.read()
+    header = parse_header(data)
+    ifds = parse_all_ifds(data, header)
+    return set(ifds[0].entries.keys())
+
+
+def _make_da(x_coords: np.ndarray, y_coords: np.ndarray) -> xr.DataArray:
+    arr = np.arange(len(y_coords) * len(x_coords), dtype=np.float32)
+    arr = arr.reshape(len(y_coords), len(x_coords))
+    return xr.DataArray(
+        arr,
+        dims=('y', 'x'),
+        coords={'y': y_coords, 'x': x_coords},
+    )
+
+
+def test_descending_x_roundtrip(tmp_path):
+    """Descending x coords survive the round trip."""
+    # x decreases left-to-right (unusual but valid)
+    x = np.array([200.0, 190.0, 180.0, 170.0, 160.0], dtype=np.float64)
+    y = np.array([50.0, 40.0, 30.0, 20.0], dtype=np.float64)  # north-up
+    da = _make_da(x, y)
+
+    out = tmp_path / 'tmp_1716_desc_x.tif'
+    to_geotiff(da, str(out), crs=4326)
+
+    loaded = open_geotiff(str(out))
+    np.testing.assert_allclose(loaded.coords['x'].values, x)
+    np.testing.assert_allclose(loaded.coords['y'].values, y)
+    np.testing.assert_array_equal(loaded.values, da.values)
+
+
+def test_ascending_y_roundtrip(tmp_path):
+    """Ascending y coords survive the round trip."""
+    x = np.array([160.0, 170.0, 180.0, 190.0, 200.0], dtype=np.float64)
+    # y increases top-to-bottom (south-up)
+    y = np.array([20.0, 30.0, 40.0, 50.0], dtype=np.float64)
+    da = _make_da(x, y)
+
+    out = tmp_path / 'tmp_1716_asc_y.tif'
+    to_geotiff(da, str(out), crs=4326)
+
+    loaded = open_geotiff(str(out))
+    np.testing.assert_allclose(loaded.coords['x'].values, x)
+    np.testing.assert_allclose(loaded.coords['y'].values, y)
+    np.testing.assert_array_equal(loaded.values, da.values)
+
+
+def test_descending_x_and_ascending_y_roundtrip(tmp_path):
+    """Both axes flipped relative to north-up."""
+    x = np.array([200.0, 190.0, 180.0, 170.0, 160.0], dtype=np.float64)
+    y = np.array([20.0, 30.0, 40.0, 50.0], dtype=np.float64)
+    da = _make_da(x, y)
+
+    out = tmp_path / 'tmp_1716_desc_x_asc_y.tif'
+    to_geotiff(da, str(out), crs=4326)
+
+    loaded = open_geotiff(str(out))
+    np.testing.assert_allclose(loaded.coords['x'].values, x)
+    np.testing.assert_allclose(loaded.coords['y'].values, y)
+    np.testing.assert_array_equal(loaded.values, da.values)
+
+
+def test_north_up_still_uses_pixel_scale_and_tiepoint(tmp_path):
+    """Standard north-up orientation keeps ModelPixelScale + ModelTiepoint."""
+    x = np.array([160.0, 170.0, 180.0, 190.0, 200.0], dtype=np.float64)
+    y = np.array([50.0, 40.0, 30.0, 20.0], dtype=np.float64)
+    da = _make_da(x, y)
+
+    out = tmp_path / 'tmp_1716_north_up.tif'
+    to_geotiff(da, str(out), crs=4326)
+
+    tag_ids = _ifd_tag_ids(str(out))
+    assert TAG_MODEL_PIXEL_SCALE in tag_ids
+    assert TAG_MODEL_TIEPOINT in tag_ids
+    assert TAG_MODEL_TRANSFORMATION not in tag_ids
+
+
+def test_descending_x_uses_transformation_tag(tmp_path):
+    """Non-standard orientation emits ModelTransformationTag and skips
+    the scale/tiepoint pair."""
+    x = np.array([200.0, 190.0, 180.0, 170.0, 160.0], dtype=np.float64)
+    y = np.array([50.0, 40.0, 30.0, 20.0], dtype=np.float64)
+    da = _make_da(x, y)
+
+    out = tmp_path / 'tmp_1716_desc_x_tags.tif'
+    to_geotiff(da, str(out), crs=4326)
+
+    tag_ids = _ifd_tag_ids(str(out))
+    assert TAG_MODEL_TRANSFORMATION in tag_ids
+    assert TAG_MODEL_PIXEL_SCALE not in tag_ids
+    assert TAG_MODEL_TIEPOINT not in tag_ids
+
+
+def test_ascending_y_uses_transformation_tag(tmp_path):
+    x = np.array([160.0, 170.0, 180.0, 190.0, 200.0], dtype=np.float64)
+    y = np.array([20.0, 30.0, 40.0, 50.0], dtype=np.float64)
+    da = _make_da(x, y)
+
+    out = tmp_path / 'tmp_1716_asc_y_tags.tif'
+    to_geotiff(da, str(out), crs=4326)
+
+    tag_ids = _ifd_tag_ids(str(out))
+    assert TAG_MODEL_TRANSFORMATION in tag_ids
+    assert TAG_MODEL_PIXEL_SCALE not in tag_ids
+    assert TAG_MODEL_TIEPOINT not in tag_ids