paleolimbot/geos

provide matrix version of relate pattern match

JosiahParry opened this issue · 1 comments

geos_relate_pattern() is a really useful function for comparing two individual geometries. It would be great to have this same functionality made available in the _matrix() format. If made available, it could be used to speed up neighbor matching in spdep.

Looking at queen contiguity with self included using geos_intersects_matrix() is very fast. Rook contiguity requires a relate which the existing functions do not provide a way to harness the tree.

library(spdep)

library(geos)

x <- sf::st_as_sf(Guerry::gfrance85) |> 
  sf::st_geometry()

bench::mark(
  spdep = include.self(poly2nb(x)),
  geos = {
    tree <- geos_strtree(x)
    geos_intersects_matrix(x, tree) 
  },
  check = FALSE
)
#> # A tibble: 2 × 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 spdep        39.6ms   40.2ms      24.9      20MB     24.9
#> 2 geos         14.5ms   15.1ms      65.8     344KB      0

It's certainly possible to provide the same interface; however, a generic version that uses the tree might require some explicit inspection of the relate pattern which I don't really know much about. The geos_strtree_query() function will get you potential intersection slightly faster than a matrix intersection and you could do lapply() + geos_relate_pattern(). Also, there is geos_basic_strtree() which would let you do something like geos_relate_pattern(vec_a[query_result[[1]]], vec_b[query_result[[2]]) (i.e., skipping the lapply).

Probably also possible to wire up a matrix version but it would probably involve skipping patterns where the "potential intersection" filter would end up generating incorrect results (e.g., disjoint).