Skip to content
Merged
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
2 changes: 2 additions & 0 deletions CHANGES.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
27 changes: 27 additions & 0 deletions tests/test_rectify.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
21 changes: 18 additions & 3 deletions xcube_resampling/rectify.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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
16 changes: 11 additions & 5 deletions xcube_resampling/reproject.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@
SpatialInterpMethodStr,
)
from .gridmapping import GridMapping
from .gridmapping.helpers import from_lon_360
from .utils import (
SourceTileIndexing,
_clip_if_needed,
Expand Down Expand Up @@ -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)

Expand Down