257 add smoothing/filter functions (#289)

docs/raster_processing/derivatives/classification.md

@@ -0,0 +1,3 @@
+# Classification
+
+::: eis_toolkit.raster_processing.derivatives.classification

docs/raster_processing/derivatives/parameters.md

@@ -0,0 +1,3 @@
+# Parameters
+
+::: eis_toolkit.raster_processing.derivatives.parameters

docs/raster_processing/derivatives/partial_derivatives.md

@@ -0,0 +1,3 @@
+# Partial derivatives
+
+::: eis_toolkit.raster_processing.derivatives.partial_derivatives

docs/raster_processing/derivatives/utilities.md

@@ -0,0 +1,3 @@
+# Utilities
+
+::: eis_toolkit.raster_processing.derivatives.utilities

docs/raster_processing/filters/focal.md

@@ -0,0 +1,3 @@
+# Focal
+
+::: eis_toolkit.raster_processing.filters.focal

docs/raster_processing/filters/kernels.md

@@ -0,0 +1,3 @@
+# Kernels
+
+::: eis_toolkit.raster_processing.filters.kernels

docs/raster_processing/filters/speckle.md

@@ -0,0 +1,3 @@
+# Speckle
+
+::: eis_toolkit.raster_processing.filters.speckle

docs/raster_processing/filters/utilities.md

@@ -0,0 +1,3 @@
+# Utilities
+
+::: eis_toolkit.raster_processing.filters.utilities

eis_toolkit/raster_processing/derivatives/parameters.py

@@ -11,9 +11,10 @@
     NonSquarePixelSizeException,
 )
 from eis_toolkit.raster_processing.derivatives.partial_derivatives import coefficients
-from eis_toolkit.raster_processing.derivatives.utilities import reduce_ndim, scale_raster, set_flat_pixels
+from eis_toolkit.raster_processing.derivatives.utilities import scale_raster, set_flat_pixels
 from eis_toolkit.utilities.checks.raster import check_quadratic_pixels
 from eis_toolkit.utilities.conversions import convert_rad_to_deg, convert_rad_to_rise
+from eis_toolkit.utilities.miscellaneous import reduce_ndim
 from eis_toolkit.utilities.nodata import nan_to_nodata, nodata_to_nan
 
 

eis_toolkit/raster_processing/derivatives/partial_derivatives.py

@@ -166,10 +166,9 @@ def coefficients(
     """Calculate the partial derivatives of a given surface.
 
     Args:
-        raster: Input raster.
+        in_array: Input array.
+        cellsize: Cellsize of the input raster.
         method: Method to calculate the partial derivatives.
-        coefficients: Coefficients to calculate from least squares equation.
-        scaling_factor: Factor to modify input raster values, e.g. for unit conversion.
 
     Returns:
         Tuple of the partial derivatives p, q, r, s, t.

eis_toolkit/raster_processing/derivatives/utilities.py

@@ -5,22 +5,6 @@
 from beartype.typing import Union
 
 
-@beartype
-def reduce_ndim(
-    data: np.ndarray,
-) -> np.ndarray:
-    """
-    Reduce the number of dimensions of a numpy array.
-
-    Args:
-        data: The input raster data as a numpy array.
-
-    Returns:
-        The reduced array.
-    """
-    return np.squeeze(data) if data.ndim >= 3 else data
-
-
 @beartype
 def scale_raster(
     data: np.ndarray,
@@ -49,9 +33,9 @@ def set_flat_pixels(
     """Treating values below a certain gradient as flat surface.
 
     Args:
-        data (np.ndarray): Input surface attribute.
-        slope_slope (np.ndarray): Input slope array in degrees.
-        slope_tolerance (Number): Value below a surface will be treated as flat surface (degrees).
+        in_array: Input array.
+        slope_gradient: Input slope array in degrees or tuple of partial derivatives p and q.
+        slope_tolerance: Value below a surface will be treated as flat surface (degrees).
         parameter: The surface attribute to modify.
 
     Returns:

eis_toolkit/raster_processing/filters/focal.py

@@ -0,0 +1,158 @@
+from numbers import Number
+
+import numpy as np
+import rasterio
+from beartype import beartype
+from beartype.typing import Literal, Optional
+
+from eis_toolkit.raster_processing.filters.kernels import basic_kernel, gaussian_kernel, mexican_hat_kernel
+from eis_toolkit.raster_processing.filters.utilities import apply_correlated_filter, apply_generic_filter, check_inputs
+from eis_toolkit.utilities.miscellaneous import cast_array_to_float, reduce_ndim
+from eis_toolkit.utilities.nodata import nan_to_nodata, nodata_to_nan
+
+
+@beartype
+def _focal_median(window: np.ndarray) -> Number:
+    weighted_value = np.nanmedian(window) if sum(np.isnan(window)) != len(window) else np.nan
+    return weighted_value
+
+
+@beartype
+def focal_filter(
+    raster: rasterio.io.DatasetReader,
+    method: Literal["mean", "median"] = "mean",
+    size: int = 3,
+    shape: Literal["square", "circle"] = "circle",
+) -> tuple[np.ndarray, dict]:
+    """
+    Apply a basic focal filter to the input raster.
+
+    Args:
+        raster: The input raster dataset.
+        method: The method to use for filtering. Can be either "mean" or "median". Default to "mean".
+        size: The size of the filter window. E.g., 3 means a 3x3 window. Default to 3.
+        shape: The shape of the filter window. Can be either "square" or "circle". Default to "circle".
+
+    Returns:
+        The filtered raster array.
+
+    Raises:
+        InvalidRasterBandException: If the input raster has more than one band.
+        InvalidParameterValueException: If the filter size is smaller than 3.
+            If the filter size is not an odd number.
+            If the shape is not "square" or "circle".
+    """
+    check_inputs(raster, size)
+
+    kernel = basic_kernel(size, shape)
+
+    raster_array = raster.read()
+    raster_array = reduce_ndim(raster_array)
+    raster_array = nodata_to_nan(raster_array, raster.nodata)
+
+    if method == "mean":
+        out_array = apply_correlated_filter(raster_array, kernel)
+    elif method == "median":
+        out_array = apply_generic_filter(raster_array, _focal_median, kernel)
+
+    out_array = nan_to_nodata(out_array, raster.nodata)
+    out_array = cast_array_to_float(out_array, cast_float=True)
+    out_meta = raster.meta.copy()
+
+    return out_array, out_meta
+
+
+@beartype
+def gaussian_filter(
+    raster: rasterio.io.DatasetReader,
+    sigma: Number = 1,
+    truncate: Number = 4,
+    size: Optional[int] = None,
+) -> tuple[np.ndarray, dict]:
+    """
+    Apply a gaussian filter to the input raster.
+
+    Args:
+        raster: The input raster dataset.
+        sigma: The standard deviation of the gaussian kernel.
+        truncate: The truncation factor for the gaussian kernel based on the sigma value.
+            Only if size is not given. Default to 4.0.
+            E.g., for sigma = 1 and truncate = 4.0, the kernel size is 9x9.
+        size: The size of the filter window. E.g., 3 means a 3x3 window.
+            If size is not None, it overrides the dynamic size calculation based on sigma and truncate.
+            Default to None.
+
+    Returns:
+        The filtered raster array.
+
+    Raises:
+        InvalidRasterBandException: If the input raster has more than one band.
+        InvalidParameterValueException: If the filter size is smaller than 3.
+            If the filter size is not an odd number.
+            If the resulting radius is smaller than 1.
+    """
+    check_inputs(raster, size, sigma, truncate)
+
+    kernel = gaussian_kernel(sigma, truncate, size)
+
+    raster_array = raster.read()
+    raster_array = reduce_ndim(raster_array)
+    raster_array = nodata_to_nan(raster_array, raster.nodata)
+
+    out_array = apply_correlated_filter(raster_array, kernel)
+    out_array = nan_to_nodata(out_array, raster.nodata)
+
+    out_array = cast_array_to_float(out_array, cast_float=True)
+    out_meta = raster.meta.copy()
+
+    return out_array, out_meta
+
+
+@beartype
+def mexican_hat_filter(
+    raster: rasterio.io.DatasetReader,
+    sigma: Number = 1,
+    truncate: Number = 4,
+    size: Optional[int] = None,
+    direction: Literal["rectangular", "circular"] = "circular",
+) -> tuple[np.ndarray, dict]:
+    """
+    Apply a mexican hat filter to the input raster.
+
+    Circular: Lowpass filter for smoothing.
+    Rectangular: Highpass filter for edge detection. Results may need further normalization.
+
+    Args:
+        raster: The input raster dataset.
+        sigma: The standard deviation.
+        truncate: The truncation factor.
+            E.g., for sigma = 1 and truncate = 4.0, the kernel size is 9x9.
+            Default to 4.0.
+        size: The size of the filter window. E.g., 3 means a 3x3 window. Default to None.
+        direction: The direction of calculating the kernel values.
+            Can be either "rectangular" or "circular". Default to "circular".
+
+    Returns:
+       The filtered raster array.
+
+    Raises:
+        InvalidRasterBandException: If the input raster has more than one band.
+        InvalidParameterValueException: If the filter size is smaller than 3.
+            If the filter size is not an odd number.
+            If the resulting radius is smaller than 1.
+    """
+    check_inputs(raster, size, sigma, truncate)
+
+    kernel = mexican_hat_kernel(sigma, truncate, size, direction)
+
+    raster_array = raster.read()
+    raster_array = reduce_ndim(raster_array)
+    raster_array = nodata_to_nan(raster_array, raster.nodata)
+
+    out_array = apply_correlated_filter(raster_array, kernel)
+    out_array = nan_to_nodata(out_array, raster.nodata)
+
+    out_array = cast_array_to_float(out_array, cast_float=True)
+    out_meta = raster.meta.copy()
+
+    return out_array, out_meta

eis_toolkit/raster_processing/filters/kernels.py

@@ -0,0 +1,125 @@
+from numbers import Number
+
+import numpy as np
+from beartype import beartype
+from beartype.typing import Literal, Optional
+from scipy.signal import ricker
+
+
+@beartype
+def _get_kernel_size(sigma: Number, truncate: Number, size: Optional[int]) -> tuple[int, int]:
+    """
+    Calculate the kernel size and radius based on the given parameters.
+
+    Args:
+        sigma: The standard deviation.
+        truncate: The truncation factor.
+        size: The optional size of the kernel.
+
+    Returns:
+        A tuple containing the calculated size and radius of the kernel.
+    """
+    if size is not None:
+        radius = int(size // 2)
+    else:
+        radius = int(float(truncate) * float(sigma) + 0.5)
+        size = int(2 * radius + 1)
+
+    return size, radius
+
+
+@beartype
+def _create_grid(radius: int, size: int) -> tuple[np.ndarray, np.ndarray]:
+    """
+    Create a grid of coordinates.
+
+    Args:
+        radius: The radius of the grid.
+        size: The size of the grid.
+
+    Returns:
+        A tuple with x and y coordinates of the grid.
+    """
+    y, x = np.ogrid[-radius : (size - radius), -radius : (size - radius)]  # noqa: E203
+    return x, y
+
+
+@beartype
+def basic_kernel(size: int, shape: Literal["square", "circle"]) -> np.ndarray:
+    """
+    Generate a basic kernel of a specified size and shape.
+
+    Args:
+        size: The size of the kernel.
+        shape: The shape of the kernel. Can be either "square" or "circle".
+
+    Returns:
+        The generated kernel.
+    """
+    if shape == "square":
+        kernel = np.ones((size, size))
+    elif shape == "circle":
+        radius = int(size // 2)
+        x, y = _create_grid(radius, size)
+        mask = x**2 + y**2 <= radius**2
+        kernel = np.zeros((size, size))
+        kernel[mask] = 1
+
+    return kernel
+
+
+@beartype
+def gaussian_kernel(sigma: Number, truncate: Number, size: Optional[int]) -> np.ndarray:
+    """
+    Generate a Gaussian kernel for image denoising.
+
+    Args:
+        sigma: Standard deviation.
+        truncate: Truncation value.
+        size: Size of the kernel. If not provided, it will be calculated based on sigma and truncate.
+
+    Returns:
+        The Gaussian kernel.
+    """
+    size, radius = _get_kernel_size(sigma, truncate, size)
+
+    x, y = _create_grid(radius, size)
+    kernel = 1 / (2 * np.pi * sigma**2) * np.exp((x**2 + y**2) / (2 * sigma**2) * -1)
+    kernel /= np.max(kernel)
+
+    return kernel
+
+
+@beartype
+def mexican_hat_kernel(
+    sigma: Number, truncate: Number, size: Optional[int], direction: Literal["rectangular", "circular"]
+) -> np.ndarray:
+    """
+    Generate a Mexican Hat kernel for denoising (circular) or edge detection (rectangular).
+
+    Args:
+        sigma: The standard deviation of the Gaussian function.
+        truncate: The truncation factor for the Gaussian function.
+        size: The size of the kernel.
+        direction: Shape of the value distribution, either "rectangular" or "circular".
+
+    Returns:
+        The Mexican Hat kernel.
+    """
+    size, radius = _get_kernel_size(sigma, truncate, size)
+
+    if direction == "rectangular":
+        ricker_wavelet = ricker(size, sigma)
+        kernel = np.outer(ricker_wavelet, ricker_wavelet)
+    elif direction == "circular":
+        x, y = _create_grid(radius, size)
+        kernel = (
+            2
+            / np.sqrt(3 * sigma)
+            * np.pi ** (1 / 4)
+            * (1 - (x**2 + y**2) / sigma**2)
+            * np.exp(-(x**2 + y**2) / (2 * sigma**2))
+        )
+    kernel /= np.max(kernel)
+
+    return kernel