Add crop operator (#112)

## Purpose

Provide a cropping operator that maintains the metadata up to date. This is required when using an operator that relies on the metadata such the regridding to an automatically computed output grid. 

## Code changes:

- Added `idpi.operators.crop` module.

src/idpi/operators/crop.py

@@ -0,0 +1,80 @@
+"""Horizontal cropping operator."""
+
+# Standard library
+import typing
+
+# Third-party
+import numpy as np
+import xarray as xr
+
+# Local
+from .. import metadata
+from . import gis
+
+
+class Bounds(typing.NamedTuple):
+    xmin: int
+    xmax: int
+    ymin: int
+    ymax: int
+
+
+def crop(field: xr.DataArray, bounds: Bounds) -> xr.DataArray:
+    """Crop the field to the given bounds.
+
+    Only fields defined on regular grids in rotlatlon coordinates,
+    without rotation nor flipped axes are supported.
+
+    Parameters
+    ----------
+    field : xarray.DataArray
+        The field to crop.
+    bounds : Bounds
+        Bounds of the cropped area given as indices of the array in following order:
+        xmin, xmax, ymin, ymax.
+        All bounds are inclusive.
+
+    Raises
+    ------
+    ValueError
+        If there are any consistency issues with the provided bounds
+        or any of the conditions on the input grid not met.
+
+    Returns
+    -------
+    xarray.DataArray
+        The data is set to cropped domain and the metadata is updated accordingly.
+
+    """
+    xmin, xmax, ymin, ymax = bounds
+
+    sizes = field.sizes
+    if (
+        xmin > xmax
+        or ymin > ymax
+        or any(v < 0 for v in bounds)
+        or xmax >= sizes["x"]
+        or ymax >= sizes["y"]
+    ):
+        raise ValueError(f"Inconsistent bounds: {bounds}")
+
+    grid = gis.get_grid(field.geography)
+    lon_min, lon_max = np.round(grid.rlon.isel(x=[xmin, xmax]).values * 1e6)
+    lat_min, lat_max = np.round(grid.rlat.isel(y=[ymin, ymax]).values * 1e6)
+    ni = xmax - xmin + 1
+    nj = ymax - ymin + 1
+    npts = ni * nj
+
+    return xr.DataArray(
+        field.isel(x=slice(xmin, xmax + 1), y=slice(ymin, ymax + 1)),
+        attrs=metadata.override(
+            field.message,
+            longitudeOfFirstGridPoint=lon_min,
+            longitudeOfLastGridPoint=lon_max,
+            Ni=ni,
+            latitudeOfFirstGridPoint=lat_min,
+            latitudeOfLastGridPoint=lat_max,
+            Nj=nj,
+            numberOfDataPoints=npts,
+        ),
+    )

src/idpi/operators/support_operators.py

@@ -105,6 +105,11 @@ def replace_vertical(items):
 def get_grid_coords(n: int, x0: float, dx: float, dim: str) -> xr.DataArray:
     """Compute coordinates for an equally spaced grid.
 
+    Values are rounded to the 6th decimal because the data representation
+    in the GRIB specification calls for integers values in microdegrees.
+    It is assumed that data input will in the GRIB format and thus subject
+    to this property.
+
     Parameters
     ----------
     n : int
@@ -122,4 +127,5 @@ def get_grid_coords(n: int, x0: float, dx: float, dim: str) -> xr.DataArray:
         A 1-D field containing the coordinates of the grid along the given dimension
 
     """
-    return xr.DataArray(np.arange(n, dtype=np.float32) * dx + x0, dims=dim)
+    values = np.arange(n) * dx + x0
+    return xr.DataArray(np.round(values, 6), dims=dim)

tests/test_idpi/test_crop.py

@@ -0,0 +1,33 @@
+# Third-party
+from numpy.testing import assert_equal
+
+# First-party
+from idpi.grib_decoder import GribReader
+from idpi.operators import crop, gis
+
+
+def test_crop(data_dir):
+    cdatafile = data_dir / "COSMO-1E/1h/const/000/lfff00000000c"
+
+    reader = GribReader.from_files([cdatafile])
+    ds = reader.load_fieldnames(["HHL"])
+    hhl = ds["HHL"]
+
+    observed = crop.crop(hhl, crop.Bounds(1, 2, 3, 5))
+
+    grid = gis.get_grid(hhl.geography)
+    cropped = hhl.assign_coords(x=grid.rlon, y=grid.rlat).isel(x=[1, 2], y=[3, 4, 5])
+
+    expected_values = cropped.values
+    expected_geography = hhl.geography | {
+        "Ni": 2,
+        "Nj": 3,
+        "longitudeOfFirstGridPointInDegrees": cropped.coords["x"].min().item(),
+        "longitudeOfLastGridPointInDegrees": cropped.coords["x"].max().item(),
+        "latitudeOfFirstGridPointInDegrees": cropped.coords["y"].min().item(),
+        "latitudeOfLastGridPointInDegrees": cropped.coords["y"].max().item(),
+    }
+
+    assert_equal(observed.values, expected_values)
+    assert observed.geography == expected_geography
+    assert observed.parameter == hhl.parameter