Merge pull request #61 from paleolimbot/geos-wk-read-crash

Fix segfault on empty project, improve plotting

DESCRIPTION

@@ -27,8 +27,7 @@ RoxygenNote: 7.1.2
 Suggests: 
     testthat (>= 2.1.0),
     vctrs,
-    sf,
-    wkutils
+    sf
 Imports: 
     libgeos (>= 3.8.1-4),
     wk (>= 0.4.1)

NEWS.md

@@ -1,5 +1,18 @@
 # geos (development version)
 
+* Add support for new features in GEOS 3.10, including
+  `geos_is_within_distance()`, `geos_prepared_is_within_distance()`,
+  `geos_read_geojson()`, `geos_write_geojson()`,
+  `geos_constrained_delaunay_triangles()`, `geos_densify()`,
+  and a new argument `flavor` to `geos_write_wkb()` and
+  `geos_write_hex()` to allow ISO WKB output instead of
+  EWKB output for geometries with Z values (#57).
+* Fix crash when EMPTY points were passed to
+  `geos_project()` or `geos_project_normalized()` (#52, #61).
+* Use `geos_simplify()`, `geos_strtree_query()`, and
+  `geos_clip_by_rect()` to speed up plotting for polygons/
+  lines with many vertices (#59, #61).
+
 # geos 0.1.2
 
 * Update tests to pass on GEOS 3.9.1 and 3.10.0 (#54).

R/geos-plot.R

@@ -2,29 +2,99 @@
 #' Plot GEOS geometries
 #'
 #' @param x A [GEOS geometry vector][as_geos_geometry]
-#' @inheritParams wkutils::wkb_plot
+#' @inheritParams wk::wk_plot
+#' @param simplify A relative tolerance to use for simplification of
+#'   geometries. Use 0 to disable simplification; use a higher number
+#'   to make simplification coarser.
+#' @param crop Use `TRUE` to crop the input to the extent of the plot.
 #'
 #' @return The input, invisibly
 #' @export
 #'
 #' @examples
-#' if (requireNamespace("wkutils")) {
-#'   plot(as_geos_geometry("LINESTRING (0 0, 1 1)"))
-#'   plot(as_geos_geometry("POINT (0.5 0.4)"), add = TRUE)
-#' }
+#' plot(as_geos_geometry("LINESTRING (0 0, 1 1)"))
+#' plot(as_geos_geometry("POINT (0.5 0.4)"), add = TRUE)
 #'
 plot.geos_geometry <- function(x, ..., asp = 1, bbox = NULL, xlab = "", ylab = "",
-                               rule = "evenodd", add = FALSE) {
-  wkutils::wkb_plot(
-    wk::as_wkb(x, include_srid = FALSE, include_z = FALSE),
-    ...,
-    asp = asp,
-    bbox = bbox,
-    xlab = xlab,
-    ylab = ylab,
-    rule = rule,
-    add = add
-  )
+                               rule = "evenodd", add = FALSE,
+                               simplify = 1, crop = TRUE) {
+  # this is too hard without vctrs (already in Suggests)
+  if (!requireNamespace("vctrs", quietly = TRUE)) {
+    stop("Package 'vctrs' is required for plot.geos_geometry()", call. = FALSE) # nocov
+  }
+
+  if (!add) {
+    plot_bbox <- if (is.null(bbox)) wk::wk_bbox(x) else bbox
+    wk::wk_plot(wk::wkb(), asp = asp, bbox = plot_bbox, xlab = xlab)
+  }
+
+  if (length(x) == 0) {
+    return(x)
+  }
+
+  # estimate resolution for simplifying
+  usr <- graphics::par("usr")
+  usr_x <- usr[1:2]
+  usr_y <- usr[3:4]
+  device_x <- graphics::grconvertX(usr_x, to = "device")
+  device_y <- graphics::grconvertY(usr_y, to = "device")
+
+  # Use resolution of 0.05 at the device level, scale to usr coords.
+  # This rarely results in simplification that will be noticed by the user.
+  scale_x <- diff(device_x) / diff(usr_x)
+  scale_y <- diff(device_y) / diff(usr_y)
+  scale <- min(abs(scale_x), abs(scale_y))
+  resolution_usr <- 0.05 / scale
+
+  x_plot <- x
+
+  # we do some subsetting of x_plot so we need to keep the dots aligned
+  dots <- list(...)
+  dots_is_vector <- vapply(dots, vctrs::vec_is, logical(1)) &
+    vapply(dots, function(x) !identical(length(x), 1L), logical(1))
+  dots_scalar <- dots[!dots_is_vector]
+  dots_vector <- dots[dots_is_vector]
+
+  # if adding, we only need features that touch the plot bbox
+  plot_area <- wk::rct(usr_x[1], usr_y[1], usr_x[2], usr_y[2], crs = wk::wk_crs(x_plot))
+  ignore_bbox <- is.null(bbox) && !add
+
+  if (!ignore_bbox) {
+    x_tree <- geos_strtree(x_plot)
+    x_touching <- geos_strtree_query(x_tree, plot_area)[[1]]
+    x_plot <- x_plot[x_touching]
+    dots_vector <- lapply(dots_vector, vctrs::vec_slice, x_touching)
+  }
+
+  if (simplify > 0) {
+    x_plot <- geos_simplify(x_plot, resolution_usr * simplify)
+  }
+
+  if (!ignore_bbox && crop) {
+    # give the crop bbox 5% so that new borders get generated outside
+    x_mid <- mean(usr_x)
+    y_mid <- mean(usr_y)
+    usr_x_expanded <- x_mid + ((usr_x - x_mid) * 1.05)
+    usr_y_expanded <- y_mid + ((usr_y - y_mid) * 1.05)
+
+    crop_area <- wk::rct(
+      usr_x_expanded[1], usr_y_expanded[1],
+      usr_x_expanded[2], usr_y_expanded[2],
+      crs = wk::wk_crs(x_plot)
+    )
+
+    x_plot <- geos_clip_by_rect(x_plot, crop_area)
+  }
+
+  # skip empties generated by simplification or clip by rect
+  x_is_empty <- geos_is_empty(x_plot)
+  x_plot <- x_plot[!x_is_empty]
+  dots_vector <- lapply(dots_vector, vctrs::vec_slice, !x_is_empty)
+
+  # pass plotting on to wk_plot()
+  if (length(x_plot) > 0) {
+    do.call(wk::wk_plot, c(list(x_plot), dots_scalar, dots_vector, list(add = TRUE)))
+  }
 
   invisible(x)
 }

src/geos-binary-atomic.c

@@ -150,7 +150,8 @@ SEXP geos_c_distance_frechet_densify(SEXP geom1, SEXP geom2, SEXP densifyFrac) {
 
 // project and project _normalized both return their result rather than use a
 // pointer arg (-1 for error, but this is undocumented)
-#define GEOS_BINARY_REAL_RETURN(_func)                                        \
+// empty points cause a segfault so we have to check for them (return NaN)
+#define GEOS_PROJECT_BINARY_REAL_RETURN(_func)                                \
   R_xlen_t size = Rf_xlength(geom1);                                          \
   SEXP result = PROTECT(Rf_allocVector(REALSXP, size));                       \
   double* pResult = REAL(result);                                             \
@@ -176,23 +177,28 @@ SEXP geos_c_distance_frechet_densify(SEXP geom1, SEXP geom2, SEXP densifyFrac) {
     geometry2 = (GEOSGeometry*) R_ExternalPtrAddr(item2);                     \
     GEOS_CHECK_GEOMETRY(geometry2, i);                                        \
                                                                               \
+    if (GEOSisEmpty_r(handle, geometry1) || GEOSisEmpty_r(handle, geometry2)) { \
+      pResult[i] = R_NaN;                                                     \
+      continue;                                                               \
+    }                                                                         \
+                                                                              \
     itemResult = _func(handle, geometry1, geometry2);                         \
     if (itemResult == -1) {                                                   \
-      Rf_error("[%d] %s", i + 1, globalErrorMessage);                                           \
+      Rf_error("[%d] %s", i + 1, globalErrorMessage);                         \
     }                                                                         \
                                                                               \
     pResult[i] = itemResult;                                                  \
   }                                                                           \
                                                                               \
-    UNPROTECT(1);                                                               \
+  UNPROTECT(1);                                                               \
   return result;
 
 SEXP geos_c_project(SEXP geom1, SEXP geom2) {
-  GEOS_BINARY_REAL_RETURN(GEOSProject_r);
+  GEOS_PROJECT_BINARY_REAL_RETURN(GEOSProject_r);
 }
 
 SEXP geos_c_project_normalized(SEXP geom1, SEXP geom2) {
-  GEOS_BINARY_REAL_RETURN(GEOSProjectNormalized_r);
+  GEOS_PROJECT_BINARY_REAL_RETURN(GEOSProjectNormalized_r);
 }
 
 

tests/testthat/test-geos-binary-atomic.R

@@ -66,13 +66,13 @@ test_that("linear referencing works", {
     c(1, NA)
   )
   expect_identical(
-    geos_project("LINESTRING (0 0, 0 10)", c("POINT (0 1)", NA)),
-    c(1, NA)
+    geos_project("LINESTRING (0 0, 0 10)", c("POINT (0 1)", "POINT EMPTY", NA)),
+    c(1, NaN, NA)
   )
 
   expect_identical(
-    geos_project_normalized("LINESTRING (0 0, 0 10)", c("POINT (0 1)", NA)),
-    c(0.1, NA)
+    geos_project_normalized("LINESTRING (0 0, 0 10)", c("POINT (0 1)", "POINT EMPTY", NA)),
+    c(0.1, NaN, NA)
   )
 })
 

tests/testthat/test-geos-plot.R

@@ -1,6 +1,18 @@
 
 test_that("ploting works", {
-  skip_if_not_installed("wkutils")
+  expect_identical(
+    plot(geos_geometry(), bbox = wk::rct(0, 0, 1, 1)),
+    geos_geometry()
+  )
+
   geom <- as_geos_geometry("LINESTRING (0 0, 1 1)")
   expect_identical(plot(geom), geom)
+
+  geom2 <- as_geos_geometry(wk::rct(0, 0, 0.4, 0.4))
+  expect_identical(plot(geom2, border = "red", col = "grey80", add = TRUE), geom2)
+
+  # make sure clip/select operations don't separate vectorized par
+  # (expecting green and brown points)
+  geom3 <- as_geos_geometry(wk::xy(c(-10, 0.5, 10, 1), c(-10, 0.5, 10, 1)))
+  plot(geom3, pch = 16, col = c("blue", "green", "purple", "brown"), add = TRUE)
 })