resolve test failures

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(

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(

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,
     )

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,
         )