Integrate GeosNG into sf

NEWS.md

@@ -1,5 +1,7 @@
 # version 1.0-13
 
+* use GEOS' overlayNG routines for (GEOS) Intersection, Difference, Union and SymDifference; #2143
+
 * added `duplicated.sf()`; #2138, #2140, thanks to @bart1
 
 * `select.sf()` allows selecting the same column twice under different names; #1886

R/sample.R

@@ -35,6 +35,8 @@ st_sample = function(x, size, ...) UseMethod("st_sample")
 #' 
 #' Fibonacci sampling see: Alvaro Gonzalez, 2010. Measurement of Areas on a Sphere Using Fibonacci and Latitude-Longitude Lattices. 
 #' Mathematical Geosciences 42(1), p. 49-64
+#' 
+#' For regular sampling on the sphere, see also \code{geosphere::regularCoordinates}.
 #'
 #' Sampling methods from package \code{spatstat} are interfaced (see examples), and need their own parameters to be set. 
 #' For instance, to use \code{spatstat.random::rThomas()}, set \code{type = "Thomas"}.

configure

@@ -3660,8 +3660,8 @@ printf "%s\n" "$as_me: Install failure: compilation and/or linkage problems." >&
    as_fn_error $? "cannot link projection code" "$LINENO" 5
 fi
 
-{ printf "%s\n" "$as_me:${as_lineno-$LINENO}: checking GDAL: checking whether PROJ is available fur running:" >&5
-printf %s "checking GDAL: checking whether PROJ is available fur running:... " >&6; }
+{ printf "%s\n" "$as_me:${as_lineno-$LINENO}: checking GDAL: checking whether PROJ is available for running:" >&5
+printf %s "checking GDAL: checking whether PROJ is available for running:... " >&6; }
 ./gdal_proj
 if test `echo $?` -ne 0 ; then
 gdal_has_proj=no

configure.ac

@@ -288,7 +288,7 @@ if test "${gdal_has_proj}" = no; then
    AC_MSG_ERROR([cannot link projection code])
 fi
 
-AC_MSG_CHECKING(GDAL: checking whether PROJ is available fur running:)
+AC_MSG_CHECKING(GDAL: checking whether PROJ is available for running:)
 ./gdal_proj
 if test `echo $?` -ne 0 ; then
 gdal_has_proj=no

src/geos.cpp

@@ -62,6 +62,7 @@ typedef char (* log_fn)(GEOSContextHandle_t, const GEOSGeometry *, const GEOSGeo
 typedef char (* log_prfn)(GEOSContextHandle_t, const GEOSPreparedGeometry *,
 	const GEOSGeometry *);
 typedef GEOSGeom (* geom_fn)(GEOSContextHandle_t, const GEOSGeom, const GEOSGeom);
+typedef GEOSGeom (* geom_fnp)(GEOSContextHandle_t, const GEOSGeom, const GEOSGeom, double grid_size);
 
 static int notice = 0; // global var to silently catch notice of illegal geoms, e.g. non-closed rings
 
@@ -170,6 +171,17 @@ static std::vector<GEOSGeometry*> to_raw(std::vector<GeomPtr> & g) {
 	return raw;
 }
 
+double geos_grid_size(Rcpp::List x) {
+	double precision = x.attr("precision");
+	if (precision != 0.0)
+		precision = 1. / precision;
+	return precision;
+}
+
+double geos_grid_size_xy(Rcpp::List x, Rcpp::List y) {
+	return std::max(geos_grid_size(x), geos_grid_size(y));
+}
+
 std::vector<GeomPtr> geometries_from_sfc(GEOSContextHandle_t hGEOSCtxt, Rcpp::List sfc, int *dim = NULL, bool stop_on_NULL = true) {
 
 	Rcpp::List sfc_cls = get_dim_sfc(sfc);
@@ -185,10 +197,8 @@ std::vector<GeomPtr> geometries_from_sfc(GEOSContextHandle_t hGEOSCtxt, Rcpp::Li
 		Rcpp::stop("GEOS does not support XYM or XYZM geometries; use st_zm() to drop M\n"); // #nocov
 
 #ifdef HAVE_390
-	double precision = sfc.attr("precision");
-	bool set_precision = precision != 0.0;
-	if (set_precision)
-		precision = 1/precision;
+	double grid_size = geos_grid_size(sfc);
+	bool set_precision = grid_size != 0.0;
 	sfc.attr("precision") = 0.0; // so that CPL_write_wkb doesn't do the rounding;
 #endif
 	Rcpp::List wkblst = CPL_write_wkb(sfc, true);
@@ -205,7 +215,7 @@ std::vector<GeomPtr> geometries_from_sfc(GEOSContextHandle_t hGEOSCtxt, Rcpp::Li
 		}
 #ifdef HAVE_390
 		else if (set_precision)
-			g[i] = geos_ptr(GEOSGeom_setPrecision_r(hGEOSCtxt, g[i].get(), precision, GEOS_PREC_VALID_OUTPUT), hGEOSCtxt);
+			g[i] = geos_ptr(GEOSGeom_setPrecision_r(hGEOSCtxt, g[i].get(), grid_size, GEOS_PREC_VALID_OUTPUT), hGEOSCtxt);
 #endif
 	}
 	GEOSWKBReader_destroy_r(hGEOSCtxt, wkb_reader);
@@ -933,6 +943,7 @@ Rcpp::List CPL_geos_op2(std::string op, Rcpp::List sfcx, Rcpp::List sfcy) {
 	std::vector<GeomPtr> out;
 	std::vector<double> index_x, index_y;
 	std::vector<size_t> items(x.size());
+	double grid_size = geos_grid_size_xy(sfcx, sfcy);
 
 	if (op == "intersection") {
 
@@ -955,7 +966,11 @@ Rcpp::List CPL_geos_op2(std::string op, Rcpp::List sfcx, Rcpp::List sfcy) {
 			std::sort(sel.begin(), sel.end());
 			for (size_t item = 0; item < sel.size(); item++) {
 				size_t j = sel[item];
+#ifndef HAVE_390
 				GeomPtr geom = geos_ptr(GEOSIntersection_r(hGEOSCtxt, x[j].get(), y[i].get()), hGEOSCtxt);
+#else
+				GeomPtr geom = geos_ptr(GEOSIntersectionPrec_r(hGEOSCtxt, x[j].get(), y[i].get(), grid_size), hGEOSCtxt);
+#endif
 				if (geom == nullptr)
 					Rcpp::stop("GEOS exception"); // #nocov
 				if (! chk_(GEOSisEmpty_r(hGEOSCtxt, geom.get()))) {
@@ -968,19 +983,33 @@ Rcpp::List CPL_geos_op2(std::string op, Rcpp::List sfcx, Rcpp::List sfcy) {
 		}
 
 	} else {
+#ifndef HAVE_390
 		geom_fn geom_function;
 		if (op == "union")
 			geom_function = (geom_fn) GEOSUnion_r;
 		else if (op == "difference")
 			geom_function = (geom_fn) GEOSDifference_r;
 		else if (op == "sym_difference")
 			geom_function = (geom_fn) GEOSSymDifference_r;
+#else
+		geom_fnp geom_function;
+		if (op == "union")
+			geom_function = (geom_fn) GEOSUnionPrec_r;
+		else if (op == "difference")
+			geom_function = (geom_fn) GEOSDifferencePrec_r;
+		else if (op == "sym_difference")
+			geom_function = (geom_fn) GEOSSymDifferencePrec_r;
+#endif
 		else
 			Rcpp::stop("invalid operation"); // #nocov
 
 		for (size_t i = 0; i < y.size(); i++) {
 			for (size_t j = 0; j < x.size(); j++) {
+#ifndef HAVE_390
 				GeomPtr geom = geos_ptr(geom_function(hGEOSCtxt, x[j].get(), y[i].get()), hGEOSCtxt);
+#else
+				GeomPtr geom = geos_ptr(geom_function(hGEOSCtxt, x[j].get(), y[i].get(), grid_size), hGEOSCtxt);
+#endif
 				if (geom == nullptr)
 					Rcpp::stop("GEOS exception"); // #nocov
 				if (! chk_(GEOSisEmpty_r(hGEOSCtxt, geom.get()))) {
@@ -1156,6 +1185,7 @@ Rcpp::List CPL_nary_difference(Rcpp::List sfc) {
 	GEOSContextHandle_t hGEOSCtxt = CPL_geos_init();
 	std::vector<GeomPtr> x = geometries_from_sfc(hGEOSCtxt, sfc, &dim);
 	std::vector<GeomPtr> out;
+	double grid_size = geos_grid_size(sfc);
 	// initialize trees to find overlapping areas quickly
 	for (size_t i = 0; i < x.size(); i++) {
 		// if i'th geometry in x is empty then skip it
@@ -1185,7 +1215,11 @@ Rcpp::List CPL_nary_difference(Rcpp::List sfc) {
 					// test if the items intersect with geom
 					if (chk_(GEOSIntersects_r(hGEOSCtxt, geom.get(), out[tree_sel[j]].get()))) {
 						// if they do then erase overlapping parts from geom
+#ifndef HAVE_390
 						geom = geos_ptr(GEOSDifference_r(hGEOSCtxt, geom.get(), out[tree_sel[j]].get()), hGEOSCtxt);
+#else
+						geom = geos_ptr(GEOSDifferencePrec_r(hGEOSCtxt, geom.get(), out[tree_sel[j]].get(), grid_size), hGEOSCtxt);
+#endif
 						if (geom == nullptr)
 							Rcpp::stop("GEOS exception"); // #nocov
 						// ensure that geom is valid
@@ -1221,6 +1255,7 @@ Rcpp::List CPL_nary_intersection(Rcpp::List sfc) {
 	std::vector<GeomPtr> x = geometries_from_sfc(hGEOSCtxt, sfc, &dim);
 	std::vector<GeomPtr> out;
 	int errors = 0;
+	double grid_size = geos_grid_size(sfc);
 #ifdef HAVE350
 	notice = 0;
 	GEOSContext_setNoticeMessageHandler_r(hGEOSCtxt,
@@ -1249,7 +1284,11 @@ Rcpp::List CPL_nary_intersection(Rcpp::List sfc) {
 				// iterate over items in query and erase overlapping areas in geom
 				for (size_t j = 0; j < tree_sel.size(); j++) {
 					size_t k = tree_sel[j];
+#ifndef HAVE_390
 					GeomPtr inters = geos_ptr(GEOSIntersection_r(hGEOSCtxt, out[k].get(), geom.get()), hGEOSCtxt);
+#else
+					GeomPtr inters = geos_ptr(GEOSIntersectionPrec_r(hGEOSCtxt, out[k].get(), geom.get(), grid_size), hGEOSCtxt);
+#endif
 					if (geom.get() != nullptr) {
 						if (inters == nullptr)
 							errors++;
@@ -1259,7 +1298,11 @@ Rcpp::List CPL_nary_intersection(Rcpp::List sfc) {
 							if (geom == nullptr)
 								Rcpp::stop("GEOS exception"); // #nocov
 							// cut out inters from out[k]:
+#ifndef HAVE_390
 							GeomPtr g = geos_ptr(GEOSDifference_r(hGEOSCtxt, out[k].get(), inters.get()), hGEOSCtxt); 
+#else
+							GeomPtr g = geos_ptr(GEOSDifferencePrec_r(hGEOSCtxt, out[k].get(), inters.get(), grid_size), hGEOSCtxt); 
+#endif
 							if (g == nullptr)
 								Rcpp::warning("GEOS difference returns NULL"); // #nocov
 							else {

src/stars.cpp

@@ -366,6 +366,7 @@ List CPL_read_gdal(CharacterVector fname, CharacterVector options, CharacterVect
 	CharacterVector descriptions(bands.size());
 	NumericMatrix ranges(bands.size(), 4);
 	IntegerMatrix blocksizes(bands.size(), 2);
+	IntegerVector colorInterp(bands.size());
 	for (int i = 0; i < bands.size(); i++) {
 		if ((poBand = poDataset->GetRasterBand(bands(i))) == NULL)
 			stop("trying to read a band that is not present");
@@ -391,6 +392,7 @@ List CPL_read_gdal(CharacterVector fname, CharacterVector options, CharacterVect
 		poBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
 		blocksizes(i, 0) = nBlockXSize;
 		blocksizes(i, 1) = nBlockYSize;
+		colorInterp(i) = (int) poBand->GetColorInterpretation();
 	}
 
 	// get metadata items:
@@ -462,7 +464,8 @@ List CPL_read_gdal(CharacterVector fname, CharacterVector options, CharacterVect
 		_["blocksizes"] = blocksizes,
 		_["descriptions"] = descriptions,
 		_["default_geotransform"] = default_geotransform,
-		_["proxy"] = LogicalVector::create(!read_data)
+		_["proxy"] = LogicalVector::create(!read_data),
+		_["colorInterp"] = colorInterp
 	);
 	if (read_data) {
 		ReturnList.attr("data") = read_gdal_data(poDataset, nodatavalue, nXOff, nYOff,