speed prospects with {wk}
Opened this issue · 2 comments
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
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)))
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 ...