begin of use with geos strtree

[mdsumner]

coords <- function(x) {
  cd <- wk::wk_coords(x)[,c("x", "y")]
  wk::xy(cd[,1L], cd[,2L], crs = wk::wk_crs(x))
}

read_geom <- function(x) {
  rga <- x$returnGeomAs
  x$returnGeomAs <- "WKB"
  out <- x$getNextFeature()
  if (is.null(out)) return(NULL)
  x$returnGeomAs <- rga
  dplyr::bind_cols(tibble::as_tibble(head(out, -1)), tibble::tibble(geometry = wk::wkb(out$geometry, crs = x$getSpatialRef())))
}

dsn <- sds::CGAZ()
v <- new(GDALVector, dsn)
#v$returnGeomAs <- "WKB"
v$resetReading()
read_geom(v)

## first time 
pts <- coords(g)
idx <- geos_strtree_query(geos_strtree(pts), pts)
idx2 <- idx[lengths(idx) > 1]
table <- tibble::tibble(id = seq_along(pts), id2 = id)
for (i in seq_along(idx2)) {
  table$id2[idx2[[i]][2]] <- idx2[[i]][1]
}
pts <- pts[sort(unique(table$id2))]

a <- 0
## now repeat 
while(!is.null(g)) {
  g <- read_geom(v)

  pts <- c(pts, coords(g))
  idx <- geos_strtree_query(geos_strtree(pts), pts)
  idx2 <- idx[lengths(idx) > 1]
  table <- dplyr::bind_rows(tibble::tibble(id = seq_along(pts) + nrow(table), id2 = id))
  for (i in seq_along(idx2)) {
    table$id2[idx2[[i]][2]] <- idx2[[i]][1]
  }
  pts <- pts[sort(unique(table$id2))]
  print(g)  
}

[mdsumner]

a bit better

coords <- function(x) {
  cd <- wk::wk_coords(x)[,c("x", "y")]
  wk::xy(cd[,1L], cd[,2L], crs = wk::wk_crs(x))
}

read_geom <- function(x, crs = NULL) {
  rga <- x$returnGeomAs
  if (is.null(crs)) crs <- x$getSpatialRef()
  x$returnGeomAs <- "WKB"
  geom <- x$getGeometryColumn()
  if (is.null(geom) || !nzchar(geom)) geom <- "geometry"
  out <- try(x$getNextFeature())
  if (is.null(out)) return(NULL)
  x$returnGeomAs <- rga
  dplyr::bind_cols(tibble::as_tibble(head(out, -1)), tibble::tibble(geometry = wk::wkb(out[[geom]], crs = crs)))
}

library(geos)
library(gdalraster)
dsn <- sds::CGAZ()
dsn <- system.file("gpkg/nc.gpkg", package = "sf", mustWork = TRUE)
v <- new(GDALVector, dsn)
v$resetReading()

table <- NULL
pts0 <- NULL
g <- read_geom(v)
## now repeat 
while(!is.null(g)) {
  if (is.null(g)) break;
  pts0 <- if (is.null(pts0))  coords(g) else  c(pts0, coords(g))
  idx <- geos_strtree_query(geos_strtree(pts0), pts0)
  idx2 <- idx[lengths(idx) > 1]
  table <- dplyr::bind_rows(table, tibble::tibble(id = seq_along(pts0) + nrow(table), id2 = id))
  for (i in seq_along(idx2)) {
    table$id2[idx2[[i]][2]] <- idx2[[i]][1]
  }
  pts0 <- pts0[sort(unique(table$id2))]

  g <- read_geom(v)

}