New transformer: s2_convex_hull

NAMESPACE

@@ -103,6 +103,7 @@ export(s2_closest_feature)
 export(s2_closest_point)
 export(s2_contains)
 export(s2_contains_matrix)
+export(s2_convex_hull_agg)
 export(s2_coverage_union_agg)
 export(s2_covered_by)
 export(s2_covered_by_matrix)

R/RcppExports.R

@@ -413,6 +413,10 @@ cpp_s2_buffer_cells <- function(geog, distance, maxCells, minLevel) {
     .Call(`_s2_cpp_s2_buffer_cells`, geog, distance, maxCells, minLevel)
 }
 
+cpp_s2_convex_hull_agg <- function(geog, s2options) {
+    .Call(`_s2_cpp_s2_convex_hull_agg`, geog, s2options)
+}
+
 s2_xptr_test <- function(size) {
     .Call(`_s2_s2_xptr_test`, size)
 }

R/s2-transformers.R

@@ -94,6 +94,16 @@
 #'   s2_options(snap = s2_snap_level(30))
 #' )
 #'
+#' # s2_convex_hull_agg builds the convex hull of a list of geometries
+#' s2_convex_hull_agg(
+#'   c(
+#'     "POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))",
+#'     "POLYGON ((5 5, 15 5, 15 15, 5 15, 5 5))"
+#'   ),
+#'   # 32 bit platforms may need to set snap rounding
+#'   s2_options(snap = s2_snap_level(30))
+#' )
+#'
 #' # use s2_union_agg() to aggregate geographies in a vector
 #' s2_coverage_union_agg(
 #'   c(
@@ -221,6 +231,11 @@ s2_union_agg <- function(x, options = s2_options(), na.rm = FALSE) {
   new_s2_xptr(cpp_s2_union_agg(s2_union(x, options = options), options, na.rm), "s2_geography")
 }
 
+#' @rdname s2_boundary
+#' @export
+s2_convex_hull_agg <- function(x, options = s2_options() ) {
+  new_s2_xptr(cpp_s2_convex_hull_agg(as_s2_geography(x), options), "s2_geography")
+}
 
 #' Linear referencing
 #'
@@ -267,3 +282,4 @@ s2_interpolate_normalized <- function(x, distance_normalized) {
     "s2_geography"
   )
 }
+

src/RcppExports.cpp

@@ -1252,6 +1252,18 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// cpp_s2_convex_hull_agg
+List cpp_s2_convex_hull_agg(List geog, List s2options);
+RcppExport SEXP _s2_cpp_s2_convex_hull_agg(SEXP geogSEXP, SEXP s2optionsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< List >::type geog(geogSEXP);
+    Rcpp::traits::input_parameter< List >::type s2options(s2optionsSEXP);
+    rcpp_result_gen = Rcpp::wrap(cpp_s2_convex_hull_agg(geog, s2options));
+    return rcpp_result_gen;
+END_RCPP
+}
 // s2_xptr_test
 List s2_xptr_test(R_xlen_t size);
 RcppExport SEXP _s2_s2_xptr_test(SEXP sizeSEXP) {
@@ -1382,6 +1394,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_s2_cpp_s2_unary_union", (DL_FUNC) &_s2_cpp_s2_unary_union, 2},
     {"_s2_cpp_s2_interpolate_normalized", (DL_FUNC) &_s2_cpp_s2_interpolate_normalized, 2},
     {"_s2_cpp_s2_buffer_cells", (DL_FUNC) &_s2_cpp_s2_buffer_cells, 4},
+    {"_s2_cpp_s2_convex_hull_agg", (DL_FUNC) &_s2_cpp_s2_convex_hull_agg, 2},
     {"_s2_s2_xptr_test", (DL_FUNC) &_s2_s2_xptr_test, 1},
     {"_s2_s2_xptr_test_op", (DL_FUNC) &_s2_s2_xptr_test_op, 1},
     {"c_s2_coord_filter_new",        (DL_FUNC) &c_s2_coord_filter_new,        4},

src/s2-transformers.cpp

@@ -14,6 +14,7 @@
 #include "s2/s2builderutil_snap_functions.h"
 #include "s2/s2shape_index_buffered_region.h"
 #include "s2/s2region_coverer.h"
+#include "s2/s2convex_hull_query.h"
 
 #include "s2-options.h"
 #include "geography-operator.h"
@@ -671,3 +672,40 @@ List cpp_s2_buffer_cells(List geog, NumericVector distance, int maxCells, int mi
   Op op(distance, maxCells, minLevel);
   return op.processVector(geog);
 }
+
+
+
+// [[Rcpp::export]]
+List cpp_s2_convex_hull_agg(List geog, List s2options) {
+  // create the convex hull query
+  S2ConvexHullQuery convexHullQuery;
+  SEXP item;
+  for (R_xlen_t i = 0; i < geog.size(); i++) {
+    item = geog[i];
+    if (item != R_NilValue) {
+      Rcpp::XPtr<Geography> feature(item);
+      if (feature->GeographyType() == Geography::Type::GEOGRAPHY_POLYGON) {
+          const S2Polygon* pol = feature->Polygon();
+          convexHullQuery.AddPolygon(*pol);
+      } else if (feature->GeographyType() == Geography::Type::GEOGRAPHY_POINT) {
+          const std::vector<S2Point>* pts = feature->Point();
+          for(const auto& pt: *pts) {
+            convexHullQuery.AddPoint(pt);
+          }
+      } else if (feature->GeographyType() == Geography::Type::GEOGRAPHY_POLYLINE) {
+          const std::vector<std::unique_ptr<S2Polyline>>* lins = feature->Polyline();
+          for(const auto& lin: *lins) {
+              convexHullQuery.AddPolyline(*lin);
+          }
+      } else if (feature->GeographyType() == Geography::Type::GEOGRAPHY_COLLECTION) {
+          if (!feature->IsEmpty()) {
+              Rcpp::stop("GeometryCollection is not supported");
+          }
+      }
+    }
+  }
+  // Builds the convex hull and returns the polygon as a geography
+  std::unique_ptr<S2Polygon> outP = absl::make_unique<S2Polygon>(convexHullQuery.GetConvexHull());
+  XPtr<Geography> outG(new PolygonGeography(std::move(outP)));
+  return List::create(outG);
+}

tests/testthat/test-s2-transformers.R

@@ -620,3 +620,35 @@ test_that("s2_interpolate() and s2_interpolate_normalized() work", {
     "must be a simple geography"
   )
 })
+
+
+test_that("s2_convex_hull_agg() works", {
+  expect_equal(
+    s2_area(s2_convex_hull_agg(c("POINT(3.6 43.2)", "POINT (0 0)", "POINT(3.61 43.21)"))),
+    s2_area("POLYGON ((0 0, 3.61 43.21, 3.6 43.2, 0 0))")
+  )
+
+  expect_equal(
+    s2_area(s2_convex_hull_agg(c(
+    "POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))",
+    "POLYGON ((5 5, 15 5, 15 15, 5 15, 5 5))"))),
+    s2_area("POLYGON ((0 0, 10 0, 15 5, 15 15, 5 15, 0 10, 0 0))")
+  )
+
+  expect_equal(
+    s2_area(s2_convex_hull_agg(c(
+     "POINT (3.6 43.2)",
+     "LINESTRING (3.49 43.05, 3.52 43.1, 3.38 43.2, 3.1 43.1)",
+     "POLYGON ((3.01 43.2, 3.4 44.01, 3.5 43.5, 3.1 43.2, 3.01 43.2))",
+     "GEOMETRYCOLLECTION EMPTY"
+    ))),
+    s2_area(
+      "POLYGON ((3.49 43.05, 3.6 43.2, 3.4 44.01, 3.01 43.2, 3.1 43.1, 3.49 43.05))"
+    )
+  )
+
+  expect_error(
+    s2_convex_hull_agg(c("GEOMETRYCOLLECTION(POLYGON ((3.01 43.2, 3.4 44.01, 3.5 43.5, 3.1 43.2, 3.01 43.2)), POINT (3.6 43.2))")),
+    "GeometryCollection is not supported"
+  )
+})