diff --git a/CHANGES.md b/CHANGES.md index 7265944..76b7517 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -10,6 +10,8 @@ `xcube_resampling.gridmapping.GridMapping`. These methods were only used by `xcube_resampling.rectify_dataset`; their updated implementations now live in `xcube_resampling.rectify._compute_source_tile_indexing`. +- Added support for rectifying datasets with decreasing x-coordinate in + `xcube_resampling.rectify_dataset`. ## Changes in 0.2.4 diff --git a/tests/test_rectify.py b/tests/test_rectify.py index 89926c2..1d41b82 100644 --- a/tests/test_rectify.py +++ b/tests/test_rectify.py @@ -89,6 +89,33 @@ def test_rectify_2x2_to_default(self): ), ) + def test_rectify_2x2_decreasing_x_coords(self): + source_ds = create_2x2_dataset_with_irregular_coords() + source_ds = source_ds.isel({"x": slice(None, None, -1)}) + + target_gm = GridMapping.regular( + size=(4, 4), xy_min=(0, 50), xy_res=2, crs=CRS_WGS84 + ) + target_ds = rectify_dataset( + source_ds, + target_gm=target_gm, + interp_methods=0, + output_indices_names=("indices_lon", "indices_lat"), + ) + + np.testing.assert_almost_equal( + target_ds.rad.values, + np.array( + [ + [nan, nan, nan, nan], + [nan, 1.0, 2.0, nan], + [3.0, 3.0, 2.0, nan], + [nan, 4.0, nan, nan], + ], + dtype=target_ds.rad.dtype, + ), + ) + def test_rectify_2x2_to_regular(self): source_ds = create_2x2_dataset_with_irregular_coords() target_ds = rectify_dataset(source_ds, interp_methods=0) diff --git a/xcube_resampling/rectify.py b/xcube_resampling/rectify.py index d873d3c..d3b40ad 100644 --- a/xcube_resampling/rectify.py +++ b/xcube_resampling/rectify.py @@ -39,6 +39,8 @@ SpatialInterpMethodStr, ) from .gridmapping import GridMapping +from .gridmapping.cfconv import get_dataset_grid_mapping_proxies +from .gridmapping.helpers import from_lon_360 from .utils import ( SourceTileIndexing, _clip_if_needed, @@ -128,6 +130,7 @@ def rectify_dataset( coordinate variables are ignored in the output. """ source_ds = _select_variables(source_ds, variables) + source_ds = _ensure_increasing_x(source_ds) source_gm = source_gm or GridMapping.from_dataset(source_ds) source_ds = normalize_grid_mapping(source_ds, source_gm) @@ -200,9 +203,11 @@ def _assemble_rectified_dataset( coords = source_ds.coords.to_dataset() coords = coords.drop_vars((src_x_name, src_y_name), errors="ignore") - target_coords = target_gm.to_coords() - coords[tgt_x_name] = target_coords[tgt_x_name] - coords[tgt_y_name] = target_coords[tgt_y_name] + tgt_x = target_gm.x_coords + if target_gm.is_lon_360: + tgt_x = from_lon_360(tgt_x) + coords[tgt_x_name] = tgt_x + coords[tgt_y_name] = target_gm.y_coords coords["spatial_ref"] = xr.DataArray(0, attrs=target_gm.crs.to_cf()) target_ds = xr.Dataset(coords=coords, attrs=source_ds.attrs) @@ -721,3 +726,13 @@ def _rectify_block( sampled = _sample_array_at_indices(data_array, ix[valid], iy[valid], interp_method) data_rectified[:, valid] = sampled return data_rectified + + +def _ensure_increasing_x(ds: xr.Dataset) -> xr.Dataset: + gm_proxy = next(iter(get_dataset_grid_mapping_proxies(ds).values())) + x_coords = gm_proxy.coords.x + x_diff = x_coords[0, 1] - x_coords[0, 0] + if -180 < x_diff < 0: + x_dim = x_coords.dims[1] + ds = ds.isel({x_dim: slice(None, None, -1)}) + return ds diff --git a/xcube_resampling/reproject.py b/xcube_resampling/reproject.py index 3bdcd69..2c09937 100644 --- a/xcube_resampling/reproject.py +++ b/xcube_resampling/reproject.py @@ -37,6 +37,7 @@ SpatialInterpMethodStr, ) from .gridmapping import GridMapping +from .gridmapping.helpers import from_lon_360 from .utils import ( SourceTileIndexing, _clip_if_needed, @@ -187,12 +188,17 @@ def _assemble_reprojected_dataset( interp_methods: SpatialInterpMethods | None = None, fill_values: FillValues | None = None, ): - x_name, y_name = source_gm.xy_var_names + src_x_name, src_y_name = source_gm.xy_var_names + tgt_x_name, tgt_y_name = target_gm.xy_var_names + coords = source_ds.coords.to_dataset() - coords = coords.drop_vars((x_name, y_name), errors="ignore") - x_name, y_name = target_gm.xy_var_names - coords[x_name] = target_gm.x_coords - coords[y_name] = target_gm.y_coords + coords = coords.drop_vars((src_x_name, src_y_name), errors="ignore") + + tgt_x = target_gm.x_coords + if target_gm.is_lon_360: + tgt_x = from_lon_360(tgt_x) + coords[tgt_x_name] = tgt_x + coords[tgt_y_name] = target_gm.y_coords coords["spatial_ref"] = xr.DataArray(0, attrs=target_gm.crs.to_cf()) target_ds = xr.Dataset(coords=coords, attrs=source_ds.attrs)