First stab at intersections extension module

[tlestang]

This wraps the findIntersectionsLinestring function into a C++ function that takes a Shapely LineString and raster metadata and compute splits according to grid intersections. The function returns a list of Shapely linestrings.
Its prototype is

std::vector<py::object> split(py::object linestring_py, int nrows, int ncols,
			    std::vector<double> transform);

The input Shapely object is explicitly converted to a C++ vector of points, i.e. std::vector<Vec2>

linestr convert_py2cpp(py::object linestring_py) {
  py::object coords = linestring_py.attr("coords");
  linestr linestring;
  for (int i = 0; i < py::len(coords); i++) {
    py::tuple xy = (py::tuple)coords[py::cast(i)];
    geo::Vec2<double> p((py::float_)xy[0], (py::float_)xy[1]);
    linestring.push_back(p);
  }
  return linestring;
}

The split calculation remains unchanged, based on the findIntersectionsLineString function. Finally the vector of split is converted back to a vector of Shapely linestrings:

std::vector<py::object> convert_cpp2py(std::vector<linestr> splits) {
  std::vector<py::object> splits_py;
  std::vector<std::vector<double>> split_py;
  for (auto split : splits) {
    /// Rebuild split as a vector that can be fed to
    /// the Shapely linestring constructor
    /// ...
    splits_py.push_back(SHPLY_LINESTR(split_py));
    split_py.clear();
  }
  return splits_py;
}

In the above, SHPLY_LINESTR is the Shapely LineString, imported in C++ via pybind11:

py::object SHPLY_LINESTR =
    py::module_::import("shapely.geometry").attr("LineString");

The C++ function is based on explicit conversion between Python and C++ types (see snail#14). I don't think this is the most efficient nor elegant approach, but this takes us a step further nonetheless. My may concern is performance: looping over a large geodataframe and converting each geometry might be slow. We could measure the time required to convert back and forth and compare it to the type needed to cmpute the splits.

Building and Using the module

The module can be built with

cmake -Bbuild .
cmake --build build/ --target intersections

For now, the module ends up in build/ and therefore this module must be added to sys.path. A next task is to update the setup script to build and install properly.

The module can be used as follows:

>>> import sys
>>> sys.path.append("build")
>>> from intersections import split
>>> from shapely.geometry import LineString
>>> 
>>> linestring = LineString([(0.5, 0.5), (0.75, 0.5), (1.5, 1.5)])
>>> ncols = 2
>>> nrows = 2
>>> transform = [1, 0, 0, 0, 1, 0]
>>> 
>>> splits = split(linestring, nrows, ncols, transform)

See tests/test_intersections.cpp.