ready for review

Author: konstntokas

CHANGES.md

@@ -1,3 +1,12 @@
+## Changes in 0.2.0 (under development)
+
+- Enhanced the function `bbox_overlap` so that it can handle bounding boxed crossing
+  the antimeridian
+- Added new class method `GridMapping.regular_from_bbox`, which allows creating a
+  regular grid mapping directly from a bounding box, spatial resolution, and CRS.
+- Bug fix: fixed grid mapping creation for irregular grids with decreasing longitude
+  along axis 1.
+
 ## Changes in 0.1.1
 
 - Improved `xcube_resampling.utils.clip_dataset_by_bbox` to support datasets with 

tests/gridmapping/test_regular.py

@@ -335,3 +335,13 @@ def test_to_regular(self):
         self.assertEqual(gm_test.crs, CRS_WGS84)
         self.assertEqual(gm_test.xy_res, (0.01, 0.01))
         self.assertTrue(gm_test.is_j_axis_up)
+
+    def test_create_regular_from_bbox(self):
+        gm_test = GridMapping.regular_from_bbox(
+            [10, 40, 20, 50], 0.1, "EPSG:4326", tile_size=10
+        )
+        self.assertEqual(gm_test.size, (100, 100))
+        self.assertEqual(gm_test.tile_size, (10, 10))
+        self.assertEqual(gm_test.crs, CRS_WGS84)
+        self.assertEqual(gm_test.xy_res, (0.1, 0.1))
+        self.assertFalse(gm_test.is_j_axis_up)

tests/test_utils.py

@@ -330,6 +330,33 @@ def test_bbox_overlap(self):
         expected = 36 / 36
         self.assertAlmostEqual(bbox_overlap(source, target), expected)
 
+        # antimeridian - identical wrapped boxes (e.g. 170°E → -170°W)
+        source = (170, -10, -170, 10)
+        target = (170, -10, -170, 10)
+        self.assertEqual(bbox_overlap(source, target), 1.0)
+
+        # antimeridian - source crosses, target fully inside one side
+        source = (170, 0, -170, 10)
+        target = (175, 0, 178, 10)
+        expected = 30 / 200
+        self.assertAlmostEqual(bbox_overlap(source, target), expected)
+
+        # antimeridian - partial overlap across both segments
+        source = (170, 0, -170, 10)
+        target = (160, 0, 175, 10)
+        expected = 50 / 200
+        self.assertAlmostEqual(bbox_overlap(source, target), expected)
+
+        # antimeridian - non-wrapped target covering entire world
+        source = (170, 0, -170, 10)
+        target = (-180, -90, 180, 90)
+        self.assertEqual(bbox_overlap(source, target), 1.0)
+
+        # antimeridian - touching at boundary only (no real overlap)
+        source = (170, 0, -170, 10)
+        target = (-170, 0, -160, 10)
+        self.assertAlmostEqual(bbox_overlap(source, target), 0.0)
+
 
 class TestClipDatasetByBBox(unittest.TestCase):
 

xcube_resampling/gridmapping/base.py

@@ -20,6 +20,7 @@
 # DEALINGS IN THE SOFTWARE.
 
 import abc
+from collections.abc import Sequence
 import copy
 import math
 import threading
@@ -61,22 +62,20 @@ class GridMapping(abc.ABC):
     a transformation from image pixel coordinates to spatial Earth coordinates
     defined in a well-known coordinate reference system (CRS).
 
-    This class cannot be instantiated directly. Use one of its factory methods
+    This class cannot be instantiated directly. Use one of its class methods
     to create instances:
 
-    * :meth:`regular`
-    * :meth:`from_dataset`
-    * :meth:`from_coords`
+    * `regular`
+    * `regular_from_bbox`
+    * `from_dataset`
+    * `from_coords`
 
     Some instance methods can be used to derive new instances:
 
-    * :meth:`derive`
-    * :meth:`scale`
-    * :meth:`transform`
-    * :meth:`to_regular`
-
-    This class is thread-safe.
-
+    * `derive`
+    * `scale`
+    * `transform`
+    * `to_regular`
     """
 
     def __init__(
@@ -737,6 +736,36 @@ def regular(
             is_j_axis_up=is_j_axis_up,
         )
 
+    @classmethod
+    def regular_from_bbox(
+        cls,
+        bbox: Sequence[FloatInt],
+        xy_res: float | tuple[float, float],
+        crs: str | pyproj.crs.CRS,
+        tile_size: int | tuple[int, int] = None,
+        is_j_axis_up: bool = False,
+    ) -> "GridMapping":
+        """Creates a regular grid mapping for a given coordinate reference system based
+        on the given bounding box and spatial resolution.
+
+        Args:
+            bbox: Bounding box coordinates in the format [west, south, east, north].
+                The values must be in the given CRS.
+            xy_res: Resolution in x- and y-directions.
+            crs: Coordinate reference system (e.g., "EPSG:4326").
+            tile_size: Chunk size for the grid; if a single int is given,
+                square chunk size is assumed. Defaults to 1024.
+            is_j_axis_up: Whether positive j-axis points up. Defaults to False.
+
+        Returns:
+            A regular grid mapping object.
+        """
+        from .regular import new_regular_grid_mapping_from_bbox
+
+        return new_regular_grid_mapping_from_bbox(
+            bbox, xy_res, crs, tile_size=tile_size, is_j_axis_up=is_j_axis_up
+        )
+
     def to_regular(
         self, tile_size: int | tuple[int, int] | None = None, is_j_axis_up: bool = False
     ) -> "GridMapping":

xcube_resampling/gridmapping/regular.py

@@ -19,7 +19,10 @@
 # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 # DEALINGS IN THE SOFTWARE.
 
+from collections.abc import Sequence
+
 import dask.array as da
+import numpy as np
 import pyproj
 import xarray as xr
 
@@ -33,6 +36,7 @@
     _normalize_number_pair,
     _to_int_or_float,
 )
+from xcube_resampling.constants import FloatInt
 
 
 class RegularGridMapping(GridMapping):
@@ -164,3 +168,38 @@ def to_regular_grid_mapping(
         tile_size=tile_size,
         is_j_axis_up=is_j_axis_up,
     )
+
+
+def new_regular_grid_mapping_from_bbox(
+    bbox: Sequence[FloatInt],
+    xy_res: FloatInt | tuple[FloatInt, FloatInt],
+    crs: str | pyproj.CRS,
+    tile_size: int | tuple[int, int] = 1024,
+    is_j_axis_up: bool = False,
+) -> GridMapping:
+    """Creates a regular grid mapping for a given coordinate reference system based on
+    the given bounding box and spatial resolution.
+
+    Args:
+        bbox: Bounding box coordinates in the format [west, south, east, north].
+            The values must be in the given CRS.
+        spatial_res: Spatial resolution of the grid.
+        crs: Coordinate reference system (e.g., "EPSG:4326").
+        tile_size: Chunk size for the grid; if a single int is given,
+            square chunk size is assumed. Defaults to 1024.
+
+    Returns:
+        A regular grid mapping object.
+    """
+    if xy_res is not tuple:
+        xy_res = (xy_res, xy_res)
+    x_size = int(np.ceil((bbox[2] - bbox[0]) / xy_res[1]))
+    y_size = int(np.ceil(abs(bbox[3] - bbox[1]) / xy_res[0]))
+    return new_regular_grid_mapping(
+        size=(x_size, y_size),
+        xy_min=(bbox[0], bbox[1]),
+        xy_res=xy_res,
+        crs=crs,
+        tile_size=tile_size,
+        is_j_axis_up=is_j_axis_up,
+    )

xcube_resampling/utils.py

@@ -236,7 +236,7 @@ def bbox_overlap(
 ) -> float:
     """
     Calculate the fraction of the source bounding box that overlaps with the target
-    bounding box.
+    bounding box relative to the area of the source bounding box.
 
     Args:
         source_bbox: (min_x, min_y, max_x, max_y)
@@ -245,19 +245,37 @@ def bbox_overlap(
     Returns:
         float in [0, 1]
     """
-    inter_min_x = max(source_bbox[0], target_bbox[0])
-    inter_min_y = max(source_bbox[1], target_bbox[1])
-    inter_max_x = min(source_bbox[2], target_bbox[2])
-    inter_max_y = min(source_bbox[3], target_bbox[3])
-
-    inter_w = max(0, inter_max_x - inter_min_x)
-    inter_h = max(0, inter_max_y - inter_min_y)
-    inter_area = inter_w * inter_h
-    area_source = (source_bbox[2] - source_bbox[0]) * (source_bbox[3] - source_bbox[1])
+    source_bboxes = _split_bbox_antimeridian(source_bbox)
+    target_bboxes = _split_bbox_antimeridian(target_bbox)
+
+    inter_area = 0.0
+    for source_bbox in source_bboxes:
+        for target_bbox in target_bboxes:
+            inter_min_x = max(source_bbox[0], target_bbox[0])
+            inter_min_y = max(source_bbox[1], target_bbox[1])
+            inter_max_x = min(source_bbox[2], target_bbox[2])
+            inter_max_y = min(source_bbox[3], target_bbox[3])
+
+            inter_w = max(0, inter_max_x - inter_min_x)
+            inter_h = max(0, inter_max_y - inter_min_y)
+            inter_area += inter_w * inter_h
+    area_source = 0.0
+    for source_bbox in source_bboxes:
+        area_source += (source_bbox[2] - source_bbox[0]) * (
+            source_bbox[3] - source_bbox[1]
+        )
 
     return inter_area / area_source
 
 
+def _split_bbox_antimeridian(bbox: Sequence[FloatInt]):
+    min_x, min_y, max_x, max_y = bbox
+
+    if max_x < min_x:
+        return [(min_x, min_y, 180, max_y), (-180, min_y, max_x, max_y)]
+    return [bbox]
+
+
 def normalize_grid_mapping(ds: xr.Dataset, gm: GridMapping) -> xr.Dataset:
     """
     Normalize the grid mapping of a dataset to use a standard "spatial_ref" coordinate.

xcube_resampling/version.py

@@ -19,4 +19,4 @@
 # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 # DEALINGS IN THE SOFTWARE.
 
-__version__ = "0.1.1"
+__version__ = "0.2.0.dev0"