Merge pull request #163 from r-spatial/convex-hull

Add `s2_convex_hull()`, add `na.rm` argument to `s2_convex_hull_agg()`

NAMESPACE

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

NEWS.md

@@ -3,6 +3,8 @@
 - Fix for s2 build on Windows with R <= 3.6.x (#142)
 - Fix for s2 build on MacOS with multiple openssl versions (#142, #145, #146)
 - Fix for s2 build on 32-bit openssl (#143, #147)
+- Added `s2_convex_hull()` and `s2_convex_hull_agg()` (@spiry34, #150, 
+  #151, #163).
 - Added `max_distance` argument to `s2_closest_edges()`, making
   distance-constrained k-nearest neighbours possible (#125, #156, #162).
 - Added a spherical `s2_point_on_surface()` implementation for polygons

R/RcppExports.R

@@ -473,8 +473,12 @@ 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)
+cpp_s2_convex_hull <- function(geog) {
+    .Call(`_s2_cpp_s2_convex_hull`, geog)
+}
+
+cpp_s2_convex_hull_agg <- function(geog, naRm) {
+    .Call(`_s2_cpp_s2_convex_hull_agg`, geog, naRm)
 }
 
 s2_xptr_test <- function(size) {

R/s2-transformers.R

@@ -52,7 +52,7 @@
 #'
 #' # s2_point_on_surface guarantees a point on surface
 #' # Note: this is not the same as st_point_on_surface
-#' s2_centroid("POLYGON ((0 0, 10 0, 1 1, 0 10, 0 0)"")
+#' s2_centroid("POLYGON ((0 0, 10 0, 1 1, 0 10, 0 0))")
 #' s2_point_on_surface("POLYGON ((0 0, 10 0, 1 1, 0 10, 0 0))")
 #'
 #' # returns the unweighted centroid of the entire input
@@ -104,9 +104,7 @@
 #'   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
@@ -213,6 +211,12 @@ s2_buffer_cells <- function(x, distance, max_cells = 1000, min_level = -1,
   new_s2_xptr(cpp_s2_buffer_cells(recycled[[1]], recycled[[2]], max_cells, min_level), "s2_geography")
 }
 
+#' @rdname s2_boundary
+#' @export
+s2_convex_hull <- function(x) {
+  new_s2_xptr(cpp_s2_convex_hull(as_s2_geography(x)), "s2_geography")
+}
+
 #' @rdname s2_boundary
 #' @export
 s2_centroid_agg <- function(x, na.rm = FALSE) {
@@ -239,8 +243,8 @@ s2_union_agg <- function(x, options = s2_options(), na.rm = FALSE) {
 
 #' @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")
+s2_convex_hull_agg <- function(x, na.rm = FALSE) {
+  new_s2_xptr(cpp_s2_convex_hull_agg(as_s2_geography(x), na.rm), "s2_geography")
 }
 
 #' Linear referencing

src/RcppExports.cpp

@@ -1436,15 +1436,26 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// cpp_s2_convex_hull
+List cpp_s2_convex_hull(List geog);
+RcppExport SEXP _s2_cpp_s2_convex_hull(SEXP geogSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< List >::type geog(geogSEXP);
+    rcpp_result_gen = Rcpp::wrap(cpp_s2_convex_hull(geog));
+    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) {
+List cpp_s2_convex_hull_agg(List geog, bool naRm);
+RcppExport SEXP _s2_cpp_s2_convex_hull_agg(SEXP geogSEXP, SEXP naRmSEXP) {
 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));
+    Rcpp::traits::input_parameter< bool >::type naRm(naRmSEXP);
+    rcpp_result_gen = Rcpp::wrap(cpp_s2_convex_hull_agg(geog, naRm));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -1593,6 +1604,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", (DL_FUNC) &_s2_cpp_s2_convex_hull, 1},
     {"_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},

src/geography-collection.h

@@ -121,6 +121,10 @@ class GeographyCollection: public Geography {
     return shapeIds;
   }
 
+  const std::vector<std::unique_ptr<Geography> >* CollectionFeatures() {
+    return &(this->features);
+  }
+
   void Export(WKGeometryHandler* handler, uint32_t partId) {
     WKGeometryMeta meta(WKGeometryType::GeometryCollection, false, false, false);
     meta.hasSize = true;

src/geography.h

@@ -78,6 +78,10 @@ class Geography {
     return nullptr;
   }
 
+  virtual const std::vector<std::unique_ptr<Geography>>* CollectionFeatures() {
+    return nullptr;
+  }
+
   // other calculations use ShapeIndex
   virtual S2ShapeIndex* ShapeIndex() {
     if (!this->hasIndex) {

src/s2-transformers.cpp

@@ -731,38 +731,71 @@ List cpp_s2_buffer_cells(List geog, NumericVector distance, int maxCells, int mi
 }
 
 
+class ConvexHullGeographyQuery: public S2ConvexHullQuery {
+public:
+
+  void AddGeography(Geography* feature) {
+    if (feature->GeographyType() == Geography::Type::GEOGRAPHY_POLYGON) {
+      const S2Polygon* pol = feature->Polygon();
+      this->AddPolygon(*pol);
+    } else if (feature->GeographyType() == Geography::Type::GEOGRAPHY_POINT) {
+      const std::vector<S2Point>* pts = feature->Point();
+      for(const auto& pt: *pts) {
+        this->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) {
+        this->AddPolyline(*lin);
+      }
+    } else if (feature->GeographyType() == Geography::Type::GEOGRAPHY_COLLECTION) {
+      const std::vector<std::unique_ptr<Geography>>* features = feature->CollectionFeatures();
+      for (const auto& feat: *features) {
+        this->AddGeography(feat.get());
+      }
+    }
+  }
+
+  std::unique_ptr<S2Polygon> GetConvexHullPolygon() {
+    return absl::make_unique<S2Polygon>(this->GetConvexHull());
+  }
+};
 
 // [[Rcpp::export]]
-List cpp_s2_convex_hull_agg(List geog, List s2options) {
+List cpp_s2_convex_hull(List geog) {
+  class Op: public UnaryGeographyOperator<List, SEXP> {
+    SEXP processFeature(XPtr<Geography> feature, R_xlen_t i) {
+      ConvexHullGeographyQuery convexHullQuery;
+      convexHullQuery.AddGeography(feature);
+
+      std::unique_ptr<S2Polygon> outP = convexHullQuery.GetConvexHullPolygon();
+      return XPtr<Geography>(new PolygonGeography(std::move(outP)));
+    }
+  };
+
+  Op op;
+  return op.processVector(geog);
+}
+
+
+// [[Rcpp::export]]
+List cpp_s2_convex_hull_agg(List geog, bool naRm) {
   // create the convex hull query
-  S2ConvexHullQuery convexHullQuery;
+  ConvexHullGeographyQuery convexHullQuery;
   SEXP item;
   for (R_xlen_t i = 0; i < geog.size(); i++) {
     item = geog[i];
+    if (item == R_NilValue && !naRm) {
+      return List::create(R_NilValue);
+    }
+
     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");
-          }
-      }
+      XPtr<Geography> feature(item);
+      convexHullQuery.AddGeography(feature);
     }
   }
-  // Builds the convex hull and returns the polygon as a geography
-  std::unique_ptr<S2Polygon> outP = absl::make_unique<S2Polygon>(convexHullQuery.GetConvexHull());
+
+  std::unique_ptr<S2Polygon> outP = convexHullQuery.GetConvexHullPolygon();
   XPtr<Geography> outG(new PolygonGeography(std::move(outP)));
   return List::create(outG);
 }