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)
  
}

[mdsumner]

I started again (maybelook at earlier attempt)

library(silicate)
library(wk)
x <- dplyr::sample_n(inlandwaters, 1)

wk <- wk::as_wkb(x)
gib <- gibble::gibble(x)
grp <- rep(seq_along(gib$nrow), gib$nrow)

verts <- tibble::tibble(vert = wk_vertices(wk))
tree <- geos::geos_strtree(verts)
id <- geos::geos_strtree_query(tree, verts)

idx <- seq_along(id)
repeats <- unique(id[lengths(id) > 1])
for (i in seq_along(repeats)) {
  mid <- repeats[[i]][1]
  for (j in 2:length(repeats[[i]])) {
     idx[repeats[[i]][j]] <- mid
  }
}

dd <- unjoin::unjoin(tibble::tibble(id = idx), id)
segs <- do.call(rbind, 
                lapply(split(dd$data$.idx0, grp), 
                       silicate:::path_to_segment))


uniq_verts <- verts[dd$.idx0$id, ]
cds <- wk_coords(uniq_verts)
plot(cds$x, cds$y, type = "p", pch = ".")
segments(cds$x[segs$.vx0], cds$y[segs$.vx0], 
         cds$x[segs$.vx1], cds$y[segs$.vx1], col = sample(hcl.colors(nrow(segs))), lwd = 2)