feat: add circular raster buffer

[jome1]

Changes proposed in this Pull Request

Add buffer_geometry argument to ExclusionContainer.add_raster()
Previously, the buffer applied to raster exclusions in shape_availability() always used scipy.ndimage.binary_dilation which produces a diamond-shaped (L1-metric) buffer: each iteration expands the excluded area by one pixel in the four cardinal directions, so the footprint after N iterations is a diamond of radius N.
This PR adds a buffer_geometry parameter to add_raster() that lets users choose between two buffer shapes:
• "diamond" (default) – retains the existing binary_dilation behaviour. No change to existing code that does not pass buffer_geometry.
• "circular" – uses scipy.ndimage.distance_transform_edt to compute exact Euclidean distances. Every pixel whose Euclidean distance to the nearest excluded pixel is ≤ buffer / res is marked as excluded, producing a geometrically accurate circular buffer.

According to Claude, the circular buffer method migth also be faster:
"Performance: The circular method runs in O(pixels) with a single linear-time scan regardless of buffer size. The diamond method runs in O(N × pixels) where N = buffer / res, meaning it becomes proportionally slower as the buffer grows. For typical use cases (e.g. a 500 m buffer at 100 m resolution, N = 5) the circular method is already several times faster, and the gap widens with larger buffers."

Additional info

I got aware of the raster buffering while testing atlite raster exclusion with buffering. I measured the distance of the buffer in the output raster and the buffer distance was only correct along the grid (vertically and horizontally) but not in diagonal directions.
I used Claude to understand the current raster buffering and to implement the circular buffering.
Before adding realease notes and add the argument to add_raster() in the documentation, please give the PR a review and let me know if it makes sense.

I used the following code in a Jupyter Notebook (separate code cells) to test the results (Let me know if you need the input files to test it)

import atlite
import os
import geopandas as gpd
import matplotlib.pyplot as plt

regionPath = os.path.join("Odense_square_EPSG4326.geojson")
region = gpd.read_file(regionPath)
region = region.to_crs(25832) #local UTM zone
print(region.crs)

# initiate Exclusion container
excluder = atlite.ExclusionContainer(crs=25832, res=10)
# add raster (single pixel) with buffer
excluder.add_raster("single_pixel.tif", codes=1, buffer=500, buffer_geometry="diamond", crs=4326)
masked, transform = excluder.compute_shape_availability(region)
# plot
fig, ax = plt.subplots()
excluder.plot_shape_availability(region)

# initiate Exclusion container
excluder = atlite.ExclusionContainer(crs=25832, res=10)
# add raster (single pixel) with buffer
excluder.add_raster("single_pixel.tif", codes=1, buffer=500, buffer_geometry="circular", crs=4326)
masked, transform = excluder.compute_shape_availability(region)
# plot
fig, ax = plt.subplots()
excluder.plot_shape_availability(region)

These are the plots with the different raster buffering methods (there is a single pixel in the center. The buffer is applied around this pixel):

Checklist

• Code changes are sufficiently documented; i.e. new functions contain docstrings and further explanations may be given in doc.
• Unit tests for new features were added (if applicable).
• Newly introduced dependencies are added to environment.yaml, environment_docs.yaml and setup.py (if applicable).
• A note for the release notes doc/release_notes.rst of the upcoming release is included.
• I consent to the release of this PR's code under the MIT license.

[maurerle]

Great idea. This looks fine in general.

I made a suggestion to improve this here:
jome1#1

and also added unit tests to verify the suggested improvements.

Regarding the performance claim, it looks as if circular would be faster, because it only does one (heavy) iteration, instead of multiple iterations.
A small test did show that diamond is still much faster:

Generating 5000x5000 mask with a diagonal exclusion line...
Buffer Size (px)   | Diamond Time (s)   | Circular Time (s)  | Speedup Ratio
---------------------------------------------------------------------------
20                 | 0.32025            | 2.02672            | 0.16
100                | 0.33293            | 2.01348            | 0.17
1000               | 0.50594            | 2.00961            | 0.25
10000              | 0.83164            | 1.98676            | 0.42
100000             | 0.82718            | 2.03129            | 0.41
1000000            | 0.82630            | 2.02098            | 0.41

Just as an info aside

benchmark source code

import time
import numpy as np
import scipy.ndimage

def run_benchmark():
    # Large 5000 x 5000 grid (25 million pixels)
    shape = (5000, 5000)
    print(f"Generating {shape[0]}x{shape[1]} mask with a diagonal exclusion line...")

    mask = np.zeros(shape, dtype=bool)
    # Add a diagonal line representing a linear exclusion (e.g. road, river, or coastline)
    for i in range(shape[0]):
        mask[i, i] = True

    buffer_sizes = [20, 100, 1000, 10000, 100000, 1000000]

    print(f"{'Buffer Size (px)':<18} | {'Diamond Time (s)':<18} | {'Circular Time (s)':<18} | {'Speedup Ratio':<15}")
    print("-" * 75)

    for N in buffer_sizes:
        # Benchmark Diamond Buffering (L1 Metric / binary dilation)
        t0 = time.perf_counter()
        res_diamond = scipy.ndimage.binary_dilation(mask, iterations=N)
        t_diamond = time.perf_counter() - t0

        # Benchmark Circular Buffering (L2 Metric / EDT)
        t0 = time.perf_counter()
        if mask.any():
            res_circular = scipy.ndimage.distance_transform_edt(~mask) <= N
        else:
            res_circular = mask
        t_circular = time.perf_counter() - t0

        speedup = t_diamond / t_circular
        print(f"{N:<18} | {t_diamond:<18.5f} | {t_circular:<18.5f} | {speedup:<15.2f}x")

if __name__ == "__main__":
    run_benchmark()

[maurerle]

A base case is missing here, instead of silently doing nothing.

Also importing old buffers would not have an entry for buffer_geometry, and do nothing.
This would break importing older data. Adding a default to "diamond", if nothing is set would be necessary

[jome1]

Ah yes, very good!
"diamond" is also the default in the add_raster() function, so old code does not need to be updated with the buffer_geometry argument.

[jome1]

This is resolved in your PR with the raising the ValueError right? Or do you want me to adjust something?

[jome1]

Amazing! Thanks for you swift review and the improvements @maurerle!
I merged your PR.

Great that you did an actual speed test. Just too bad that the "circular" buffer is much slower. I am running a big resource potential study for Kenya at 10m resolution. I have not test the "circular" buffer yet, but so far the runs always needed around 2h. I want to use the "circular" buffer because it makes more sense, but then it will take even longer...

[maurerle]

@jome1 I tested on a different computer where the differences between the two was much lower.
YMMV so I'd recommend you to run the provided benchmark script (in my comment above) on your hardware.