diff --git a/test/core/resampling/test_rectify.py b/test/core/resampling/test_rectify.py index dcc097455..7f5093684 100644 --- a/test/core/resampling/test_rectify.py +++ b/test/core/resampling/test_rectify.py @@ -868,7 +868,8 @@ def _assert_compute_and_extract_source_pixels( target_rad = np.full((13, 13), np.nan, dtype=np.float64) - compute_var_image(source_ds.rad.values, dst_src_ij, target_rad, 0) + src_bbox = [0, 0, source_ds.rad.shape[-1], source_ds.rad.shape[-2]] + compute_var_image(source_ds.rad.values, dst_src_ij, target_rad, src_bbox, 0) if not is_j_axis_up: np.testing.assert_almost_equal( diff --git a/xcube/core/gridmapping/coords.py b/xcube/core/gridmapping/coords.py index 84a9b287c..0eaa44d5b 100644 --- a/xcube/core/gridmapping/coords.py +++ b/xcube/core/gridmapping/coords.py @@ -4,7 +4,7 @@ import abc import math -from typing import Tuple, Union, Dict +from typing import Union import dask.array as da import numpy as np @@ -85,6 +85,7 @@ def new_grid_mapping_from_coords( crs: Union[str, pyproj.crs.CRS], *, xy_res: Union[Number, tuple[Number, Number]] = None, + xy_bbox: tuple[Number, Number, Number, Number] = None, tile_size: Union[int, tuple[int, int]] = None, tolerance: float = DEFAULT_TOLERANCE, ) -> GridMapping: @@ -139,8 +140,8 @@ def new_grid_mapping_from_coords( x_diff, y_res, atol=tolerance ) if is_regular: - x_res = round_to_fraction(x_res, 5, 0.25) - y_res = round_to_fraction(y_res, 5, 0.25) + x_res = round_to_fraction(float(x_res), 5, 0.25) + y_res = round_to_fraction(float(y_res), 5, 0.25) else: x_res = round_to_fraction(float(np.nanmedian(x_diff)), 2, 0.5) y_res = round_to_fraction(float(np.nanmedian(y_diff)), 2, 0.5) @@ -176,10 +177,10 @@ def new_grid_mapping_from_coords( x = da.asarray(x_coords) y = da.asarray(y_coords) - x_x_diff = _abs_no_nan(da.diff(x[0, :])) - x_y_diff = _abs_no_nan(da.diff(x[:, 0])) - y_x_diff = _abs_no_nan(da.diff(y[0, :])) - y_y_diff = _abs_no_nan(da.diff(y[:, 0])) + x_x_diff = _abs_no_nan(da.diff(x[0, : x.chunksize[1]])) + x_y_diff = _abs_no_nan(da.diff(x[: x.chunksize[0], 0])) + y_x_diff = _abs_no_nan(da.diff(y[0, : x.chunksize[0]])) + y_y_diff = _abs_no_nan(da.diff(y[: x.chunksize[1], 0])) if not is_lon_360 and crs.is_geographic: is_anti_meridian_crossed = da.any(da.max(x_x_diff) > 180) or da.any( @@ -201,8 +202,8 @@ def new_grid_mapping_from_coords( is_regular = ( da.allclose(x_x_diff, x_res, atol=tolerance) and da.allclose(y_y_diff, y_res, atol=tolerance) - and da.all(x_y_diff == 0) - and da.all(y_x_diff == 0) + and da.allclose(x_y_diff, 0, atol=tolerance) + and da.allclose(y_x_diff, 0, atol=tolerance) ) if not is_regular and xy_res is None: @@ -276,7 +277,9 @@ def new_grid_mapping_from_coords( ) # Guess j axis direction - is_j_axis_up = da.all(y_coords[0, :] < y_coords[-1, :]) or None + is_j_axis_up = da.all( + y_coords[0, : y.chunksize[1]] < y_coords[-1, : y.chunksize[1]] + ) assert_true( x_res > 0 and y_res > 0, @@ -285,11 +288,17 @@ def new_grid_mapping_from_coords( ) x_res, y_res = _to_int_or_float(x_res), _to_int_or_float(y_res) - x_res_05, y_res_05 = x_res / 2, y_res / 2 - x_min = _to_int_or_float(x_coords.min() - x_res_05) - y_min = _to_int_or_float(y_coords.min() - y_res_05) - x_max = _to_int_or_float(x_coords.max() + x_res_05) - y_max = _to_int_or_float(y_coords.max() + y_res_05) + if xy_bbox is None: + x_res_05, y_res_05 = x_res / 2, y_res / 2 + x_min = _to_int_or_float(x_coords[..., 0].min() - x_res_05) + x_max = _to_int_or_float(x_coords[..., -1].max() + x_res_05) + if is_j_axis_up: + y_min = _to_int_or_float(float(y_coords[0, ...].min()) - y_res_05) + y_max = _to_int_or_float(float(y_coords[-1, ...].max()) + y_res_05) + else: + y_min = _to_int_or_float(float(y_coords[-1, ...].min()) - y_res_05) + y_max = _to_int_or_float(float(y_coords[0, ...].max()) + y_res_05) + xy_bbox = (x_min, y_min, x_max, y_max) if cls is Coords1DGridMapping and is_regular: from .regular import RegularGridMapping @@ -302,7 +311,7 @@ def new_grid_mapping_from_coords( crs=crs, size=size, tile_size=tile_size, - xy_bbox=(x_min, y_min, x_max, y_max), + xy_bbox=xy_bbox, xy_res=(x_res, y_res), xy_var_names=xy_var_names, xy_dim_names=(str(x_dim), str(y_dim)), @@ -313,13 +322,13 @@ def new_grid_mapping_from_coords( def _abs_no_zero(array: Union[xr.DataArray, da.Array, np.ndarray]): - array = np.fabs(array) - return np.where(np.isclose(array, 0), np.nan, array) + array = da.fabs(array) + return da.where(da.isclose(array, 0), np.nan, array) -def _abs_no_nan(array: Union[xr.DataArray, da.Array, np.ndarray]): - array = np.fabs(array) - return np.where(np.logical_or(np.isnan(array), np.isclose(array, 0)), 0, array) +def _abs_no_nan(array: Union[da.Array, np.ndarray]): + array = da.fabs(array) + return da.where(da.logical_or(da.isnan(array), da.isclose(array, 0)), 0, array) def grid_mapping_to_coords( diff --git a/xcube/core/gridmapping/transform.py b/xcube/core/gridmapping/transform.py index 11eea404a..39bff2456 100644 --- a/xcube/core/gridmapping/transform.py +++ b/xcube/core/gridmapping/transform.py @@ -2,7 +2,7 @@ # Permissions are hereby granted under the terms of the MIT License: # https://opensource.org/licenses/MIT. -from typing import Union, Tuple +from typing import Union import numpy as np import pyproj @@ -70,6 +70,9 @@ def _transform(block: np.ndarray) -> np.ndarray: output_dtypes=[np.float64], dask="parallelized", ) + xy_min = transformer.transform(*grid_mapping.xy_bbox[:2]) + xy_max = transformer.transform(*grid_mapping.xy_bbox[2:]) + xy_bbox = xy_min + xy_max xy_var_names = xy_var_names or ("transformed_x", "transformed_y") @@ -90,6 +93,7 @@ def _transform(block: np.ndarray) -> np.ndarray: y_coords=xy_coords[1].rename(xy_var_names[1]), crs=target_crs, xy_res=xy_res, + xy_bbox=xy_bbox, tile_size=tile_size, tolerance=tolerance, ) diff --git a/xcube/core/resampling/rectify.py b/xcube/core/resampling/rectify.py index f94ec585e..df58fe555 100644 --- a/xcube/core/resampling/rectify.py +++ b/xcube/core/resampling/rectify.py @@ -642,8 +642,9 @@ def _compute_var_image_numpy( dst_height = dst_src_ij_images.shape[-2] dst_shape = src_var.shape[:-2] + (dst_height, dst_width) dst_values = np.full(dst_shape, fill_value, dtype=src_var.dtype) + src_bbox = [0, 0, src_var.shape[-2], src_var.shape[-1]] _compute_var_image_numpy_parallel( - src_var, dst_src_ij_images, dst_values, interpolation + src_var, dst_src_ij_images, dst_values, src_bbox, interpolation ) return dst_values @@ -686,6 +687,7 @@ def _compute_var_image_numpy_parallel( src_var_image: np.ndarray, dst_src_ij_images: np.ndarray, dst_var_image: np.ndarray, + src_bbox: list[int], interpolation: int, ): """Extract source pixels from np.ndarray source @@ -698,6 +700,7 @@ def _compute_var_image_numpy_parallel( src_var_image, dst_src_ij_images, dst_var_image, + src_bbox, interpolation, )