fusion

[mdsumner]

break and re-mesh by combining two together

[mdsumner]

Fakey for now, need to "interpolate fused vertex attributes" for the elevation.

library(rgdal)
dp <- "."
dsn1 <- file.path(dp, "LIST_PARCELS_HOBART/list_parcels_hobart.gdb")
dsn2 <- file.path(dp, "LIST_5M_CONTOURS_HOBART/list_5m_contours_hobart.gdb")
layer1 <- "list_parcels_hobart"
layer2 <- "list_5m_contours_hobart"
cad0 <- readOGR(dsn1, layer1)
cont0 <- readOGR(dsn2, layer2 )

library(raster)
# p <- as(extent(cad), "SpatialPolygons")
# projection(p) <- projection(cad)
# library(dismo)
# 
# gm <- gmap(p, type = "satellite")

e <- new("Extent"
         , xmin = 523391.060663871
         , xmax = 524529.793420587
         , ymin = 5249817.21778203
         , ymax = 5250452.78955322
)

cad <- crop(cad0, e)
cont <- crop(cont0, e)

library(rangl)

x1<- mesh(cad)
x2 <- mesh(cont)

x1 <- rangl:::mesh.SpatialLines(cad)
x2 <- mesh(cont)
vv <- x2$v
library(dplyr)
vv <- right_join(vv, 
                 inner_join(inner_join(x2$o[, c("ELEVATION", "object_")], x2$l), x2$lXv)[, c("ELEVATION", "vertex_")])

vv <- distinct_(vv, c("vertex_"), .keep_all = TRUE)
vv <- mutate(vv, z_ = ELEVATION)
vv$ELEVATION <- NULL
## combine all vertices and edges and triangulate
library(dplyr)
v <- 
  bind_rows(x1$v, vv)
v$z_[is.na(v$z_)] <- 50
v$countingIndex <- seq(nrow(v))
lXv <- bind_rows(x1$lXv, x2$lXv)
ps <- RTriangle::pslg(P = cbind(v$x_, v$y_), 
                      S = matrix(inner_join(lXv, v)$countingIndex, ncol = 2, byrow = TRUE), 
                      PA = cbind(v$z_))

tr <- RTriangle::triangulate(ps)

## rebuild
x <- structure(list(o = tibble(object_ = 1), 
               t = tibble(object_ = 1, triangle_ = seq(nrow(tr$T))), 
               tXv = tibble(triangle_ = rep(seq(nrow(tr$T)), each = 3), vertex_ = as.vector(t(tr$T))), 
               v = tibble(x_ = tr$P[,1], y_ = tr$P[,2], z_ = tr$PA[,1], vertex_=  seq(nrow(tr$P)))
               ), class = "trimesh")

plot(x, col = viridis::viridis(nrow(x$v)))

[mdsumner]

This version does a nn match for Z, not sure if a trace-edge approach or a neighbourhood approach is best here. Also need to mesh in each object separately, so we can keep track of their identity. Somehow a generalization of mesh line and mesh poly is key here. This interpolation is craap, but where it goes whack is where the most impact of detailed cadastral parcels is meshed in - and it's not the biggest problem still - keeping the parcel objs is the crux, and generalizing to 3 or more inputs, and maintaining those links to sources . . .

library(rgdal)
dp <- "."
dsn1 <- file.path(dp, "LIST_PARCELS_HOBART/list_parcels_hobart.gdb")
dsn2 <- file.path(dp, "LIST_5M_CONTOURS_HOBART/list_5m_contours_hobart.gdb")
layer1 <- "list_parcels_hobart"
layer2 <- "list_5m_contours_hobart"
cad0 <- readOGR(dsn1, layer1)
cont0 <- readOGR(dsn2, layer2 )

library(raster)
# p <- as(extent(cad), "SpatialPolygons")
# projection(p) <- projection(cad)
# library(dismo)
# 
# gm <- gmap(p, type = "satellite")

e <- new("Extent"
         , xmin = 523391.060663871
         , xmax = 524529.793420587
         , ymin = 5249817.21778203
         , ymax = 5250452.78955322
)

cad <- crop(cad0, e)
cont <- crop(cont0, e)

library(rangl)

x1<- mesh(cad)
x2 <- mesh(cont)

x1 <- rangl:::mesh.SpatialLines(cad)
x2 <- mesh(cont)
vv <- x2$v
library(dplyr)
vv <- right_join(vv, 
                 inner_join(inner_join(x2$o[, c("ELEVATION", "object_")], x2$l), x2$lXv)[, c("ELEVATION", "vertex_")])

vv <- distinct_(vv, c("vertex_"), .keep_all = TRUE)
vv <- mutate(vv, z_ = ELEVATION)
vv$ELEVATION <- NULL
## combine all vertices and edges and triangulate
## this part really needs to be done per object, so we can
## keep track of what belongs to what
library(dplyr)
v <- 
  bind_rows(x1$v, vv)
#v$z_[is.na(v$z_)] <- 50
v$countingIndex <- seq(nrow(v))
lXv <- bind_rows(x1$lXv, x2$lXv)
ps <- RTriangle::pslg(P = cbind(v$x_, v$y_), 
                      S = matrix(inner_join(lXv, v)$countingIndex, ncol = 2, byrow = TRUE))

tr <- RTriangle::triangulate(ps)

## rebuild
x <- structure(list(o = tibble(object_ = 1), 
               t = tibble(object_ = 1, triangle_ = seq(nrow(tr$T))), 
               tXv = tibble(triangle_ = rep(seq(nrow(tr$T)), each = 3), vertex_ = as.vector(t(tr$T))), 
               v = tibble(x_ = tr$P[,1], y_ = tr$P[,2],  vertex_=  seq(nrow(tr$P)))
               ), class = "trimesh")

nn <- nabor::knn(vv[, c("x_", "y_")], x$v[, c("x_", "y_")], k = 1, eps = 0)
x$v$z_ <- vv$z_[nn$nn.idx]
#x$v$z_ <- apply(matrix(vv$z_[nn$nn.idx], nrow = 5), 2, median, na.rm = TRUE)

scl <- function(x) (x - min(x)) / diff(range(x))

#tab <- x$t %>% inner_join(x$tXv) %>% inner_join(x$v) %>% group_by(triangle_) %>% summarize(z = mean(z_))
zm <- apply(matrix(a$vb[3,a$it], nrow = 3), 2, mean)
a <- plot(x, col = rep(viridis::inferno(29)[scl(zm) * 28 + 1], each = 3))

## put the lines on and add to the plot
x2$v <- vv
x2$o$color_ = "black"
plot(x2)

[mdsumner]

Topology overlay (identity), crux is to triangulate all edges and then filter triangles based on pip of their centroids.

library(raster)
library(spex)
p <- shapefile(system.file("external/lux.shp", package="raster"))
b <- spex(extent(6, 6.4, 49.75, 50), projection(p))

library(rangl)

x1<- mesh(as(p, "SpatialLinesDataFrame"))
x2 <- mesh(as(b, "SpatialLinesDataFrame"))

## combine all vertices and edges and triangulate
## this part really needs to be done per object, so we can
## keep track of what belongs to what
library(dplyr)
v <- 
  bind_rows(x1$v, x2$v)

v$countingIndex <- seq(nrow(v))
lXv <- bind_rows(x1$lXv, x2$lXv)
ps <- RTriangle::pslg(P = cbind(v$x_, v$y_), 
                      S = matrix(inner_join(lXv, v)$countingIndex, ncol = 2, byrow = TRUE))

tr <- RTriangle::triangulate(ps)

centroids <- t(apply(tr$T, 1, function(x) apply(tr$P[x, ], 2, mean)))
## x1 has this triangle
x1_pip <- !is.na(extract(p, centroids)$poly.ID)
## x2 has this triangle
x2_pip <- !is.na(extract(b, centroids)$poly.ID)

shade3d(structure(list(vb = t(cbind(tr$P, 0, 1)), it = t(tr$T[x1_pip, ])), class = c("mesh3d", "shape3d")), col = "firebrick")
shade3d(structure(list(vb = t(cbind(tr$P, 1, 1)), it = t(tr$T[x2_pip, ])), class = c("mesh3d", "shape3d")), border = "black", col = "dodgerblue")
lines3d(cbind(tr$P, 0.5)[t(tr$E), ])
aspect3d(1, 1, 0.1)

[mdsumner]

Previously done here, which will be worth comparing with the stuff above https://github.com/r-gris/cgalgris/blob/master/inst/topoverlay.R

[mdsumner]

This can only work with edge-based triangulations, which is handy because we can keep RTriangle here and let silicate stay pure with decido.

[mdsumner]

A fuse() function will operate like a generalized copy_down(), taking any number of objects and converting all to the lowest level model required.

It could take multiple objects but also a form with expressions would be good,il at least for constant object value copy down.

[mdsumner]

A much firmer basis for this now, key things are to combine two SC objects and set as one "object_", then re-normalize the edge-vertex map. Triangle sorts out the intersections.

library(raster)
library(anglr)
library(silicate)

library(raster)
library(spex)
p1 <- SC(shapefile(system.file("external/lux.shp", package="raster")))
p2 <- SC(spex::spex(extent(6, 6.4, 49.75, 50), p1$meta$proj))

#' Basic fusion
#' 
#' This just gives the plane-partition, unique triangles that include
#' every input edge. Still requires each triangle to be identified in
#' each input layer so we can differentiate them (or determine intersection etc.)
fuse0 <- function(p1, p2) {
  c_SC <- function(x, y) {
    structure(purrr::map2(x, y, dplyr::bind_rows), class = c("SC", "sc"))
  }

  x <- c_SC(p1, p2)
  x$object <- x$object[1, ]
  x$object_link_edge$object_ <- x$object$object_[1]
  ## re-normalize
  eXv <- x$edge %>% tidyr::gather(end_point, vertex_, -edge_) %>% 
    dplyr::inner_join(x$vertex)
  idx <- dplyr::group_indices(eXv, x_, y_)
  eXv$vertex_ <- sc_uid(length(unique(idx)))[idx]
  x$edge <- eXv %>% dplyr::select(edge_, end_point, vertex_) %>% tidyr::spread(end_point, vertex_)
  x$vertex <- dplyr::distinct(eXv, x_, y_, vertex_)
  DEL(x)  
}

x <- fuse0(p1, p2)

plot(x)  ## only one triangle, but might belong to both inputs

[mdsumner]

This is getting closer, taken from accepted answer https://gis.stackexchange.com/questions/140504/extracting-intersection-areas-in-r

library(sp)
library(raster)
library(rgdal)
library(rgeos)
library(maptools)

# Example data from raster package
p1 <- shapefile(system.file("external/lux.shp", package="raster"))
# Remove attribute data
p1 <- as(p1, 'SpatialPolygons')
# Add in some fake soil type data
soil <- SpatialPolygonsDataFrame(p1, data.frame(soil=LETTERS[1:12]), match.ID=F)

# Field polygons
p2 <- union(as(extent(6, 6.4, 49.75, 50), 'SpatialPolygons'),
            as(extent(5.8, 6.2, 49.5, 49.7), 'SpatialPolygons'))
field <- SpatialPolygonsDataFrame(p2, data.frame(field=c('x','y')), match.ID=F)
projection(field) <- projection(soil)

# intersect from raster package
pi <- intersect(soil, field)
plot(soil, axes=T); plot(field, add=T); plot(pi, add=T, col='red')

# Extract areas from polygon objects then attach as attribute
pi$area <- area(pi) / 1000000

# For each field, get area per soil type
aggregate(area~field + soil, data=pi, FUN=sum)

############################################################
############################################################

library(raster)
library(anglr)
library(silicate)

library(raster)
library(spex)
p1 <- SC(soil)
p2 <- SC(field)

#' Basic fusion
#' 
#' fuse0 just gives the plane-partition, unique triangles that include
#' every input edge. Still requires each triangle to be identified in
#' each input layer so we can differentiate them (or determine intersection etc.)
fuse0 <- function(p1, p2) {
  c_SC <- function(x, y) {
    structure(purrr::map2(x, y, dplyr::bind_rows), class = c("SC", "sc"))
  }
  x <- c_SC(p1, p2)
  x$object <- x$object[1, ]
  x$object_link_edge$object_ <- x$object$object_[1]
  ## re-normalize
  eXv <- x$edge %>% tidyr::gather(end_point, vertex_, -edge_) %>% 
    dplyr::inner_join(x$vertex)
  idx <- dplyr::group_indices(eXv, x_, y_)
  eXv$vertex_ <- sc_uid(length(unique(idx)))[idx]
  x$edge <- eXv %>% dplyr::select(edge_, end_point, vertex_) %>% tidyr::spread(end_point, vertex_)
  x$vertex <- dplyr::distinct(eXv, x_, y_, vertex_)
  DEL(x)  
}
triverts <- function(model, ...) {
  tidyr::gather(model$triangle, node, vertex_, -triangle_, -visible) %>% 
    dplyr::inner_join(model$vertex)
}
## unnecessarily slow, needs bbox-filter and point-in-triangle lookup
triindex <- function(triv, centroids) {
  idx <- triv %>% split(.$triangle_) %>% 
    purrr::map(~which(sp::point.in.polygon(centroids$x_, centroids$y_, 
                                           .x$x_, .x$y_) == 1L)) 
  tibble::tibble(tri = rep(names(idx), lengths(idx)), 
                 mesh = centroids$triangle_[unlist(idx)])
}
fuse <- function(p1, p2) {
  mrg <- fuse0(p1, p2)
  t1 <- DEL(p1)
  t2 <- DEL(p2)
  mesh_triverts <- triverts(mrg)
  mesh_centroids <- mesh_triverts %>% 
    dplyr::group_by(triangle_) %>% 
    dplyr::summarize(x_ = mean(x_), y_ = mean(y_))  ## not the fastest triangle centre calculator ...
  ## t1 triangles
  t1_triverts <- triverts(t1)
  t2_triverts <- triverts(t2)
  ## layer triangles (tri) with member mesh triangles (mesh)
  t2_link <- triindex(t2_triverts, mesh_centroids) %>% 
      dplyr::inner_join(t2$object_link_triangle, c("tri" = "triangle_")) %>% 
      dplyr::transmute(triangle_ = .data$mesh, object_ = .data$object_) %>% 
      dplyr::distinct(triangle_, object_)
  t1_link <- triindex(t1_triverts, mesh_centroids) %>% 
    dplyr::inner_join(t1$object_link_triangle, c("tri" = "triangle_")) %>% 
  dplyr::transmute(triangle_ = .data$mesh, object_ = .data$object_) %>% 
    dplyr::distinct(triangle_, object_)

  structure(list(object1 = t1$object, object2 = t2$object, object_link_triangle = rbind(t1_link, t2_link), 
            triangle = mrg$triangle, vertex = mrg$vertex, meta = mrg$meta), class = c("FUSION", "DEL", "TRI", "sc"))
}

x <- fuse(p1, p2)

## simulate activate
x$object <- x$object1

par(mfrow   = c(1, 2))
plot(x$vertex[c("x_", "y_")], type = "n", asp = 1)
plot(x, add = TRUE)
x$object <- x$object2
plot(x$vertex[c("x_", "y_")], type = "n", asp = 1)
plot(x, add = TRUE) 

[mdsumner]

This is probably super-easy in an all-sf context, but we need collections of geometries (same as GC) to respect all eges as a set, not per-feature - then we could just use sfdct with sfg-only, use centroids and grouping to link back up.

[mdsumner]

This should be easy from the SC perspective. Need to have a way to merge SC models, deal with implicit z as missing (not 0) or prioritize some model z over others. RTriangle can deal with the interpolation of fused vertex geometry, maybe a DEL(....) space bucket kinda thing

[mdsumner]

Mesh fusion, with leftover #78

dp <- "/rdsi/PUBLIC/raad/data/listdata.thelist.tas.gov.au/opendata/data"
dsn1 <- file.path(dp, "list_parcels_hobart.gdb")
dsn2 <- file.path(dp, "list_5m_contours_hobart.gdb")
layer1 <- "list_parcels_hobart"
layer2 <- "list_5m_contours_hobart"
library(sf)
cad0 <- read_sf(dsn1, layer1)
cont0 <- read_sf(dsn2, layer2 )

library(raster)
# p <- as(extent(cad), "SpatialPolygons")
# projection(p) <- projection(cad)
# library(dismo)
# 
# gm <- gmap(p, type = "satellite")

e <- new("Extent"
         , xmin = 523391.060663871
         , xmax = 524529.793420587
         , ymin = 5249817.21778203
         , ymax = 5250452.78955322
)

bb <- st_bbox(e)
cad <- st_crop(cad0, bb)
cont <- st_cast(st_crop(cont0,bb), "MULTILINESTRING")

#rownames(cad) <- sc_uid(nrow(cad))
#rownames(cont) <- sc_uid(nrow(cont))
library(silicate)
library(anglr)
cad_sc <- SC(cad)
cont_sc <- copy_down(SC(cont), "ELEVATION")

## merge edge sets
library(dplyr)
x <- cad_sc
x$object <- bind_rows(x$object["object_"], cont_sc$object["object_"])
x$object_link_edge <- bind_rows(x$object_link_edge, cont_sc$object_link_edge)
x$edge <- bind_rows(x$edge, cont_sc$edge)
x$vertex <- bind_rows(x$vertex, cont_sc$vertex)
## remap vertices
uj <- unjoin::unjoin(x$vertex, x_, y_, key_col = "old")

newv <- uj$old %>% mutate(vertex_ = sc_uid(uj$old))
mapp <- inner_join(uj$old, uj$data) %>% arrange(old)

x$edge$.vx0 <- newv$vertex_[match(mapp$old[match(x$edge$.vx0, mapp$vertex_)], uj$old$old)]
x$edge$.vx1 <- newv$vertex_[match(mapp$old[match(x$edge$.vx1, mapp$vertex_)], uj$old$old)]
newv$old <- NULL
x$vertex <- newv
plot(x) ## check

## all the empty set
# str(setdiff(x$object$object_, x$object_link_edge$object_))
# str(setdiff(x$object_link_edge$edge_, x$edge$edge_))
# str(setdiff(c(x$edge$.vx0, x$edge$.vx1), x$vertex$vertex_))

y <- DEL(x)
# Warning messages:
#   1: In max(TRIS[[i]]) : no non-missing arguments to max; returning -Inf
#...

plot(y)
#Error in left_join_impl(x, y, by_x, by_y, aux_x, aux_y, na_matches) : 
#  Can't join on 'vertex_' x 'vertex_' because of incompatible types (character / logical)

[mdsumner]

Oh, we can't just rbind objects together - the fused data doesn't have objects properly, it has links to the original objects. We can produce new objects, but we have to commit - so the mesh is joined by pure vertex+edges, and vertex de-duplication, and then we have to decide what object means. In this plot, it's clear that the property boundaries are winning, because the line objects from the input have no say in the triangle's interior.

[mdsumner]

Ah no, only the poly mesh is all there:

plot(y$vertex$x_, y$vertex$y_, pch = ".")
polygon(y$vertex[match(t(cbind(as.matrix(y$triangle[c(".vx0", ".vx1", ".vx2")]), NA)), y$vertex$vertex_), c("x_", "y_")])

this is where DEL is unsuitable, we need to drop all objects.

[mdsumner]

MOOO, just assign as a single object!

library(silicate)
library(anglr)
cad_sc <- SC(cad)
cont_sc <- copy_down(SC(cont), "ELEVATION")

## merge edge sets
library(dplyr)
x <- cad_sc
x$object <- tibble(object_ = "a_object")
x$object_link_edge <- bind_rows(x$object_link_edge, cont_sc$object_link_edge)
x$object_link_edge$object_ <- "a_object"
x$edge <- bind_rows(x$edge, cont_sc$edge)
x$vertex <- bind_rows(x$vertex, cont_sc$vertex)
## remap vertices
uj <- unjoin::unjoin(x$vertex, x_, y_, key_col = "old")

newv <- uj$old %>% mutate(vertex_ = sc_uid(uj$old))
mapp <- inner_join(uj$old, uj$data) %>% arrange(old)

x$edge$.vx0 <- newv$vertex_[match(mapp$old[match(x$edge$.vx0, mapp$vertex_)], uj$old$old)]
x$edge$.vx1 <- newv$vertex_[match(mapp$old[match(x$edge$.vx1, mapp$vertex_)], uj$old$old)]
newv$old <- NULL
x$vertex <- newv
plot(x) ## check

# str(setdiff(x$object$object_, x$object_link_edge$object_))
# str(setdiff(x$object_link_edge$edge_, x$edge$edge_))
# str(setdiff(c(x$edge$.vx0, x$edge$.vx1), x$vertex$vertex_))

y <- DEL(x)

plot(y)

This is now the spacebucket, come at from the silicate approach. We can't maintain objects at this level, because some are polygon surfaces and some are lines and some points. The mesh contains everything, and we can query it for primitives that properly link to the inputs

[mdsumner]

related to #117, together will need general handling of the capture-z-copy-over thing with RTriangle P/PA