feat: add circular raster buffer

RELEASE_NOTES.rst

@@ -17,6 +17,13 @@ Upcoming Release
    To use the features already you have to install the ``master`` branch, e.g.
    ``pip install git+https://github.com/pypsa/atlite``.
 
+**Features**
+
+* Add ``buffer_geometry`` argument to ``ExclusionContainer.add_raster``. Choose
+  between the historic ``"diamond"`` buffer (default) and a new ``"circular"``
+  buffer that computes a geometrically accurate Euclidean buffer, correct in
+  diagonal directions.
+
 **Bug fixes**
 
 * Fix ``Cutout.line_rating`` passing line azimuth in radians while

atlite/gis.py

@@ -27,6 +27,7 @@
 from rasterio.plot import show
 from rasterio.warp import transform_bounds
 from scipy.ndimage import binary_dilation as dilation
+from scipy.ndimage import distance_transform_edt
 from shapely.ops import transform
 from shapely.strtree import STRtree
 from tqdm import tqdm
@@ -430,6 +431,10 @@ def shape_availability(
     transform : rasterio.Affine
         Affine transform of the mask.
 
+    Raises
+    ------
+    ValueError
+        If an unsupported buffer geometry is encountered.
     """
     if not excluder.all_open:
         excluder.open_files()
@@ -462,8 +467,23 @@ def shape_availability(
         if d["invert"]:
             masked_ = ~masked_
         if d["buffer"]:
-            iterations = int(d["buffer"] / excluder.res) + 1
-            masked_ = dilation(masked_, iterations=iterations)
+            buffer_geom = d.get("buffer_geometry", "diamond")
+            if buffer_geom == "circular":
+                if masked_.any():
+                    radius = d["buffer"] / excluder.res
+                    masked_ = distance_transform_edt(~masked_) <= radius
+                else:
+                    # If mask is empty, distance_transform_edt treats the boundary
+                    # as background, which would incorrectly exclude border pixels.
+                    pass
+            elif buffer_geom == "diamond":
+                iterations = int(d["buffer"] / excluder.res) + 1
+                masked_ = dilation(masked_, iterations=iterations)
+            else:
+                raise ValueError(
+                    f"Unsupported buffer geometry: '{buffer_geom}'. "
+                    "Must be one of 'diamond' or 'circular'."
+                )
 
         exclusions = exclusions | masked_
 
@@ -557,6 +577,7 @@ def add_raster(
         | Callable[[NDArray], NDArray]
         | None = None,
         buffer: float = 0,
+        buffer_geometry: str = "diamond",
         invert: bool = False,
         nodata: int = 255,
         allow_no_overlap: bool = False,
@@ -579,6 +600,21 @@ def add_raster(
             Buffer around the excluded areas in units of ExclusionContainer.crs.
             Use this to create a buffer around the excluded/included area.
             The default is 0.
+        buffer_geometry : {"diamond", "circular"}, optional
+            Shape of the buffer applied around raster exclusions. Only used
+            when ``buffer > 0``. The default is ``"diamond"``.
+
+            * ``"diamond"`` – dilates the mask with
+              :func:`scipy.ndimage.binary_dilation`, expanding it by one pixel
+              in the four cardinal directions per iteration. The footprint is a
+              diamond (L1 metric), so the buffer distance is only accurate along
+              the grid axes and overreaches diagonally. This is the historic
+              behaviour.
+            * ``"circular"`` – marks every pixel whose exact Euclidean distance
+              to the mask (via :func:`scipy.ndimage.distance_transform_edt`) is
+              within ``buffer``, giving a geometrically accurate circular
+              (L2 metric) buffer. Prefer this when diagonal accuracy matters;
+              note that in practice it is typically slower than ``"diamond"``.
         nodata : int, optional
             Value to use for nodata pixels. The default is 255.
         invert : bool, optional
@@ -591,11 +627,22 @@ def add_raster(
         crs : rasterio.CRS/EPSG
             CRS of the raster. Specify this if the raster has invalid crs.
 
+        Raises
+        ------
+        ValueError
+            If ``buffer_geometry`` is not one of 'diamond' or 'circular'.
         """
+        if buffer_geometry not in ("diamond", "circular"):
+            raise ValueError(
+                f"Invalid buffer_geometry: '{buffer_geometry}'. "
+                "Must be one of 'diamond' or 'circular'."
+            )
+
         d: dict[str, Any] = {
             "raster": raster,
             "codes": codes,
             "buffer": buffer,
+            "buffer_geometry": buffer_geometry,
             "invert": invert,
             "nodata": nodata,
             "allow_no_overlap": allow_no_overlap,

test/test_gis.py

@@ -12,6 +12,7 @@
 # IDEAS for tests
 
 import functools
+import pickle
 import warnings
 
 import geopandas as gpd
@@ -646,3 +647,60 @@ def test_plot_shape_availability(ref, raster):
     with warnings.catch_warnings():
         warnings.simplefilter("ignore")
         excluder.plot_shape_availability(shapes)
+
+
+class TestRasterBufferGeometry:
+    """Diamond vs. circular raster buffering: validation, serialization, edge cases."""
+
+    res = 0.01
+
+    @pytest.fixture
+    def shapes(self, ref):
+        return gpd.GeoSeries([box(X0, Y0, X1, Y1)], crs=ref.crs)
+
+    def test_serialization_roundtrip_preserves_geometry(self, shapes, raster):
+        excluder = ExclusionContainer(shapes.crs, res=self.res)
+        excluder.add_raster(raster, buffer=self.res, buffer_geometry="circular")
+
+        restored = pickle.loads(pickle.dumps(excluder))
+        assert restored.rasters[0]["buffer_geometry"] == "circular"
+
+    def test_missing_geometry_key_defaults_to_diamond(self, shapes, raster):
+        legacy = ExclusionContainer(shapes.crs, res=self.res)
+        legacy.add_raster(raster, buffer=self.res, buffer_geometry="circular")
+        del legacy.rasters[0]["buffer_geometry"]
+
+        diamond = ExclusionContainer(shapes.crs, res=self.res)
+        diamond.add_raster(raster, buffer=self.res, buffer_geometry="diamond")
+
+        masked_legacy, _ = shape_availability(shapes, legacy)
+        masked_diamond, _ = shape_availability(shapes, diamond)
+        assert (masked_legacy == masked_diamond).all()
+
+    def test_circular_buffer_on_empty_mask_keeps_borders(self, shapes, raster):
+        excluder = ExclusionContainer(shapes.crs, res=self.res)
+        excluder.add_raster(
+            raster, codes=[-1], buffer=self.res * 5, buffer_geometry="circular"
+        )
+        masked, _ = shape_availability(shapes, excluder)
+        assert masked.all()
+
+    def test_invalid_geometry_rejected_at_registration(self, shapes, raster):
+        excluder = ExclusionContainer(shapes.crs, res=self.res)
+        with pytest.raises(ValueError, match="Invalid buffer_geometry"):
+            excluder.add_raster(raster, buffer=self.res, buffer_geometry="hexagon")
+
+    def test_invalid_geometry_rejected_at_computation(self, shapes, raster):
+        excluder = ExclusionContainer(shapes.crs, res=self.res)
+        excluder.rasters.append({
+            "raster": raster,
+            "codes": None,
+            "buffer": self.res,
+            "buffer_geometry": "hexagon",
+            "invert": False,
+            "nodata": 255,
+            "allow_no_overlap": False,
+            "crs": None,
+        })
+        with pytest.raises(ValueError, match="Unsupported buffer geometry"):
+            shape_availability(shapes, excluder)