build sf

@param gm a gibble, geom map @param vertex_pool the vertices to use, described by the gibble geom map `gm``

build_sf <- function(gm, coords_in, crs = NULL) {
  glist <- vector("list", length(unique(gm$object)))
  coords_in <- gm %>% dplyr::select(-type, -ncol, -nrow) %>%
    dplyr::slice(rep(seq_len(nrow(gm)), gm$nrow)) %>% dplyr::bind_cols(coords_in)
  ufeature <- unique(gm$object)
  for (ifeature in seq_along(ufeature)) {
    gm0 <- gm %>% dplyr::filter(object == ufeature[ifeature])
    type <- gm0$type[1]
    coord0 <- coords_in %>% dplyr::filter(object == ifeature)
    ## object becomes sub-feature element (not a hole, that is "part")
    coord0$object <- rep(seq_len(nrow(gm0)), gm0$nrow)
    glist[[ifeature]] <- switch(type,
                                POINT = sf::st_point(unlist(coord0 %>% dplyr::select(x_, y_))),
                                MULTIPOINT = sf::st_multipoint(as.matrix(coord0 %>% dplyr::select(x_, y_))),
                                LINESTRING = sf::st_linestring(as.matrix(coord0 %>% dplyr::select(x_, y_))),
                                MULTILINESTRING = sf::st_multilinestring(lapply(split(coord0 %>% dplyr::select(x_, y_), coord0$object), as.matrix)),
                                POLYGON = sf::st_polygon(lapply(split(coord0 %>% dplyr::select(x_, y_), coord0$object), as.matrix)),
                                MULTIPOLYGON = sf::st_multipolygon(lapply(split(coord0 %>% dplyr::select(x_, y_, part), coord0$object),
                                                                          function(part) lapply(split(part %>% select(x_, y_), part$part), as.matrix)))
    )
  }
   if (is.null(crs)) crs <- NA_crs_
  out <-   st_sfc(glist, crs = crs)
  out
}


prepare_sf_ct <- function(x) {
  tabs <- sc::PRIMITIVE(x)
  
  segment <-  tibble::tibble(vertex_ = c(t(as.matrix(tabs$segment %>% 
    dplyr::select(.vertex0, .vertex1))))) %>%
    dplyr::inner_join(tabs$vertex %>% 
    dplyr::mutate(vertex = row_number() - 1)) %>% 
    dplyr::mutate(segment = (row_number() + 1) %/% 2)
  segs <- split(segment$vertex, segment$segment)
  
  list(coords = cbind(tabs$vertex$x_, tabs$vertex$y_), segs = distinct_uord_segments(segs))
}

distinct_uord_segments <- function(segs) {
  x <- dplyr::distinct(tibble::as_tibble(do.call(rbind, segs)))
  usort <- do.call(rbind, lapply(segs, sort))
  bad <- duplicated(usort)
  x <- x[!bad, ]
  lapply(split(x, seq_len(nrow(x))), unlist)
}

st_line_from_segment <- function(segs, coords) {
  sf::st_sfc(lapply(segs, function(a) sf::st_linestring(coords[a + 1, ])))
}

build_triangles <- function(coords, segments, ...) {
  RTriangle::triangulate(RTriangle::pslg(P = coords,S = do.call(rbind, segments)+1 ), ...)
}
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(scsf)
## Loading required package: sc
psf <- prepare_sf_ct(rmapshaper::ms_simplify(rnaturalearth::ne_countries(returnclass = "sf")))
## Warning in eval(call(as_fun, df[[n]])): NAs introduced by coercion

## Warning in eval(call(as_fun, df[[n]])): NAs introduced by coercion

## Warning in eval(call(as_fun, df[[n]])): NAs introduced by coercion

## Warning in eval(call(as_fun, df[[n]])): NAs introduced by coercion

## Warning in eval(call(as_fun, df[[n]])): NAs introduced by coercion

## Warning in eval(call(as_fun, df[[n]])): NAs introduced by coercion

## Warning in eval(call(as_fun, df[[n]])): NAs introduced by coercion

## Warning in eval(call(as_fun, df[[n]])): NAs introduced by coercion

## Warning in eval(call(as_fun, df[[n]])): NAs introduced by coercion

## Warning in eval(call(as_fun, df[[n]])): NAs introduced by coercion

## Warning in eval(call(as_fun, df[[n]])): NAs introduced by coercion

## Warning in eval(call(as_fun, df[[n]])): NAs introduced by coercion
## Joining, by = "vertex_"
#data("wrld_simpl", package = "maptools")
library(sf)
## Linking to GEOS 3.5.1, GDAL 2.2.1, proj.4 4.9.3
#psf <- prepare_sf_ct(rmapshaper::ms_simplify(st_as_sf(wrld_simpl)))
tri <- build_triangles(psf$coords, psf$segs, a = 3)

library(rgl)
## Warning in rgl.init(initValue, onlyNULL): RGL: unable to open X11 display
## Warning: 'rgl_init' failed, running with rgl.useNULL = TRUE
xyz <- proj4::ptransform(cbind(tri$P, 0) * pi/180, src.proj  = "+init=epsg:4326", 
                         dst.proj = "+proj=geocent")
dim(xyz)
## [1] 5980    3
rgl.clear()
rgl.triangles(xyz[t(tri$T), ])
rglwidget()
library(marmap)
## 
## Attaching package: 'marmap'
## The following object is masked from 'package:grDevices':
## 
##     as.raster
data("hawaii", package = "marmap")
qm <- quadmesh::quadmesh(raster::aggregate(as.raster(hawaii), fact = 5))
qm$vb[1:2, ] <- t(proj4::project(t(qm$vb[1:2, ]), "+proj=laea +lon_0=180"))
rgl.clear()
shade3d(qm, col = "white")
aspect3d(1, 1, .5)
axes3d()
rglwidget()
library(gibble)
nc <- read_sf(system.file("shape/nc.shp", package="sf"))
nc_gmap <- gibble(nc)
library(dplyr)
coords <- sf::st_coordinates(nc) %>% as_tibble() %>% dplyr::transmute(x_ = X, y_ = Y)
anc <- build_sf(nc_gmap, coords, crs = st_crs(nc))
anc
## Geometry set for 108 features 
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## epsg (SRID):    4267
## proj4string:    +proj=longlat +datum=NAD27 +no_defs
## First 5 geometries:
## MULTIPOLYGON (((-81.4727554321289 36.2343559265...
## MULTIPOLYGON (((-81.2398910522461 36.3653640747...
## MULTIPOLYGON (((-80.4563446044922 36.2425575256...
## MULTIPOLYGON (((-76.0089721679688 36.3195953369...
## MULTIPOLYGON (((-76.0271682739258 36.5567169189...