Merge pull request #23 from lukasValentin/dev

Patch EOdal v0.1.1: Clip to feature geometry and mask pixels not overlapping it

.github/workflows/test-release-candidate.yaml

@@ -48,10 +48,8 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, windows-latest]
-        python-version: ["3.8", "3.9", "3.10"]
+        python-version: ["3.10"]
         include:
-          - os: macos-latest
-            python-version: "3.8"
           - os: macos-latest
             python-version: "3.10"
     env:

.github/workflows/test.yaml

@@ -23,7 +23,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest]
-        python-version: ["3.8", "3.10"]
+        python-version: ["3.10"]
         include:
           - os: windows-latest
             python-version: "3.10"

CHANGELOG.rst

@@ -10,6 +10,15 @@ The format is based on `Keep a Changelog`_, and this project adheres to `Semanti
 
 Categories for changes are: Added, Changed, Deprecated, Removed, Fixed, Security.
 
+Version `0.1.1 < https://github.com/EOA-team/eodal/releases/tag/v0.1.1>`__
+--------------------------------------------------------------------------------
+
+Release date: 2022-12-13
+
+- Fixed: Band.clip() now has a optional keyword "full_bounding_box_only" (default: False) allowing to mask pixels outside the feature geometry.
+- Fixed: Sentinel.mask_clouds_and_shadows() default SCL classes were updated so that everything but SCL classes 4 and 5 are masked by default.
+- Added: SceneCollection now also allows clipping to a feature geometry (SceneCollection.clip_scenes())
+
 Version `0.1.0 < https://github.com/EOA-team/eodal/releases/tag/v0.1.0>`__
 --------------------------------------------------------------------------------
 

eodal/core/algorithms.py

@@ -0,0 +1,208 @@
+'''
+Collection of algorithms working with EOdal core objects such as Bands,
+RasterCollections, Scenes and SceneCollections.
+
+Copyright (C) 2022 Lukas Valentin Graf
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+'''
+
+import os
+import geopandas as gpd
+import uuid
+
+from pathlib import Path
+from typing import Any, Dict, List, Optional, Tuple, Union
+from rasterio.merge import merge
+from rasterio.crs import CRS
+
+from eodal.config import get_settings
+from eodal.core.band import Band, GeoInfo
+from eodal.core.raster import RasterCollection, SceneProperties
+
+Settings = get_settings()
+
+def _get_crs_and_attribs(
+    in_file: Path, **kwargs
+) -> Tuple[GeoInfo, List[Dict[str, Any]]]:
+    """
+    Returns the ``GeoInfo`` from a multi-band raster dataset
+
+    :param in_file:
+        raster datasets from which to extract the ``GeoInfo`` and
+        attributes
+    :param kwargs:
+        optional keyword-arguments to pass to
+        `~eodal.core.raster.RasterCollection.from_multi_band_raster`
+    :returns:
+        ``GeoInfo`` and metadata attributes of the raster dataset
+    """
+
+    ds = RasterCollection.from_multi_band_raster(fpath_raster=in_file, **kwargs)
+    geo_info = ds[ds.band_names[0]].geo_info
+    attrs = [ds[x].get_attributes() for x in ds.band_names]
+    return geo_info, attrs
+
+def merge_datasets(
+    datasets: List[Path],
+    out_file: Optional[Path] = None,
+    target_crs: Optional[int | CRS] = None,
+    vector_features: Optional[Path | gpd.GeoDataFrame] = None,
+    scene_properties: Optional[SceneProperties] = None,
+    band_options: Optional[Dict[str, Any]] = None,
+    sensor: Optional[str] = None,
+    **kwargs,
+) -> Union[None, RasterCollection]:
+    """
+    Merges a list of raster datasets using the ``rasterio.merge`` module.
+
+    The function can handle datasets in different coordinate systems by resampling
+    the data into a common spatial reference system either provided in the function
+    call or infered from the first dataset in the list.
+
+    ATTENTION:
+        All datasets must have the same number of bands and data type!
+
+    :param datasets:
+        list of datasets (as path-like objects or opened raster datasets)
+        to merge into a single raster
+    :param out_file:
+        name of the resulting raster dataset (optional). If None (default)
+        returns a new ``SatDataHandler`` instance otherwise writes the data
+        to disk as new raster dataset.
+    :param target_crs:
+        optional target spatial coordinate reference system in which the output
+        product shall be generated. Must be passed as integer EPSG code or CRS
+        instance.
+    :param vector_features:
+        optional vector features to clip the merged dataset to (full bounding box).
+    :param scene_properties:
+        optional scene properties to set to the resulting merged dataset
+    :param band_options:
+        optional sensor-specific band options to pass to the sensor's
+        ``RasterCollection`` constructor
+    :param sensor:
+        if the data is from a sensor explicitly supported by eodal such as
+        Sentinel-2 the raster data is loaded into a sensor-specific collection
+    :param kwargs:
+        kwargs to pass to ``rasterio.warp.reproject``
+    :returns:
+        depending on the kwargs passed either `None` (if output is written to file directly)
+        or a `RasterCollection` instance if the operation takes place in memory
+    """
+    # check the CRS and attributes of the datasets first
+    crs_list = []
+    attrs_list = []
+    for dataset in datasets:
+        geo_info, attrs = _get_crs_and_attribs(in_file=dataset)
+        crs_list.append(geo_info.epsg)
+        attrs_list.append(attrs)
+
+    if target_crs is None:
+        # use CRS from first dataset
+        target_crs = crs_list[0]
+    # coordinate systems are not the same -> re-projection of raster datasets
+    if len(set(crs_list)) > 1:
+        pass
+    # all datasets have one coordinate system, check if it is the desired one
+    else:
+        if target_crs is not None:
+            if crs_list[0] != target_crs:
+                # re-projection into target CRS required
+                pass
+
+    # use rasterio merge to get a new raster dataset
+    dst_kwds = {"QUALITY": "100", "REVERSIBLE": "YES"}
+    try:
+        res = merge(datasets=datasets, dst_path=out_file, dst_kwds=dst_kwds, **kwargs)
+        if res is not None:
+            out_ds, out_transform = res[0], res[1]
+    except Exception as e:
+        raise Exception(f"Could not merge datasets: {e}")
+
+    # when out_file was provided the merged data is written to file directly
+    if out_file is not None:
+        return
+    # otherwise, create new SatDataHandler instance from merged datasets
+    # add scene properties if available
+    if sensor is None:
+        raster = RasterCollection(scene_properties=scene_properties)
+    else:
+        raster = eval(f'eodal.core.sensors.{sensor.lower()}.{sensor[0].upper() + sensor[1::]}' \
+                      + '(scene_properties=scene_properties)')
+    n_bands = out_ds.shape[0]
+    # take attributes of the first dataset
+    attrs = attrs_list[0]
+    geo_info = GeoInfo(
+        epsg=target_crs,
+        ulx=out_transform.c,
+        uly=out_transform.f,
+        pixres_x=out_transform.a,
+        pixres_y=out_transform.e,
+    )
+    for idx in range(n_bands):
+        band_attrs = attrs[idx]
+        nodata = band_attrs.get("nodatavals")
+        if isinstance(nodata, tuple):
+            nodata = nodata[0]
+        is_tiled = band_attrs.get("is_tiled")
+        scale = band_attrs.get("scales")
+        if isinstance(scale, tuple):
+            scale = scale[0]
+        offset = band_attrs.get("offsets")
+        if isinstance(offset, tuple):
+            offset = offset[0]
+        description = band_attrs.get("descriptions")
+        if isinstance(description, tuple):
+            description = description[0]
+        unit = band_attrs.get("units")
+        if isinstance(unit, tuple):
+            unit = unit[0]
+        raster.add_band(
+            band_constructor=Band,
+            band_name=f"B{idx+1}",
+            values=out_ds[idx, :, :],
+            geo_info=geo_info,
+            is_tiled=is_tiled,
+            scale=scale,
+            offset=offset,
+            band_alias=description,
+            unit=unit,
+        )
+
+    # clip raster collection if required to vector_features
+    if vector_features is not None:
+        tmp_dir = Settings.TEMP_WORKING_DIR
+        fname_tmp = tmp_dir.joinpath(f"{uuid.uuid4()}.tif")
+        raster.to_rasterio(fname_tmp)
+        if sensor is None:
+            expr = "RasterCollection"
+        else:
+            expr = f"eodal.core.sensors.{sensor.lower()}.{sensor[0].upper() + sensor[1::]}()"
+        if band_options is None:
+            band_options = {}
+        raster = eval(
+            f"""{expr}.from_multi_band_raster(
+            fpath_raster=fname_tmp,
+            vector_features=vector_features,
+            full_bounding_box_only=False,
+            **band_options
+        )"""
+        )
+        # set scene properties if available
+        if scene_properties is not None:
+            raster.scene_properties = scene_properties
+        os.remove(fname_tmp)
+
+    return raster

eodal/core/band.py

@@ -1118,19 +1118,31 @@ def copy(self):
 
     def clip(
         self,
-        clipping_bounds: Path | gpd.GeoDataFrame | Tuple[float,float,float,float] | Polygon,
+        clipping_bounds: Path | gpd.GeoDataFrame | Tuple[float,float,float,float] | Polygon  | MultiPolygon,
+        full_bounding_box_only: Optional[bool] = False,
         inplace: Optional[bool] = False
     ):
         """
-        Clip a band object to a spatial extent.
+        Clip a band object to a geometry or the bounding box of one or more
+        geometries. By default, pixel values outside the geometry are masked.
+        The spatial extent of the returned `Band` instance is **always** cropped
+        to the bounding box of the geomtry/ geometries.
+
+        NOTE:
+            When passing a `GeoDataFrame` with more than one feature, the single
+            feature geometries are dissolved into a single one!
 
         :param clipping_bounds:
-            spatial bounds to clip the Band to. Only clipping to rectangular shapes
-            is supported. Can be either a vector file, a shapely `Polygon`, a
-            `GeoDataFrame` or a coordinate tuple with (xmin, ymin, xmax, ymax).
+            spatial bounds to clip the Band to. Can be either a vector file, a shapely
+            `Polygon` or `MultiPolygon`, a `GeoDataFrame` or a coordinate tuple with
+            (xmin, ymin, xmax, ymax).
             Vector files and `GeoDataFrame` are reprojected into the bands' coordinate
             system if required, while the coordinate tuple and shapely geometry **MUST**
             be provided in the CRS of the band.
+        :param full_bounding_box_only:
+            if False (default), clips to the bounding box of the geometry and masks values
+            outside the actual geometry boundaries. To obtain all values within the
+            bounding box set to True.
         :param inplace:
             if False (default) returns a copy of the ``Band`` instance
             with the changes applied. If True overwrites the values
@@ -1140,18 +1152,24 @@ def clip(
         """
         if isinstance(clipping_bounds, Path):
             clipping_bounds = gpd.read_file(clipping_bounds)
+
         # check inputs
         if isinstance(clipping_bounds, tuple):
             if len(clipping_bounds) != 4:
                 raise ValueError('Expected four coordinates (xmin, ymin, xmax, ymax)')
             xmin, ymin, xmax, ymax = clipping_bounds
+            actual_geom = box(*clipping_bounds)
         elif isinstance(clipping_bounds, gpd.GeoDataFrame):
             # get the bounding box of the FIRST feature
             _clipping_bounds = clipping_bounds.copy()
             _clipping_bounds = _clipping_bounds.bounds
-            xmin, ymin, xmax, ymax = list(clipping_bounds)
-        elif isinstance(clipping_bounds, Polygon):
+            xmin, ymin, xmax, ymax = _clipping_bounds.values[0]
+            # the actual geometries are dissolved in case there is more than one record
+            # and converted to a shapely object
+            actual_geom = clipping_bounds.dissolve().geometry.values[0]
+        elif (isinstance(clipping_bounds, Polygon) or isinstance(clipping_bounds, MultiPolygon)):
             xmin, ymin, xmax, ymax = clipping_bounds.bounds
+            actual_geom = clipping_bounds
         else:
             raise TypeError(f'{type(clipping_bounds)} is not supported')
         # actual clipping operation. Calculate the rows and columns where to clip
@@ -1205,6 +1223,26 @@ def clip(
         values = self.values.copy()
         new_values = values[min_row:max_row,min_col:max_col]
 
+        # if full_bounding_box is False, mask out pixels not overlapping the
+        # geometry (but located within the bounding box)
+        if not full_bounding_box_only:
+            # to mask pixels outside the geometry we need to rasterize it
+            # Rasterize vector using the shape and coordinate system of the raster
+            mask = features.rasterize(
+                [actual_geom],
+                out_shape=new_values.shape,
+                fill=1,
+                out=None,
+                transform=new_geo_info.as_affine(),
+                all_touched=True,
+                default_value=0,
+                dtype='uint8'
+            ).astype(bool)
+            new_values = np.ma.MaskedArray(
+                data=new_values,
+                mask=mask
+            )
+
         if inplace:
             object.__setattr__(self, "values", new_values)
             object.__setattr__(self, "geo_info", new_geo_info)

eodal/core/raster.py

@@ -957,6 +957,16 @@ def clip_bands(
     ):
         """
         Clip bands in RasterCollection to a user-defined spatial bounds.
+
+        :param band_selection:
+            optional list of bands to clip. If not provided takes all available
+            bands.
+        :param inplace:
+            if False (default) returns a copy of the ``RasterCollection`` instance
+            with the changes applied. If True overwrites the values
+            in the current instance.
+        :param **kwargs:
+            key-word arguments to pass to `Band.clip` method.
         """
         if band_selection is None:
             band_selection = self.band_names

eodal/core/scene.py

@@ -404,6 +404,39 @@ def apply(self, func: Callable, *args, **kwargs) -> Any:
         except Exception as e:
             raise ValueError from e
 
+    def clip_scenes(
+        self,
+        inplace: Optional[bool] = False,
+        **kwargs
+    ):
+        """
+        Clip scenes in a SceneCollection to a user-defined spatial bounds.
+
+        :param inplace:
+            if False (default) returns a copy of the ``SceneCollection`` instance
+            with the changes applied. If True overwrites the values
+            in the current instance.
+        :param **kwargs:
+            key-word arguments to pass to `RasterCollection.clip_bands` method.
+        """
+        # loop over Scenes and try to subset them spatially
+        # initialize a new SceneCollection if inplace is False
+        scoll_new = None
+        if inplace:
+            kwargs.update({'inplace': True})
+        if not inplace:
+            scoll_new = SceneCollection(indexed_by_timestamps=self.indexed_by_timestamps)
+        # loop over band reproject the selected ones
+        for timestamp, scene in self:
+            if inplace:
+                self[timestamp].clip_bands(**kwargs)
+            else:
+                scene = self[timestamp].copy()
+                scoll_new.add_scene(scene_constructor=scene.clip_bands, **kwargs)
+
+        if not inplace:
+            return scoll_new
+
     def copy(self):
         """returns a true copy of the SceneCollection"""
         return deepcopy(self)