Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
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
40 changes: 40 additions & 0 deletions xrspatial/reproject/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,8 @@
'ellipsoidal': 4979, # WGS 84 (3D, ellipsoidal height)
}

_SPATIAL_COORD_RTOL = 1e-6


def _find_spatial_dims(raster):
"""Find the y and x dimension names, handling multi-band rasters.
Expand All @@ -88,6 +90,41 @@ def _find_spatial_dims(raster):
return dims[-2], dims[-1]


def _validate_regular_spatial_coords(raster, context):
"""Validate that spatial coordinates describe a rectilinear grid."""
ydim, xdim = _find_spatial_dims(raster)
_validate_regular_axis(raster.coords[ydim].values, 'y', context)
_validate_regular_axis(raster.coords[xdim].values, 'x', context)


def _validate_regular_axis(values, axis_name, context):
values = np.asarray(values, dtype=np.float64)
if values.size < 2:
return

diffs = np.diff(values)
if not np.all(np.isfinite(diffs)):
raise ValueError(
f"{context}: {axis_name} coordinates must contain finite values"
)

if not (np.all(diffs > 0) or np.all(diffs < 0)):
raise ValueError(
f"{context}: {axis_name} coordinates must be strictly monotonic "
"(all ascending or all descending)"
)

median_step = float(np.median(diffs))
worst_step_deviation = float(np.max(np.abs(diffs - median_step)))
tolerance = _SPATIAL_COORD_RTOL * abs(median_step)
if worst_step_deviation > tolerance:
raise ValueError(
f"{context}: {axis_name} coordinates must be regularly spaced; "
f"worst step deviation is {worst_step_deviation:g} "
f"(median step {median_step:g}, rtol {_SPATIAL_COORD_RTOL:g})"
)


def _source_bounds(raster):
"""Extract (left, bottom, right, top) from a DataArray's coordinates."""
ydim, xdim = _find_spatial_dims(raster)
Expand Down Expand Up @@ -642,6 +679,8 @@ def reproject(

_validate_resampling(resampling)

_validate_regular_spatial_coords(raster, "reproject()")

# Reject unknown vertical-datum tokens at the API boundary so we never
# write None into attrs['vertical_crs'] for typos / unsupported values.
for _name, _val in (('src_vertical_crs', src_vertical_crs),
Expand Down Expand Up @@ -1751,6 +1790,7 @@ def merge(
# (#2027). Reject 3-D up front so callers get a clear error.
_validate_raster(r, func_name='merge', name=f'rasters[{i}]',
ndim=(2,))
_validate_regular_spatial_coords(r, f"merge(): rasters[{i}]")

_validate_grid_params(
resolution=resolution,
Expand Down
41 changes: 41 additions & 0 deletions xrspatial/tests/test_reproject.py
Original file line number Diff line number Diff line change
Expand Up @@ -505,6 +505,33 @@ def test_missing_crs_raises(self):
with pytest.raises(ValueError, match="source CRS"):
reproject(raster, 'EPSG:3857')

def test_reproject_rejects_irregular_x_before_crs_resolution(self):
from xrspatial.reproject import reproject
raster = xr.DataArray(
np.zeros((4, 4)), dims=['y', 'x'],
coords={'y': [3.0, 2.0, 1.0, 0.0],
'x': [0.0, 1.0, 2.3, 3.0]},
)
with pytest.raises(
ValueError,
match=r"reproject\(\): x coordinates.*regular.*worst step deviation",
):
reproject(raster, 'EPSG:4326')

def test_reproject_rejects_non_monotonic_y(self):
from xrspatial.reproject import reproject
raster = xr.DataArray(
np.zeros((4, 4)), dims=['y', 'x'],
coords={'y': [3.0, 2.0, 2.5, 0.0],
'x': [0.0, 1.0, 2.0, 3.0]},
attrs={'crs': 'EPSG:4326'},
)
with pytest.raises(
ValueError,
match=r"reproject\(\): y coordinates must be strictly monotonic",
):
reproject(raster, 'EPSG:4326')

def test_non_dataarray_raises(self):
from xrspatial.reproject import reproject
with pytest.raises(TypeError, match="xarray.DataArray"):
Expand Down Expand Up @@ -653,6 +680,20 @@ def test_merge_invalid_strategy(self):
with pytest.raises(ValueError, match="strategy"):
merge([raster], strategy='median')

def test_merge_rejects_irregular_input_coords(self):
from xrspatial.reproject import merge
raster = xr.DataArray(
np.zeros((4, 4)), dims=['y', 'x'],
coords={'y': [3.0, 2.0, 1.0, 0.0],
'x': [0.0, 1.0, 2.3, 3.0]},
attrs={'crs': 'EPSG:4326'},
)
with pytest.raises(
ValueError,
match=r"merge\(\): rasters\[0\].*x coordinates.*regular.*worst step deviation",
):
merge([raster], target_crs='EPSG:4326')

def test_merge_strategy_last(self):
"""merge() with strategy='last' uses the last valid value."""
from xrspatial.reproject import merge
Expand Down
Loading