Split linestrings within bounds (behind "bounded" flag)

Author: tomalrussell

extension/src/intersections.cpp

@@ -53,14 +53,15 @@ std::vector<py::object> convert_cpp2py(std::vector<linestr> splits) {
 
 std::vector<py::object> splitLineString(py::object linestring_py, int nrows,
                                         int ncols,
-                                        std::vector<double> transform) {
+                                        std::vector<double> transform,
+                                        bool bounded = false) {
   linestr linestring = convert_py2cpp(linestring_py);
   transform::Affine affine(transform[0], transform[1], transform[2],
                            transform[3], transform[4], transform[5]);
   grid::Grid grid(ncols, nrows, affine);
   geometry::LineString line(linestring);
   std::vector<linestr> splits =
-      operations::findIntersectionsLineString(line, grid);
+      operations::findIntersectionsLineString(line, grid, bounded);
   return convert_cpp2py(splits);
 }
 
@@ -129,6 +130,9 @@ PYBIND11_MODULE(intersections, m) {
   m.doc() = "Vector geometry to grid intersections";
 
   m.def("split_linestring", &snail::splitLineString,
+        pybind11::arg("linestring_py"), pybind11::arg("nrows"),
+        pybind11::arg("ncols"), pybind11::arg("transform"),
+        pybind11::arg("bounded") = false,
         "Split LineString along a grid");
   m.def("get_cell_indices", &snail::get_cell_indices,
         "Get LineString cell indices in a grid");

extension/src/operations.cpp

@@ -10,6 +10,52 @@ namespace operations {
 
 using linestr = std::vector<geometry::Coord>;
 
+bool segmentIntersectsGridBounds(const geometry::Line &line,
+                                 const grid::Grid &raster) {
+
+  auto xy0 = line.start;
+  auto xy1 = line.end;
+  const double x0 = xy0.x;
+  const double y0 = xy0.y;
+  const double x1 = xy1.x;
+  const double y1 = xy1.y;
+
+  auto xymin = raster.grid_to_world * geometry::Coord(0.0, 0.0);
+  auto xymax =
+      raster.grid_to_world * geometry::Coord(raster.ncols, raster.nrows);
+  const double xmin = xymin.x;
+  const double ymin = xymin.y;
+  const double xmax = xymax.x;
+  const double ymax = xymax.y;
+
+  auto inside = [&](const double &x, const double &y) {
+    return x >= xmin && x <= xmax && y >= ymin && y <= ymax;
+  };
+
+  if (inside(x0, y0) || inside(x1, y1)) {
+    return true;
+  }
+
+  auto crosses_x = [&](const double &xtest) {
+    if (x0 <= xtest && xtest <= x1) {
+      const double ycross = y0 + (xtest - x0) * (y1 - y0) / (x1 - x0);
+      return ymin <= ycross && ycross <= ymax;
+    }
+    return false;
+  };
+
+  auto crosses_y = [&](const double &ytest) {
+    if (y0 <= ytest && ytest <= y1) {
+      const double xcross = x0 + (ytest - y0) * (x1 - x0) / (y1 - y0);
+      return xmin <= xcross && xcross <= xmax;
+    }
+    return false;
+  };
+
+  return crosses_x(xmin) || crosses_x(xmax) || crosses_y(ymin) ||
+         crosses_y(ymax);
+}
+
 /// Piecewise decomposition of a linestring according to intersection points
 std::vector<linestr> split_linestr(linestr linestring, linestr intersections) {
   // Add line start point
@@ -28,27 +74,49 @@ std::vector<linestr> split_linestr(linestr linestring, linestr intersections) {
 
 /// Find intersection points of a linestring with a raster grid
 std::vector<linestr>
-findIntersectionsLineString(geometry::LineString linestring,
-                            grid::Grid raster) {
+findIntersectionsLineString(geometry::LineString linestring, grid::Grid raster,
+                            bool bounded) {
   linestr coords = linestring.coordinates;
 
   std::vector<linestr> allsplits;
   linestr linestr_piece;
   for (std::size_t i = 0; i < coords.size() - 1; i++) {
     geometry::Line line(coords.at(i), coords.at(i + 1));
 
-    // If the line starts and ends in different cells, it needs to be cleaned.
-    if (raster.cellIndex(line.start) != raster.cellIndex(line.end)) {
+    bool single_cell =
+        raster.cellIndices(line.start) == raster.cellIndices(line.end);
+
+    // If the line starts and ends in the same cell,
+    // or (bounded and (segment does not intersect overall grid bounds))
+    if (single_cell ||
+        (bounded && !segmentIntersectsGridBounds(line, raster))) {
+      // then don't split, just push back the coordinate
+      linestr_piece.push_back(coords.at(i));
+    } else {
+      // otherwise do split this straight-line segment
       linestr intersections = raster.findIntersections(line);
+      // if only splitting within grid bounds, filter the intersections
+      if (bounded) {
+        linestr filtered;
+        filtered.reserve(intersections.size());
+        for (std::size_t idx = 0; idx < intersections.size(); ++idx) {
+          const auto &coordinate = intersections[idx];
+          bool endpoint = (idx == 0) || (idx == intersections.size() - 1);
+          // keep the endpoints as original coordinates in the linestring
+          // and keep any split intersections from within the grid bounds
+          if (endpoint || pointInBounds(coordinate, raster)) {
+            filtered.push_back(coordinate);
+          }
+        }
+        intersections = std::move(filtered);
+      }
       std::vector<linestr> splits = split_linestr(linestr_piece, intersections);
       allsplits.insert(allsplits.end(), splits.begin(), splits.end());
       if (line.end == intersections.back()) {
         linestr_piece = {};
       } else {
         linestr_piece = {intersections.back()};
       }
-    } else {
-      linestr_piece.push_back(coords.at(i));
     }
   }
 
@@ -57,9 +125,16 @@ findIntersectionsLineString(geometry::LineString linestring,
     allsplits.push_back(linestr_piece);
   }
 
-  return (allsplits);
+  return allsplits;
 }
 
+bool pointInBounds(const geometry::Coord &pt, const grid::Grid &raster) {
+  auto ll = raster.grid_to_world * geometry::Coord(0.0, 0.0);
+  auto ur = raster.grid_to_world * geometry::Coord(raster.ncols, raster.nrows);
+
+  return pt.x >= ll.x && pt.x <= ur.x && pt.y >= ll.y && pt.y <= ur.y;
+};
+
 bool isOnGridLine(geometry::Coord point, Direction direction, double level,
                   double cellSize) {
   switch (direction) {
@@ -81,7 +156,7 @@ bool isOnGridLine(geometry::Coord point, Direction direction, double level,
 //  >>-----x----o-----o-----  (don't include x)
 //        /.\   |..../
 //
-// TODO figure out what to do when some portion of the boundary is already
+// figure out what to do when some portion of the boundary is already
 // along the grid line. This is a legitimate case for odd number of crossings:
 //
 //              |......|

extension/src/operations.hpp

@@ -11,11 +11,14 @@ namespace operations {
 enum class Direction { horizontal, vertical };
 
 std::vector<std::vector<geometry::Coord>>
-    findIntersectionsLineString(geometry::LineString, grid::Grid);
+findIntersectionsLineString(geometry::LineString, grid::Grid,
+                            bool bounded = false);
 std::vector<std::vector<geometry::Coord>>
 splitAlongGridlines(std::vector<geometry::Coord>, int, int, Direction,
                     grid::Grid);
 
+bool pointInBounds(const geometry::Coord &pt, const grid::Grid &raster);
+
 } // namespace operations
 } // namespace snail
 

extension/tests/tests_intersections.cpp

@@ -258,7 +258,7 @@ TEST_CASE("LineString outside grid remains unchanged", "[decomposition]") {
   snail::geometry::LineString line(coordinates);
 
   auto splits =
-      snail::operations::findIntersectionsLineString(line, test_raster);
+      snail::operations::findIntersectionsLineString(line, test_raster, true);
 
   REQUIRE(splits.size() == 1);
   REQUIRE(splits[0].size() == coordinates.size());
@@ -281,7 +281,7 @@ TEST_CASE("LineString partially overlapping grid splits correctly",
   };
 
   auto splits =
-      snail::operations::findIntersectionsLineString(line, test_raster);
+      snail::operations::findIntersectionsLineString(line, test_raster, true);
 
   REQUIRE(splits.size() == expected.size());
   for (std::size_t i = 0; i < expected.size(); ++i) {
@@ -293,6 +293,45 @@ TEST_CASE("LineString partially overlapping grid splits correctly",
   }
 }
 
+TEST_CASE("Bounded splits keep interior vertices within grid extents",
+          "[decomposition]") {
+  snail::grid::Grid test_raster(2, 2, snail::transform::Affine());
+  linestr coordinates = {{-2.0, 0.5}, {1.5, 0.5}};
+  snail::geometry::LineString line(coordinates);
+
+  auto splits =
+      snail::operations::findIntersectionsLineString(line, test_raster, true);
+
+  REQUIRE(splits.size() == 3);
+  for (std::size_t segment_idx = 0; segment_idx < splits.size();
+       ++segment_idx) {
+    const auto &segment = splits[segment_idx];
+    REQUIRE(segment.size() >= 2);
+    for (std::size_t point_idx = 0; point_idx < segment.size(); ++point_idx) {
+      bool is_start = (segment_idx == 0 && point_idx == 0);
+      bool is_end =
+          (segment_idx == splits.size() - 1 && point_idx == segment.size() - 1);
+      if (is_start || is_end) {
+        continue;
+      }
+      REQUIRE(
+          snail::operations::pointInBounds(segment[point_idx], test_raster));
+    }
+  }
+}
+
+TEST_CASE("pointInBounds treats boundary as inside", "[bounds]") {
+  snail::grid::Grid test_raster(2, 2, snail::transform::Affine());
+  snail::geometry::Coord on_edge_x(2.0, 1.0);
+  REQUIRE(snail::operations::pointInBounds(on_edge_x, test_raster));
+  snail::geometry::Coord on_edge_y(1.0, 2.0);
+  REQUIRE(snail::operations::pointInBounds(on_edge_y, test_raster));
+  snail::geometry::Coord inner(0.9, 1.1);
+  REQUIRE(snail::operations::pointInBounds(inner, test_raster));
+  snail::geometry::Coord outer(0.9, 11.1);
+  REQUIRE(!snail::operations::pointInBounds(outer, test_raster));
+}
+
 TEST_CASE("Split with different grid", "[decomposition]") {
   Config case1;
   case1.expected_splits = {

src/snail/intersection.py

@@ -204,7 +204,9 @@ def split_points(
 
 
 def split_linestrings(
-    linestring_features: geopandas.GeoDataFrame, grid: GridDefinition
+    linestring_features: geopandas.GeoDataFrame,
+    grid: GridDefinition,
+    bounded=False,
 ) -> geopandas.GeoDataFrame:
     """Split linestrings along a grid"""
     # TODO check for MultiLineString
@@ -217,21 +219,13 @@ def split_linestrings(
             grid.width,
             grid.height,
             grid.transform,
+            bounded=bounded,
         )
         for j, s in enumerate(geom_splits):
-            # splitting sometimes returns zero-length linestrings on edge of raster
-            # see below for example linestring on eastern (lon=70W) extent of box
-            # (Pdb) geometry.coords.xy
-            # (array('d', [-70.0, -70.0]), array('d', [18.445832920952196, 18.445832920952196]))
-            # this split geometry has: j = raster_width
-            # however j should be in range: 0 <= j < raster_width
-            # as a hacky workaround, drop any splits with length 0
-            # do we need a nudge off a cell boundary somewhere when performing the splits?
-            if s.length != 0:
-                new_row = linestring_features.iloc[i].copy()
-                new_row.geometry = s
-                new_row["split"] = j
-                pieces.append(new_row)
+            new_row = linestring_features.iloc[i].copy()
+            new_row.geometry = s
+            new_row["split"] = j
+            pieces.append(new_row)
     logger.info(f"Split {len(linestring_features)} edges into {len(pieces)} pieces")
     splits_df = geopandas.GeoDataFrame(pieces, crs=grid.crs, geometry="geometry")
     return splits_df

tests/test_intersection.py

@@ -279,19 +279,24 @@ def test_split_linestrings_outside_grid_returns_geometry(grid):
     outside = gpd.GeoDataFrame(
         {"id": [1]}, geometry=[LineString([(5.0, 0.5), (16.0, 1.5)])]
     )
-    splits = split_linestrings(outside, grid)
+    # must pass bounded=True for non-splitting behaviour
+    splits = split_linestrings(outside, grid, bounded=True)
     assert len(splits) == 1
     assert splits.geometry.iloc[0].equals(outside.geometry.iloc[0])
     with_indices = apply_indices(splits, grid)
     assert (with_indices["index_i"] == -1).all()
     assert (with_indices["index_j"] == -1).all()
 
+    # default behaviour would split even outside grid bounds
+    splits = split_linestrings(outside, grid, bounded=False)
+    assert len(splits) > 1
+
 
 def test_split_linestrings_partial_overlap(grid):
     line = gpd.GeoDataFrame(
         {"id": [1]}, geometry=[LineString([(-2.0, 0.5), (1.5, 0.5)])]
     )
-    splits = split_linestrings(line, grid)
+    splits = split_linestrings(line, grid, bounded=True)
     coords = [list(geom.coords) for geom in splits.geometry]
     expected_coords = [
         [(-2.0, 0.5), (0.0, 0.5)],
@@ -307,6 +312,17 @@ def test_split_linestrings_partial_overlap(grid):
         )
     ) == [(-1, -1), (0, 0), (1, 0)]
 
+    # default behaviour would split even outside grid bounds
+    splits = split_linestrings(line, grid, bounded=False)
+    coords = [list(geom.coords) for geom in splits.geometry]
+    expected_coords = [
+        [(-2.0, 0.5), (-1.0, 0.5)],
+        [(-1.0, 0.5), (0.0, 0.5)],
+        [(0.0, 0.5), (1.0, 0.5)],
+        [(1.0, 0.5), (1.5, 0.5)],
+    ]
+    assert coords == expected_coords
+
 
 def test_box_geom_bounds():
     """Values take from tests/integration/range.tif"""