diff --git a/xrspatial/geotiff/_geotags.py b/xrspatial/geotiff/_geotags.py index ff94dcef..972817e1 100644 --- a/xrspatial/geotiff/_geotags.py +++ b/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 = [] diff --git a/xrspatial/geotiff/_vrt.py b/xrspatial/geotiff/_vrt.py index 6cc002a0..090cff0b 100644 --- a/xrspatial/geotiff/_vrt.py +++ b/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' diff --git a/xrspatial/geotiff/_writer.py b/xrspatial/geotiff/_writer.py index 7758adcb..f3423114 100644 --- a/xrspatial/geotiff/_writer.py +++ b/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: diff --git a/xrspatial/geotiff/tests/test_descending_coords_1716.py b/xrspatial/geotiff/tests/test_descending_coords_1716.py new file mode 100644 index 00000000..3903e487 --- /dev/null +++ b/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