Merge pull request #5 from lukasValentin/dev

Several bugfixes and improvements

eodal/core/band.py

@@ -1775,6 +1775,7 @@ def reduce(
         self,
         method: Union[str, List[str]],
         by: Optional[Union[Path, gpd.GeoDataFrame]] = None,
+        method_args: Optional[Dict[str, Any]] = None
     ) -> Dict[str, Union[int, float]]:
         """
         Reduces the raster data to scalar values.
@@ -1783,6 +1784,12 @@ def reduce(
             any ``numpy`` function taking a two-dimensional array as input
             and returning a single scalar. Can be a single function name
             (e.g., "mean") or a list of function names (e.g., ["mean", "median"])
+        :param by:
+            define by what to reduce the band values (not implemented yet!!)
+        :param method_args:
+            optional dictionary with arguments to pass to the single methods in
+            case the reducer method requires extra arguments to function properly
+            (e.g., `np.quantile`)
         :returns:
             a dictionary with scalar results
         """
@@ -1813,7 +1820,11 @@ def reduce(
             try:
                 # get function object and use its __call__ method
                 numpy_function = eval(expression)
-                stats[operator] = numpy_function.__call__(self.values)
+                # check if there are any function arguments
+                args = []
+                if method_args is not None:
+                    args = method_args.get(method, None)
+                stats[operator] = numpy_function.__call__(self.values, *args)
             except TypeError:
                 raise Exception(f"Unknown function name for {numpy_prefix}: {operator}")
 

eodal/core/sensors/sentinel1.py

@@ -36,10 +36,11 @@
 from eodal.core.band import Band
 from eodal.core.raster import RasterCollection
 from eodal.core.scene import SceneProperties
+from eodal.utils.decorators import prepare_point_features
 from eodal.utils.sentinel1 import get_S1_platform_from_safe, \
     get_S1_acquistion_time_from_safe, _url_to_safe_name, \
     get_s1_imaging_mode_from_safe
-from utils.exceptions import DataNotFoundError
+from eodal.utils.exceptions import DataNotFoundError
 
 Settings = get_settings()
 
@@ -151,6 +152,7 @@ def from_safe(
         return sentinel1
 
     @classmethod
+    @prepare_point_features
     def read_pixels_from_safe(
             cls,
             vector_features: Path | gpd.GeoDataFrame,

eodal/core/sensors/sentinel2.py

@@ -49,6 +49,7 @@
     s2_gain_factor,
     SCL_Classes,
 )
+from eodal.utils.decorators import prepare_point_features
 from eodal.utils.exceptions import BandNotFoundError
 from eodal.utils.sentinel2 import (
     get_S2_bandfiles_with_res,
@@ -471,6 +472,7 @@ def from_safe(
         return sentinel2
 
     @classmethod
+    @prepare_point_features
     def read_pixels_from_safe(
         cls,
         vector_features: Union[Path, gpd.GeoDataFrame],
@@ -565,23 +567,27 @@ def read_pixels_from_safe(
                 if apply_scaling:
                     gdf_scaled = gdf_band.copy()
                     gdf_scaled[band_name] = 0.0
-                    gdf_scaled[band_name] = (
-                        offset + gdf_band[band_name].loc[gdf_band[band_name] != 0]
-                    ) * gain
+                    # use only pixel values were reflectance is != 0
+                    gdf_scaled[band_name] = gdf_band[band_name].apply(
+                        lambda x, offset=offset, gain=gain:
+                        (offset + x) * gain if x != 0 else 0 
+                    )
                     band_gdfs.append(gdf_scaled)
                     continue
             band_gdfs.append(gdf_band)
 
-        # concat the single GeoDataFrames with the band data
+        # concatenate the single GeoDataFrames with the band data
         gdf = pd.concat(band_gdfs, axis=1)
         # clean the dataframe and remove duplicate column names after merging
         # to avoid (large) redundancies
         gdf = gdf.loc[:, ~gdf.columns.duplicated()]
         # skip all pixels with zero reflectance (either blackfilled or outside of the
         # scene extent); in case of dtype float check for NaNs
-        if (gdf.dtypes[gdf.columns.str.startswith("B")] == "float64").all():
+        band_names = gdf.columns[gdf.columns.str.startswith("B")]
+        if (gdf.dtypes[band_names] == "float64").all():
+            gdf[band_names] = gdf[band_names].replace({0., np.nan})
             gdf.dropna(axis=0, inplace=True)
-        elif (gdf.dtypes[gdf.columns.str.startswith("B")] == "int16").all():
+        elif (gdf.dtypes[band_names] == "int16").all():
             gdf = gdf.loc[~(gdf[band_df_safe.band_name] == 0).all(axis=1)]
 
         return gdf

eodal/core/utils/geometry.py

@@ -48,7 +48,6 @@ def read_geometries(in_dataset: Union[Path, gpd.GeoDataFrame]) -> gpd.GeoDataFra
             f"Could not read geometries of input type {type(in_dataset)}"
         )
 
-
 def check_geometry_types(
     in_dataset: Union[Path, gpd.GeoDataFrame],
     allowed_geometry_types: List[str],
@@ -96,7 +95,6 @@ def check_geometry_types(
         )
     return gdf
 
-
 def convert_3D_2D(geometry: gpd.GeoSeries) -> gpd.GeoSeries:
     """
     Takes a GeoSeries of 3D Multi/Polygons (has_z) and returns a list of 2D Multi/Polygons.
@@ -127,3 +125,22 @@ def convert_3D_2D(geometry: gpd.GeoSeries) -> gpd.GeoSeries:
             new_geo = geometry
             break
     return new_geo
+
+def multi_to_single_points(point_features: gpd.GeoDataFrame | Path) -> gpd.GeoDataFrame:
+    """
+    Casts MultiPoint geometries to single point geometries by calling
+    `gpd.GeoDataFrame.explode()`
+
+    :param point_features:
+        point features to cast
+    :returns:
+        casted point features or input if all geometries are already single parted
+    """
+    gdf = check_geometry_types(
+        in_dataset=point_features,
+        allowed_geometry_types=['Point', 'MultiPoint']
+    )
+    if (gdf.geometry.type == 'MultiPoint').any():
+            gdf = gdf.explode(index_parts=False)
+    return gdf
+

eodal/metadata/database/querying.py

@@ -20,12 +20,11 @@
 import geopandas as gpd
 
 from sqlalchemy import create_engine
-
-from eodal.metadata.database import Regions
-from eodal.utils.exceptions import RegionNotFoundError
-from eodal.config import get_settings
 from sqlalchemy.orm.session import sessionmaker
 
+from eodal.config import get_settings
+from eodal.metadata.database import Regions, S2_Raw_Metadata
+from eodal.utils.exceptions import DataNotFoundError, RegionNotFoundError
 
 Settings = get_settings()
 logger = Settings.logger
@@ -44,17 +43,37 @@ def get_region(region: str) -> gpd.GeoDataFrame:
         unique region identifier
 
     :return:
-        geodataframe with the geometry of the queried region
+        `GeoDataFrame` with the geometry of the queried region
     """
-
     query_statement = (
         session.query(Regions.geom, Regions.region_uid)
         .filter(Regions.region_uid == region)
         .statement
     )
-
     try:
         return gpd.read_postgis(query_statement, session.bind)
-
     except Exception as e:
         raise RegionNotFoundError(f"{region} not found: {e}")
+
+def get_s2_tile_footprint(tile_name: str) -> gpd.GeoDataFrame:
+    """
+    Queries the geographic extent of a Sentinel-2 tile
+
+    :param sensor:
+        name of the sensor the tiling scheme belongs to (e.g.,
+        'sentinel2')
+    :param tile_name:
+        name of the tile in the tiling scheme (e.g., 'T32TMT')
+    :returns:
+        extent of the tile in geographic coordinates (WGS84)
+    """
+    query_statement = (
+        session.query(S2_Raw_Metadata.geom)
+        .filter(S2_Raw_Metadata.tile_id == tile_name)
+        .distinct()
+        .statement
+    )
+    try:
+        return gpd.read_postgis(query_statement, session.bind)
+    except Exception as e:
+        raise DataNotFoundError(f"{tile_name} not found: {e}")

eodal/operational/mapping/mapper.py

@@ -539,7 +539,8 @@ def _get_observation(
                 in_dir = scenes_date["real_path"].iloc[0]
             # if there is only one scene all we have to do is to read
             # read pixels in case the feature's dtype is point
-            if feature_dict["features"][0]["geometry"]["type"] == "Point":
+            feature_geom = feature_dict["features"][0]["geometry"]["type"]
+            if feature_geom in ["Point", "MultiPoint"]:
                 if sensor.lower() == 'sentinel1':
                     res = Sentinel1.read_pixels_from_safe(
                         in_dir=in_dir,
@@ -552,9 +553,14 @@ def _get_observation(
                         band_selection=self.mapper_configs.band_names,
                         **kwargs,
                     )
-                res["sensing_date"] = scenes_date["sensing_date"].values
-                res["scene_id"] = scenes_date["scene_id"].values
-                res['sensing_time'] = scenes_date['sensing_time'].values
+                res["sensing_date"] = scenes_date["sensing_date"].values[0]
+                res["scene_id"] = scenes_date["scene_id"].values[0]
+                res['sensing_time'] = scenes_date['sensing_time'].values[0]
+                # make sure the result is projected into the target EPSG code, otherwise
+                # the resulting GeoDataFrame contains wrong coordinates, i.e., the
+                # coordinates were not projected into the target CRS
+                if res.crs.to_epsg() != scenes_date.target_crs.unique()[0]:
+                    res.to_crs(epsg=scenes_date.target_crs.unique()[0], inplace=True)
             # or the feature
             else:
                 if sensor.lower() == 'sentinel1':

eodal/operational/mapping/sentinel2.py

@@ -189,7 +189,8 @@ def _read_multiple_scenes(
         # if the feature is a point we take the data set that is not blackfilled.
         # If more than one data set is not blackfilled  we simply take the
         # first data set
-        if feature_gdf["geometry"].iloc[0].type == "Point":
+        feature_geom = feature_gdf["geometry"].iloc[0].type
+        if feature_geom in ["Point", "MultiPoint"]:
             for _, candidate_scene in scenes_date.iterrows():
                 if settings.USE_STAC:
                     in_dir = candidate_scene["assets"]
@@ -204,6 +205,10 @@ def _read_multiple_scenes(
                 if feature_gdf.empty:
                     continue
                 res = feature_gdf
+                # make sure the coordinates are reprojected if necessary
+                if res.crs.to_epsg() != scenes_date.target_crs.unique()[0]:
+                    res.to_crs(epsg=scenes_date.target_crs.unique()[0], inplace=True)
+
                 if isinstance(candidate_scene, (pd.Series, gpd.GeoSeries)):
                     res["sensing_date"] = candidate_scene["sensing_date"]
                     res['sensing_time'] = candidate_scene['sensing_time']

eodal/utils/decorators.py

@@ -24,8 +24,10 @@
 from rasterio.coords import BoundingBox
 
 from eodal.config import get_settings
+from eodal.core.utils.geometry import multi_to_single_points
 from eodal.utils.exceptions import UnknownProcessingLevel, BandNotFoundError
 from eodal.utils.geometry import box_to_geojson
+from eodal.core.utils.geometry import multi_to_single_points
 
 Settings = get_settings()
 
@@ -34,7 +36,7 @@ def prepare_bbox(f):
     @wraps(f)
     def wrapper(**kwargs):
         # a bounding box (vector features) is required
-        vector_features = kwargs.get('vector_features', None)
+        vector_features = kwargs.get('bounding_box', None)
         if vector_features is None:
             raise ValueError('A bounding box must be specified')
         if isinstance(vector_features, Path):
@@ -44,10 +46,33 @@ def wrapper(**kwargs):
         # and provide bounds as geojson (required by STAC)
         bbox = box_to_geojson(gdf=vector_features)
         kwargs.update({'bounding_box': bbox})
-        del kwargs['vector_features']
         return f(**kwargs)
     return wrapper
 
+def prepare_point_features(f):
+    """
+    casts MultiPoint geometries to single parts before calling pixel extraction methods
+    """
+    @wraps(f)
+    def wrapper(*args, **kwargs):
+        if len(args) >= 2:
+            vector_features = args[1]
+        else:
+            vector_features = kwargs.get('vector_features')
+        # cast to single point geometries
+        try:
+            vector_features_updated = multi_to_single_points(vector_features)
+        except Exception as e:
+            print(e)
+        if len(args) >= 2:
+            arg_list = list(args)
+            arg_list[1] = vector_features_updated
+            args = tuple(arg_list)
+        else:
+            kwargs.update({'vector_features': vector_features_updated})
+        return f(*args, **kwargs)
+    return wrapper
+
 def check_processing_level(f):
     @wraps(f)
     def wrapper(*args, **kwargs):