Skip to content

Improve detection of gdal ordering of proj:transform #1570

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
soxofaan wants to merge 2 commits into
base: master
Choose a base branch
from
+89 −25
Open

Improve detection of gdal ordering of proj:transform #1570

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7bcfd14
Improve detection of gdal ordering of proj:transform
soxofaan Mar 2, 2026
9156049
Merge branch 'master' into issue1557-gdal-transpose-detection
soxofaan Mar 4, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
51 changes: 32 additions & 19 deletions openeogeotrellis/load_stac.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1742,32 +1742,45 @@ def get(field):
)

@staticmethod
def _fix_gdal_ordered_transform(transform: Sequence[float]) -> Sequence[float]:
def _fix_gdal_ordered_transform(
transform: Sequence[float], also_check_yx_transposed: bool = False
) -> Sequence[float]:
"""
Detect and fix a proj:transform that is in GDAL GetGeoTransform order
instead of the expected rasterio/affine order.
Detect and fix a proj:transform that is wrongly in GDAL GetGeoTransform order

GDAL GetGeoTransform: [xOrigin, xPixelSize, xSkew, yOrigin, ySkew, yPixelSize]
Rasterio/affine order: [xPixelSize, xSkew, xOrigin, ySkew, yPixelSize, yOrigin]
[xOrigin, xPixelSize, xSkew, yOrigin, ySkew, yPixelSize]

Detection heuristic (for the common no-rotation case):
- In GDAL order: positions 2,4 are zero (skew) and positions 1,5 are non-zero (pixel sizes)
- In rasterio order: positions 1,3 are zero (skew) and positions 0,4 are non-zero (pixel sizes)
Additionally, origin values (large) should be larger than pixel sizes (small).
instead of the expected order (rasterio/affine style):

[xPixelSize, xSkew, xOrigin, ySkew, yPixelSize, yOrigin]

Detection heuristic (for the common XY oriented case):
- In rasterio order: (non-zero) pixel sizes are at positions 0 and 4,
skew values at positions 1 and 3 are small/zero
- In GDAL order: (non-zero) pixel sizes are at positions 1 and 5,
skew values at positions 2 and 4 are small/zero

:param also_check_yx_transposed: whether to also consider the possibility of YX oriented data
"""
if len(transform) < 6:
return transform
t0, t1, t2, t3, t4, t5 = transform[:6]
if (
t2 == 0.0
and t4 == 0.0
and t1 != 0.0
and t5 != 0.0
and abs(t0) > abs(t1)
and abs(t3) > abs(t5)
):
fixed = [t1, t2, t0, t4, t5, t3] + list(transform[6:])
return fixed

# Data consistent with rasterio order?
rasterio_xy_consistent = (abs(t0) > abs(t1)) and (abs(t3) < abs(t4))
rasterio_yx_consistent = also_check_yx_transposed and (abs(t0) < abs(t1)) and (abs(t3) > abs(t4))

# Data consistent with GDAL order?
gdal_xy_consistent = (abs(t1) > abs(t2)) and (abs(t4) < abs(t5))
gdal_yx_consistent = also_check_yx_transposed and (abs(t1) < abs(t2)) and (abs(t4) > abs(t5))

if rasterio_xy_consistent or rasterio_yx_consistent:
pass
elif gdal_xy_consistent or gdal_yx_consistent:
transform = [t1, t2, t0, t4, t5, t3] + list(transform[6:])
else:
logger.debug(f"Failed to detect proj:transform order from {transform=}, leaving as-is.")

return transform

@functools.lru_cache
Expand Down
63 changes: 57 additions & 6 deletions tests/test_load_stac.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1027,16 +1027,67 @@ def test_get_proj_metadata_from_asset(item_properties, asset_extra_fields, expec
class TestFixGdalOrderedTransform:
"""Tests for _ProjectionMetadata._fix_gdal_ordered_transform"""

def test_already_valid_rasterio_order(self):
@pytest.mark.parametrize(
"transform",
[
[0.001, 0.0, 3.0, 0.0, -0.001, 51.0],
[0.001, 0.0, 0.0, 0.0, -0.001, 51.0],
[0.001, 0.0, 0.0, 0.0, -0.001, 0.0],
[0.001, 0.0, 3.0, 0.0, -0.001, 0.0],
],
)
def test_already_valid_rasterio_order(self, transform):
"""Valid rasterio/affine order should not be changed."""
transform = [0.00025, 0.0, 3.0, 0.0, -0.00025, 51.0]
assert _ProjectionMetadata._fix_gdal_ordered_transform(transform) == transform

def test_fix_gdal_order(self):
@pytest.mark.parametrize(
"transform",
[
[0.001, 0.0, 3.0, 0.0, -0.001, 51.0],
[0.001, 0.0, 0.0, 0.0, -0.001, 51.0],
[0.001, 0.0, 0.0, 0.0, -0.001, 0.0],
[0.001, 0.0, 3.0, 0.0, -0.001, 0.0],
[0.0, 0.001, 3.0, -0.001, 0.0, 51.0],
[0.0, 0.001, 0.0, -0.001, 0.0, 51.0],
[0.0, 0.001, 0.0, -0.001, 0.0, 0.0],
[0.0, 0.001, 3.0, -0.001, 0.0, 0.0],
],
)
def test_valid_rasterio_order_including_yx_transpose(self, transform):
"""Valid rasterio/affine order should not be changed."""
assert _ProjectionMetadata._fix_gdal_ordered_transform(transform, also_check_yx_transposed=True) == transform

@pytest.mark.parametrize(
["given", "expected"],
[
([3, 0.001, 0.0, 51, 0.0, -0.001], [0.001, 0.0, 3, 0.0, -0.001, 51]),
([0, 0.001, 0.0, 51, 0.0, -0.001], [0.001, 0.0, 0, 0.0, -0.001, 51]),
([0, 0.001, 0.0, 0, 0.0, -0.001], [0.001, 0.0, 0, 0.0, -0.001, 0]),
([3, 0.001, 0.0, 0, 0.0, -0.001], [0.001, 0.0, 3, 0.0, -0.001, 0]),
],
)
def test_fix_gdal_order(self, given, expected):
"""GDAL GetGeoTransform order should be reshuffled to rasterio/affine order."""
gdal_order = [3.0, 0.00025, 0.0, 51.0, 0.0, -0.00025]
expected = [0.00025, 0.0, 3.0, 0.0, -0.00025, 51.0]
assert _ProjectionMetadata._fix_gdal_ordered_transform(gdal_order) == expected
assert _ProjectionMetadata._fix_gdal_ordered_transform(given) == expected

@pytest.mark.parametrize(
["given", "expected"],
[
([3, 0.001, 0.0, 51, 0.0, -0.001], [0.001, 0.0, 3, 0.0, -0.001, 51]),
# TODO Possible to detect for this case as well? Is consistent with rasterio order
# ([0, 0.001, 0.0, 51, 0.0, -0.001], [0.001, 0.0, 0, 0.0, -0.001, 51]),
([0, 0.001, 0.0, 0, 0.0, -0.001], [0.001, 0.0, 0, 0.0, -0.001, 0]),
([3, 0.001, 0.0, 0, 0.0, -0.001], [0.001, 0.0, 3, 0.0, -0.001, 0]),
([3, 0.0, 0.001, 51, -0.001, 0.0], [0.0, 0.001, 3, -0.001, 0.0, 51]),
([0, 0.0, 0.001, 51, -0.001, 0.0], [0.0, 0.001, 0, -0.001, 0.0, 51]),
([0, 0.0, 0.001, 0, -0.001, 0.0], [0.0, 0.001, 0, -0.001, 0.0, 0]),
# TODO Possible to detect for this case as well? Is consistent with rasterio order
# ([3, 0.0, 0.001, 0, -0.001, 0.0], [0.0, 0.001, 3, -0.001, 0.0, 0]),
],
)
def test_fix_gdal_order_including_yx_transpose(self, given, expected):
"""GDAL GetGeoTransform order should be reshuffled to rasterio/affine order."""
assert _ProjectionMetadata._fix_gdal_ordered_transform(given, also_check_yx_transposed=True) == expected

def test_fix_gdal_order_9_elements(self):
"""GDAL order with 9 elements (full 3x3 matrix) should preserve trailing elements."""
Expand Down