Filters

Filters#

GeoUtils implements filtering capabilities for both rasters and point clouds .

All support scalable execution on large datasets and geospatial awareness to define filter distances and skip invalid values.

Filters are not limited to smoothing data or removing outliers (e.g. Gaussian or median filters). They can perform any computation on a window (for rasters) or neighborhood (for points), such as deriving a local derivative or a local autocorrelation index.

Raster filters#

ds.rst.filter() or Raster.filter()

A raster filter computes each output pixel value from nearby input pixels, returning a new raster with the same shape and metadata. Neighborhoods are defined with a window size and shape, and optionally support fractional cell coverage.

For example, apply a median filter in a 15 × 15 pixel window:

import geoutils as gu
filename_rast = gu.examples.get_path("exploradores_aster_dem")
ds = gu.open_raster(filename_rast)
# Filter the raster with a median filter of size 15
ds_filtered = ds.rst.filter(method="median", backend="scipy", window=15)

_images/119cc98fe2a54e3c0b8d29bc165ef7197a37c4842f0d898c065115890ad889ea.png

Use window to set the window size in pixels (odd, three by default), sigma sets the Gaussian standard deviation in pixels. Filters use SciPy by default, use backend="numba" to use Numba.

Point filters#

points.pc.filter() or PointCloud.filter()

A point filter computes each output point value from nearby input values, returning a new point cloud with the same coordinates and metadata. Neighborhoods are defined by a radius in coordinate units and maximum closest number of points.

For example, to compute the mean in a 1000-m radius, using at most the closest 20 points:

# Sample the projected DEM, then average up to 20 measurements within 1000 m
points = ds.rst.subsample(1500, random_state=42)
smoothed_points = points.pc.filter(
    method="mean", backend="scipy", radius=1000, max_points=20
)

_images/50f96edb5dd655c72269a912a57aaa35ff7f32634a7a5299559e4a66f2e1ba2c.png

Use include_self=False to exclude the point original value at this coordinate, or min_points to require enough valid neighbors. Lower batch_size to reduce temporary memory for dense neighborhoods.

Available methods#

GeoUtils supports the following raster and point filters. See Reducers for illustrations of supported reducers.

Customizing filters#

Nodata propagation#

Filters can customize how they handle missing values. By default, we ignore them, but they can also be propagated.

The following example covers both behaviors by applying a mean filter to a DEM with gaps:

Custom point reducers#

You can define a custom reducer to customize the calculation on neighboring values, while GeoUtils handles their selection and missing data.

Below, we show how define a custom trimmed mean, which discards the lowest and highest values before averaging. We apply it to a point cloud with a few added outliers to illustrate how it reduces their influence:

from geoutils.operators import LocalData, Reducer


class TrimmedMean(Reducer):
    """Custom trimmed mean reducer."""

    # We define the minimum values required, below that a NaN is returned
    minimum_inputs = 3

    # The custom reducer simply needs to implement "reduce", which receives a local data object
    def reduce(self, data: LocalData) -> float:
        # Remove the two extreme values
        values = data.values
        trimmed_sum = np.sum(values) - np.min(values) - np.max(values)
        return float(trimmed_sum / (values.size - 2))


# Add some outliers to the input raster
noisy_values = np.array(points.pc.data, dtype=float, copy=True)
noisy_values[::100] += 2000
noisy_points = points.pc.copy(new_array=noisy_values)
trimmed_points = noisy_points.pc.filter(
    TrimmedMean(), radius=1000, max_points=20, batch_size=256
)

_images/53ade693f044ea0ec2abe5c5b2f5adebf36d876ab8a3ea271681edd850427bd7.png

Custom raster functions#

Raster filters also accept custom functions operating on a NumPy array. Below, we derive terrain slope by combining the horizontal and vertical elevation gradients:

import numpy as np

# Convert the elevation gradient in each map direction to slope in degrees
pixel_width, pixel_height = ds.rst.res


def slope_filter(arr: np.ndarray) -> np.ndarray:
    """Calculate terrain slope in degrees from the elevation gradients."""
    gradient_y, gradient_x = np.gradient(arr, pixel_height, pixel_width)
    return np.degrees(np.arctan(np.hypot(gradient_x, gradient_y)))


ds_slope = ds.rst.filter(slope_filter, window=3)

_images/c0f51a74c9d39d1f90ecbb22d244b87308cc0f5621403f12badeb887282114b3.png

Method	Description	Raster	Point
`count`	Number of valid neighbors	✓	✓
`sum`	Sum of neighboring values	✓	✓
`mean`, `average`	Mean of neighboring values	✓	✓
`min`, `minimum`	Minimum value	✓	✓
`max`, `maximum`	Maximum value	✓	✓
`range`	Difference between maximum and minimum	✓	✓
`rms`, `rmse`	Root mean square of neighboring values	✓	✓
`stdev`, `std`, `standard_deviation`	Population standard deviation	✓	✓
`median`, `med`	Median of neighboring values	✓	✓
`mode`	Most frequent value	✓	✓
`q1`, `q3`	First or third quartile	✓	✓
`quantile`	Quantile at the supplied fraction	✓	✓
`average_distance`	Mean distance from neighbors to the target	✓	✓
`average_distance_pts`	Mean distance between pairs of neighbors	✓	✓
`gaussian`	Smooth with a Gaussian kernel	✓	—
`distance`	Mask outliers relative to a Gaussian-smoothed average	✓	—
`Reducer`	Custom calculation on each neighborhood	✓	✓
Callable	Custom function taking and returning a 2D NumPy array	✓	—