From ad34f461b4085aabfd342e1850c3d0c32f352518 Mon Sep 17 00:00:00 2001 From: Edzer Pebesma Date: Fri, 7 Apr 2023 22:03:21 +0200 Subject: [PATCH 1/3] address #2143 --- NEWS.md | 2 ++ src/geos.cpp | 53 +++++++++++++++++++++++++++++++++++++++++++++++----- 2 files changed, 50 insertions(+), 5 deletions(-) diff --git a/NEWS.md b/NEWS.md index f81a8182a..a3d439674 100644 --- a/NEWS.md +++ b/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 diff --git a/src/geos.cpp b/src/geos.cpp index 3fb24d7ed..da73006b3 100644 --- a/src/geos.cpp +++ b/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 to_raw(std::vector & 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 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 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 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 out; std::vector index_x, index_y; std::vector 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,6 +983,7 @@ 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; @@ -975,12 +991,25 @@ Rcpp::List CPL_geos_op2(std::string op, Rcpp::List sfcx, Rcpp::List sfcy) { 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 x = geometries_from_sfc(hGEOSCtxt, sfc, &dim); std::vector 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 x = geometries_from_sfc(hGEOSCtxt, sfc, &dim); std::vector 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 { From 01c549ccdd8ae31da5c42399826313e96dda85a7 Mon Sep 17 00:00:00 2001 From: Edzer Pebesma Date: Sun, 9 Apr 2023 20:59:22 +0200 Subject: [PATCH 2/3] fix typo see https://github.com/rspatial/terra/issues/1106 --- configure | 4 ++-- configure.ac | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/configure b/configure index 262b5bc58..81443f5cd 100755 --- a/configure +++ b/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 diff --git a/configure.ac b/configure.ac index 819cd5430..aabdfd1e1 100644 --- a/configure.ac +++ b/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 From 4d4f2c667ca4baa1732f880884bc296e1d52fdd7 Mon Sep 17 00:00:00 2001 From: Edzer Pebesma Date: Mon, 10 Apr 2023 21:58:13 +0200 Subject: [PATCH 3/3] add ref to geosphere::regularCoordinates --- R/sample.R | 2 ++ man/st_sample.Rd | 2 ++ src/stars.cpp | 5 ++++- 3 files changed, 8 insertions(+), 1 deletion(-) diff --git a/R/sample.R b/R/sample.R index 68592f511..7e889cc03 100644 --- a/R/sample.R +++ b/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"}. diff --git a/man/st_sample.Rd b/man/st_sample.Rd index 6424bfebd..3fc76907b 100644 --- a/man/st_sample.Rd +++ b/man/st_sample.Rd @@ -69,6 +69,8 @@ As parameter called \code{offset} can be passed to control ("fix") regular or he 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"}. } diff --git a/src/stars.cpp b/src/stars.cpp index 6b5f2b0e7..c2a5189da 100644 --- a/src/stars.cpp +++ b/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,