Buffer road lines in the raster's own CRS

combine_vector_targets reprojected everything to EPSG:3857, buffered the
line features there, then reprojected back to the raster CRS. The detour
existed to get a CRS whose units are metres - but the raster's own CRS
normally already is one, since Sentinel-2 is UTM.

Reprojecting straight to it does in one pass what took two. Measured on a
60 km Perth box carrying 918,399 Overture features (256,906 roads and
661,493 buildings), the function's peak drops from +2.67 GB to +1.26 GB;
the two reprojections were 2.1 GB of the 2.7 GB it spent. A full tile
covers ~3.4x the area.

It also buffers by the distance actually asked for. EPSG:3857's scale
factor is 1/cos(latitude), so buffer(distance=5) there laid down about
4.24 m on the ground at Perth's latitude, and a different amount at every
other one - the same code gave a Kununurra scene and a Hobart scene
different real-world road widths. Roads are now buffered a true 5 m,
which widens them by ~18% and sets 1.7% more negative-target pixels
(7,470,943 -> 7,598,563 on the box measured).

A geographic raster CRS cannot buffer in metres, so it keeps the detour
through EPSG:3857 and is reprojected back as before.

make_valid is unchanged.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Author: wrignj08

omniwatermask/target_builders.py

@@ -156,8 +156,20 @@ def combine_vector_targets(
     if all_targets.empty:
         return None
 
-    # re-project in 3857 to buffer and then to raster crs
-    all_targets = gpd.GeoDataFrame(all_targets).to_crs("EPSG:3857")
+    # Buffering needs a projected CRS so the distance is in metres, and the
+    # raster's own CRS normally is one - Sentinel-2 is UTM. Reprojecting
+    # straight to it does in one pass what used to take two, a round trip out
+    # to EPSG:3857 and back, which cost 2.1 GB of the 2.7 GB this function
+    # spent on a Perth tile's 918k features.
+    #
+    # It also buffers by the distance actually asked for. EPSG:3857's scale
+    # factor is 1/cos(latitude), so buffering 5 m there lays down ~4.2 m on the
+    # ground at Perth, and a different amount at every other latitude. Only a
+    # geographic raster CRS still needs the detour.
+    raster_crs = CRS.from_user_input(raster_src.crs)
+    buffer_crs = raster_crs if raster_crs.is_projected else CRS.from_epsg(3857)
+
+    all_targets = gpd.GeoDataFrame(all_targets).to_crs(buffer_crs)
     if all_targets is None:
         return None
     # remove points
@@ -174,7 +186,8 @@ def combine_vector_targets(
 
     all_targets["geometry"] = all_targets.geometry.make_valid()
 
-    all_targets = gpd.GeoDataFrame(all_targets).to_crs(raster_src.crs)
+    if buffer_crs != raster_crs:
+        all_targets = gpd.GeoDataFrame(all_targets).to_crs(raster_crs)
 
     return gpd.GeoDataFrame(all_targets)
 

tests/test_target_builders.py

@@ -1,3 +1,6 @@
+from rasterio.transform import from_bounds
+import rasterio as rio
+import numpy as np
 import logging
 from queue import Queue
 from unittest.mock import patch
@@ -582,3 +585,67 @@ def test_records_source_and_ocean_in_cache_db(
         assert len(df) == 1
         assert df["source"].iloc[0] == "overture"
         assert bool(df["ocean"].iloc[0]) is True
+
+
+class TestCombineVectorTargetsBufferCrs:
+    """Lines are buffered in a projected CRS so the distance means metres."""
+
+    @staticmethod
+    def _raster(tmp_path, crs, transform, name):
+        path = tmp_path / name
+        with rio.open(
+            path,
+            "w",
+            driver="GTiff",
+            height=100,
+            width=100,
+            count=1,
+            dtype="uint8",
+            crs=crs,
+            transform=transform,
+        ) as dst:
+            dst.write(np.zeros((1, 100, 100), "uint8"))
+        return rio.open(path)
+
+    def test_projected_raster_buffers_by_a_true_five_metres(self, tmp_path):
+        """EPSG:3857 would lay down ~4.2 m here, its scale factor at -32 lat.
+
+        Buffering in the raster's own UTM avoids both that error and the
+        reprojection round trip it used to require.
+        """
+        src = self._raster(
+            tmp_path,
+            "EPSG:32750",
+            from_bounds(390000, 6460000, 391000, 6461000, 100, 100),
+            "utm.tif",
+        )
+        line = gpd.GeoDataFrame(
+            geometry=[LineString([(390200, 6460500), (390800, 6460500)])],
+            crs="EPSG:32750",
+        ).to_crs("EPSG:4326")
+
+        out = combine_vector_targets([line], src)
+        assert out.crs == src.crs
+        minx, miny, maxx, maxy = out.total_bounds
+        # a 5 m buffer either side of a horizontal line: 10 m tall
+        assert (maxy - miny) == pytest.approx(10.0, abs=0.2)
+
+    def test_geographic_raster_still_buffers_in_metres(self, tmp_path):
+        """A degrees-based raster CRS cannot buffer in metres, so it detours."""
+        src = self._raster(
+            tmp_path,
+            "EPSG:4326",
+            from_bounds(115.80, -32.00, 115.90, -31.90, 100, 100),
+            "geo.tif",
+        )
+        line = gpd.GeoDataFrame(
+            geometry=[LineString([(115.82, -31.95), (115.88, -31.95)])],
+            crs="EPSG:4326",
+        )
+
+        out = combine_vector_targets([line], src)
+        assert out.crs == src.crs
+        assert (out.geometry.type == "Polygon").all()
+        # a few metres expressed in degrees of latitude is small but non-zero
+        minx, miny, maxx, maxy = out.total_bounds
+        assert 0 < (maxy - miny) < 0.01