hypertidy / anglr

Mesh creation and topology for spatial data (and not just geographic)
https://hypertidy.github.io/anglr/
83 stars 10 forks source link

lines, and refactor #14

Closed mdsumner closed 6 years ago

mdsumner commented 8 years ago

include ability to "copy down" object attributes to vertices, for e.g. contour lines

need to think about how maintaining polygon branches sits with decomposing lines generally (probably always do it?)

mdsumner commented 8 years ago

basics of lines in devel branch, with mesh() generic. Probably use "mesh()" rather than "tri_mesh" and other specifics.

replace "globe" as a plotting type to a vertex conversion that just converts projected or longlat to XYZ. Use plotglobe as an extra function?

map_table should add branch topology to output, or it can be inferred from branch "island" property?

mdsumner commented 7 years ago

This suggests a speedup for line primitives, so all the decompositions can be done from Spatial to lines in one hit. It might make sense to use over to classify triangles to branches/objects, rather than loop at all.

f <- "C:/data/rangl/trails.rds"
library(tibble)
library(dplyr)
trails <- readRDS(f)

x <- trails %>% filter(index < 20) %>% 
  transmute(order_ = row_number(), x_ = x, y_ = y, branch_ = id, id = as.integer(factor(paste(x_, y_, sep = "-"))))

v <- x %>% select(x_, y_, id) %>% mutate(vertex_ = id) %>% distinct(vertex_, .keep_all = TRUE)
bXv <- x %>% transmute(vertex_ = id, branch_ = branch_)
b <- bXv %>% distinct(branch_) %>% mutate(object_ = branch_)
o <- tibble(object_ = b$object_)
library(rangl)
rangl0 <- function(tabs, ...) {
  pr4 <- "+proj=stere +lat_0=-90 +lat_ts=-60 +lon_0=180 +x_0=0 +y_0=0 +a=6367470 +b=6367470 +units=m +no_defs "
  ll <- vector("list", nrow(tabs$o))
  for (i_obj in seq(nrow(tabs$o))) {
    tabs_i <- tabs; tabs_i$o <- tabs_i$o[i_obj, ]
    tabs_i <- rangl:::semi_cascade(tabs_i)
    tt_i <- rangl:::line_mesh_map_table1(tabs_i)
    ll[[i_obj]] <- tt_i
  }

  outlist <- vector("list", length(ll[[1]]))
  nms <- names(ll[[1]])
  names(outlist) <- nms
  for (i in seq_along(outlist)) {
    outlist[[i]] <- dplyr::bind_rows(lapply(ll, "[[", nms[i]))
  }

  ## renormalize the vertices
  allverts <- dplyr::inner_join(outlist$lXv, outlist$v, "vertex_")
  allverts$uvert <- as.integer(factor(paste(allverts$x_, allverts$y_, sep = "_")))
  allverts$vertex_ <- spbabel:::id_n(length(unique(allverts$uvert)))[allverts$uvert]
  outlist$lXv <- allverts[, c("segment_", "vertex_")]
  outlist$v <- dplyr::distinct_(allverts, "x_", "y_", "vertex_")
  ## finally add longitude and latitude
  outlist$meta <- tibble::tibble(proj = pr4, x = "x_", y = "y_", ctime = format(Sys.time(), tz = "UTC"))
  class(outlist) <- "linemesh"
  outlist
}

aa <- rangl0(list(o = o, b = b, bXv = bXv, v = v))
plot(aa)
mdsumner commented 7 years ago

I tried this one-hit approach, and I messed up the object colours. (Why??)

library(maptools)
library(dplyr)
data(wrld_simpl)

x <- wrld_simpl[1:100, ]
pr4 <- proj4string(x)
tabs <- spbabel::map_table(x)

tabs$v <- tabs$v %>% mutate(countingIndex = row_number())
## join the vertex-instances to the vertices table
## so, i.e. expand out the duplicated x/y coordinates
nonuq <- dplyr::inner_join(tabs$bXv, tabs$v, "vertex_")

## create Triangle's Planar Straight Line Graph
## which is an index matrix S of every vertex pair P
ps <- RTriangle::pslg(P = as.matrix(tabs$v[, c("x_", "y_")]),
                      S = do.call(rbind, lapply(split(nonuq, nonuq$branch_),
                                                function(x) rangl:::path2seg(x$countingIndex))))
max_area <- NULL
## build the triangulation, with input max_area (defaults to NULL)
tr <- RTriangle::triangulate(ps, a = max_area)
## centroid of every triangle
centroids <- matrix(unlist(lapply(split(tr$P[t(tr$T), ], rep(seq(nrow(tr$T)), each = 3)), .colMeans, 3, 2)), 
                    ncol = 2, byrow = TRUE)

x$countingIndex <- seq(nrow(x))
#tri_index <- over(SpatialPoints(centroids, proj4string = CRS(proj4string(x))), x)$countingIndex
## add triangle topology
tabs$t <- tibble::tibble(triangle_ = spbabel:::id_n(nrow(tr$T)), 
                         object_ = tabs$o$object_[over(SpatialPoints(centroids, proj4string = CRS(proj4string(x))), x)$countingIndex])

#%>% 
#  filter(!is.na(object_))

tabs$v <- tibble::tibble(x_ = tr$P[,1], y_ = tr$P[,2], vertex_ = spbabel:::id_n(nrow(tr$P)))
tabs$b <- tabs$bXv <- NULL

tabs$tXv <- tibble::tibble(triangle_ = rep(tabs$t$triangle_, each = 3), 
                           vertex_ = tabs$v$vertex_[as.vector(t(tr$T))])

tabs$tXv < tabs$tXv[tabs$tXv$triangle_ %in% tabs$t$triangle_, ]

tabs$t <- tabs$t %>% filter(!is.na(object_))
tabs$meta <- tibble::tibble(proj = pr4, x = "x_", y = "y_", ctime = format(Sys.time(), tz = "UTC"))
class(tabs) <- "trimesh"
plot(tabs)
mdsumner commented 6 years ago

silicate is the place for this