ENH: add support to generate sar rgb images, S1 despeckling with enha…

…nced lee filter and dB images (theroggy#157)

CHANGELOG.md

@@ -34,11 +34,13 @@
 - Filter away rasterio logging for extrasamples (#127)
 - Save openeo images with int16 bands as int16 locally as well (again) (#131, #135)
 - Linting improvements: add isort, sync rules with geofileops (#133, #134)
-- Add support for zonalstats calculation with ExactExctract (#139)
+- Add support to calculate zonalstats/spatial aggregations with
+  [ExactExtract](https://github.com/isciences/exactextract) (#139)
 - Make the number of top classes of the classification to retain configurable (#145)
 - Add support to configure pixel type for indices (#150)
-- Add support for rvi index (#146)
-- Add support for vvdvh index (#151)
+- Add support to generate rvi index (#146)
+- Add support to generate vvdvh index (#151)
+- Add support to generate sarrgb images and apply enhanced lee despeckling (#157)
 - General small improvements, e.g. save randomforest models compressed,.. (#144)
 
 ### Bugs fixed

bin_util/calc_periodic_mosaic.py

@@ -12,8 +12,8 @@ def main():
     # Init some variables
     roi_crs = 31370
     # BEFL
-    start_date = datetime(2023, 7, 1)
-    end_date = datetime(2023, 12, 31)
+    start_date = datetime(2024, 7, 1)
+    end_date = datetime(2024, 9, 30)
     roi_bounds = [20_000, 150_000, 260_000, 245_000]
     images_periodic_dir = Path("//dg3.be/alp/Datagis/satellite_periodic/BEFL")
 
@@ -25,9 +25,19 @@ def main():
     images_periodic_dir = Path("c:/temp/periodic_mosaic/roi_test")
     """
 
-    imageprofiles_to_get = ["s1-rvi-asc-weekly", "s1-rvi-desc-weekly"]
-    # imageprofiles_to_get = ["s1-dprvi-asc-weekly", "s1-dprvi-desc-weekly"]
-    # imageprofiles_to_get = ["s2-agri-weekly"]
+    imageprofiles_to_get = [
+        # "s2-agri-weekly",
+        # "s2-ndvi-weekly",
+        # "s1-grd-sigma0-asc-weekly",
+        # "s1-grd-sigma0-desc-weekly",
+        # "s1-coh-weekly",
+        # "s1-grd-sigma0-vvdvh-asc-weekly",
+        # "s1-grd-sigma0-vvdvh-desc-weekly",
+        # "s1-sarrgbdb-asc-weekly",
+        "s1-sarrgbdb-desc-weekly",
+        # "s1-rvi-asc-weekly", "s1-rvi-desc-weekly",
+        # "s1-dprvi-asc-weekly", "s1-dprvi-desc-weekly",
+    ]
 
     image_profiles_path = (
         Path(__file__).resolve().parent.parent

ci/envs/latest.yml

@@ -5,20 +5,20 @@ dependencies:
   - python
   - pip
   # required
-  - fiona
+  - dask-core
   - geofileops
-  - geopandas>=0.11
+  - geopandas-base>=0.11
   - h5py
-  - numpy
+  - numexpr
+  - numpy<2
   - openeo
   - pandas
   - psutil
   - pyproj
-  - scikit-learn
   - rasterio
   - rasterstats
   - rioxarray
-  - numpy<2
+  - scikit-learn
   # Optional
   - qgis
   # testing

ci/envs/minimal.yml

@@ -3,22 +3,22 @@ channels:
   - conda-forge
 dependencies:
   - python
+  - pip
   # required
-  - fiona
+  - dask-core
   - geofileops
-  - geopandas>=0.11
+  - geopandas-base>=0.11
   - h5py
-  - numpy
+  - numexpr
+  - numpy<2
   - openeo
   - pandas
   - psutil
   - pyproj
-  - scikit-learn
   - rasterio
   - rasterstats
   - rioxarray
-  - numpy<2
-  - pip
+  - scikit-learn
   # Optional
   # - qgis
   # testing

cropclassification/util/lee_enhanced.py

@@ -0,0 +1,111 @@
+"""Implementation of the Enhanced Lee Filter for speckle reduction in SAR images.
+
+References:
+  - https://stackoverflow.com/questions/4959171/improving-memory-usage-in-an-array-wide-filter-to-avoid-block-processing
+  - https://catalyst.earth/catalyst-system-files/help/concepts/orthoengine_c/Chapter_825.html
+"""
+
+import warnings
+
+import numexpr as ne
+import numpy as np
+import rasterio as rio
+import scipy.ndimage
+import scipy.signal
+
+from . import raster_util
+
+
+def _moving_average(image, size):
+    Im = np.empty(image.shape, dtype=np.float32)
+    # scipy.ndimage.filters.uniform_filter(image, filtsize, output=Im)
+    # scipy.ndimage.generic_filter(image, function=np.nanmean, size=size, output=Im)
+    Im = _filter_nanmean(image, size=size)
+    return Im
+
+
+def _moving_stddev(image, size):
+    Im = np.empty(image.shape, dtype=np.float32)
+    # scipy.ndimage.filters.uniform_filter(image, filtersize, output=Im)
+    # scipy.ndimage.generic_filter(image, function=np.nanmean, size=size, output=Im)
+    Im = _filter_nanmean(image, size=size)
+    Im = ne.evaluate("((image-Im) ** 2)")
+    # scipy.ndimage.filters.uniform_filter(Im, filtersize, output=Im)
+    # scipy.ndimage.generic_filter(Im, function=np.nanmean, size=size, output=Im)
+    Im = _filter_nanmean(Im, size=size)
+    return ne.evaluate("sqrt(Im)")
+
+
+def _filter_nanmean(image, size):
+    kernel = np.ones((size, size))
+    kernel[1, 1] = 0
+
+    neighbor_sum = scipy.signal.convolve2d(
+        image, kernel, mode="same", boundary="fill", fillvalue=0
+    )
+
+    num_neighbor = scipy.signal.convolve2d(
+        np.ones(image.shape), kernel, mode="same", boundary="fill", fillvalue=0
+    )
+
+    return neighbor_sum / num_neighbor
+
+
+def lee_enhanced(
+    image, filtersize: int = 5, nlooks: float = 10.0, dfactor: float = 10.0
+):
+    # Implementation based on PCI Geomatimagea's FELEE function documentation
+    with warnings.catch_warnings():
+        warnings.filterwarnings(
+            "ignore", category=RuntimeWarning, message="Mean of empty slice"
+        )
+        Ci = _moving_stddev(image, filtersize)
+        Im = _moving_average(image, filtersize)
+
+    Ci /= Im
+
+    Cu = np.sqrt(1 / nlooks).astype(np.float32)  # noqa: F841
+    Cmax = np.sqrt(1 + (2 * nlooks)).astype(np.float32)  # noqa: F841
+
+    W = ne.evaluate("exp(-dfactor * (Ci - Cu) / (Cmax - Ci))")
+    If = ne.evaluate("Im * W + image * (1 - W)")
+    del W
+
+    out = ne.evaluate("where(Ci <= Cu, Im, If)")
+    del Im
+    del If
+
+    out = ne.evaluate("where(Ci >= Cmax, image, out)")
+    return out
+
+
+def lee_enhanced_file(
+    input_path,
+    output_path,
+    filtersize: int = 5,
+    nlooks: float = 10.0,
+    dfactor: float = 10.0,
+    force: bool = False,
+):
+    if output_path.exists():
+        if force:
+            output_path.unlink()
+        else:
+            return
+
+    with rio.open(input_path) as input:
+        profile = input.profile
+        band1 = input.read(1)
+        band1_lee = lee_enhanced(
+            band1, filtersize=filtersize, nlooks=nlooks, dfactor=dfactor
+        )
+        band2 = input.read(2)
+        band2_lee = lee_enhanced(
+            band2, filtersize=filtersize, nlooks=nlooks, dfactor=dfactor
+        )
+    band_descriptions = raster_util.get_band_descriptions(input_path)
+
+    with rio.open(output_path, "w", **profile) as dst:
+        dst.write(band1_lee, 1)
+        dst.write(band2_lee, 2)
+    raster_util.set_band_descriptions(output_path, band_descriptions.keys())

cropclassification/util/mosaic_util.py

@@ -263,6 +263,7 @@ def calc_periodic_mosaic(
                 image_local["path"],
                 index=image_local["index_type"],
                 pixel_type=image_local["pixel_type"],
+                process_options=image_local["process_options"],
                 force=force,
             )
 

cropclassification/util/openeo_util.py

@@ -555,6 +555,7 @@ def get_date_score(day) -> openeo.rest.datacube.DataCube:
 
     kernel_size = int(150 * 20 / spatial_resolution) + 1
     gaussian_1d = gaussian(M=kernel_size, std=1)
+    assert isinstance(gaussian_1d, np.ndarray)
     gaussian_kernel = np.outer(gaussian_1d, gaussian_1d)
     gaussian_kernel /= gaussian_kernel.sum()