speed prospects with {wk}

[mdsumner]

SC0_sf <- function(x) {
  coords <- wk::wk_coords(x)
  udata <- unjoin::unjoin(coords[c("x", "y", "ring_id", "feature_id")], .data$x, .data$y)
  
  vertex <- setNames(udata$.idx0[c("x", "y", ".idx0")], c("x_", "y_", "vertex_"))
  udata$data$coord <- 1:nrow(udata$data)
  segs <- udata$data %>% dplyr::select(.data$ring_id, .data$coord, .data$feature_id)  %>%
    dplyr::mutate(.cx0 = .data$coord,   ## specify in segment terms
                  .cx1 = .data$coord + 1L) %>%
    dplyr::group_by(.data$ring_id) %>% dplyr::slice(-dplyr::n()) %>% dplyr::ungroup() %>%
    dplyr::transmute(.data$.cx0, .data$.cx1, path_ = .data$ring_id, object = .data$feature_id)
  
  object <- sc_object(x)
  segs[[".vx0"]] <- udata$data$.idx0[match(segs$.cx0, udata$data$coord)]
  segs[[".vx1"]] <- udata$data$.idx0[match(segs$.cx1, udata$data$coord)]
  ## but udata$.idx0 has the vertices, with .idx0 as the mapping
  object$topology_ <- split(segs[c(".vx0", ".vx1", "path_")], segs$object)
  
  meta <- tibble::tibble(proj = crsmeta::crs_wkt(x), ctime = Sys.time())
  vertex <- vertex %>%
    dplyr::arrange(.data$vertex_)
  
  bad <- !names(vertex) %in% c("x_", "y_", "z_", "m_", "t_")
  if (any(bad)) vertex <- vertex[!bad]
  structure(list(object = object, vertex = vertex,
                 meta = meta),
            class = c("SC0", "sc"))
}

above use 'parcels_hobart' i.e.

f <- raadfiles::thelist_files(pattern = "parcels_hobart")$fullname
x <- sf::read_sf(f)
 SC0(x)
class       : SC0
type        : Structural
vertices    : 195295 (2-space)
primitives  : 463141 (2-space)
crs         : GDA94 / MGA zone 55

also see for a more exhaustive exploration: dcooley/geometries#4

[mdsumner]

here's some sf helpers, segments and triangles (to avoid closing vertex dupe) as linestrings

sf_segment <- function(x) {
  sc <- silicate::SC0(x)
  topol <- silicate::sc_object(sc)$topology_
  idx <- t(as.matrix(do.call(rbind, topol)[c(".vx0", ".vx1")]))
  d <- silicate::sc_vertex(sc)[as.vector(idx), c("x_", "y_")]
  d$feature_id <- rep(seq_len(dim(idx)[2L]), each = 2L)
  out <- sfheaders::sf_linestring(d, x = "x_", y = "y_", linestring_id = "feature_id")
  out$feature_id <- rep(seq_along(topol), unlist(lapply(topol, nrow)))
  out
}


sf_triangle <- function(x, verts = c("x_", "y_")) {
  if (inherits(x, "TRI0")) {
    ## might alread be sc DEL0 or TRI0
    sc <- x
    } else {
      sc <- silicate::TRI0(x)
  }
  topol <- silicate::sc_object(sc)$topology_
  idx <- t(as.matrix(do.call(rbind, topol)[c(".vx0", ".vx1", ".vx2")]))
  d <- silicate::sc_vertex(sc)[as.vector(idx), verts]
  d$feature_id <- rep(seq_len(dim(idx)[2L]), each = 3L)
  out <- sfheaders::sf_linestring(d, x = verts[1], y = verts[2L], linestring_id = "feature_id")
  out$feature_id <- rep(seq_along(topol), unlist(lapply(topol, nrow)))
  out
}
example(read_sf)
plot(sf_segment(nc))
plot(sf_triangle(minimal_mesh))
plot(sf_triangle(anglr::DEL0(nc, max_area = .02)))

[mdsumner]

what are the key pathways? for example, turning sf polygons or lines into segments does not require vertex deduplication, so rather than go through SC0 every time we really only need the map of paths (number of coordinates in each), to use as a grouping

I think fleshing that list out, in light of traipse, rayshader, mapdeck, and various discussions around the place is worthwhile - like, when do you want segments? sometimes you want them all and to know what polygon they belong to, other times you need an actual mesh of segments for planar partition ... this is where sf fails us, we don't want to copy out every vertex instance ...