geovctrs for TRI

A simple prototype to take silicate/anglr to geovctrs

  library(geovctrs)
  library(vctrs)
## triangles, taken from https://paleolimbot.github.io/geovctrs/articles/extending-geovctrs.html
# base constructor for developers
new_triangle <- function(x = list(x0 = geo_xy(),
                                  x1 = geo_xy(),
                                  x2 = geo_xy())) {
  vec_assert(x$x0, geo_xy())
  vec_assert(x$x1, geo_xy())
  vec_assert(x$x2, geo_xy())

  vctrs::new_rcrd(x, class = "triangle")
}

# convenient constructor for users
triangle <- function(x0 = geo_xy(),
                     x1 = geo_xy(),
                     x2 = geo_xy()) {
  new_triangle(
    vec_recycle_common(
      x0 = vec_cast(x0, geo_xy()),
      x1 = vec_cast(x1, geo_xy()),
      x2 = vec_cast(x2, geo_xy())
    ) )
}

triangle_TRI <- function(x) {
  if (inherits(x, "DEL")) {
    x$triangle <- dplyr::filter(x$triangle, .data$visible)
  }
  idx <- do.call(rbind, sc_object(x)$topology_)[c(".vx0", ".vx1", ".vx2")]
  XX <- sc_vertex(x)$x_
  YY <- sc_vertex(x)$y_
  triangle(x0 = geo_xy(XX[idx$.vx0], YY[idx$.vx0]),
           x1 = geo_xy(XX[idx$.vx1], YY[idx$.vx1]),
           x2 = geo_xy(XX[idx$.vx2], YY[idx$.vx2]))
}
# rcrd vectors require a format method
format.triangle <- function(x, ...) {
  sprintf("<triangle: %s->%s->%s>",
          format(field(x, "x0")),
          format(field(x, "x1")),
          format(field(x, "x2")))
}


triangle_coordinates <- function(x) {
   c(
      field(x, "x0"),
      field(x, "x1"),
      field(x, "x2")
    )
}

#The easiest way to construct a geovctr is using the geo_collection() family of functions, all of which take coordinates and return a geometry vector. For a single circle, we could construct polygons like so:

triangle_polygon <- function(x) {
    geo_polygon(triangle_coordinates(x)[c(1, 2, 3, 1)])
}


as_geovctr.triangle <- function(x, ...) {
    x_data <- x
    features <- lapply(
      seq_len(vec_size(x)),
      function(i) triangle_polygon(x_data[i])
    )

    vec_c(!!!features)
}

# create a vector of triangles
x <- triangle(geo_xy(0:5, 0:5),
              geo_xy(2:7, 1:6),
              geo_xy(1:6, 2:7))
x
#> <triangle[6]>
#> [1] <triangle: (0 0)->(2 1)->(1 2)> <triangle: (1 1)->(3 2)->(2 3)>
#> [3] <triangle: (2 2)->(4 3)->(3 4)> <triangle: (3 3)->(5 4)->(4 5)>
#> [5] <triangle: (4 4)->(6 5)->(5 6)> <triangle: (5 5)->(7 6)->(6 7)>

triangle_coordinates(x)
#> <geovctrs_xy[18]>
#>  [1] (0 0) (1 1) (2 2) (3 3) (4 4) (5 5) (2 1) (3 2) (4 3) (5 4) (6 5) (7 6)
#> [13] (1 2) (2 3) (3 4) (4 5) (5 6) (6 7)

triangle_polygon(x)
#> <geovctrs_collection[1]>
#> [1] POLYGON (0 0)...+3


plot(as_geovctr(x))

geo_envelope(x)
#> <geovctrs_rect[6]>
#> [1] (0 0...2 2) (1 1...3 3) (2 2...4 4) (3 3...5 5) (4 4...6 6) (5 5...7 7)
geo_plot(x, lty = 1:6)

## draft conversion from silicate/anglr
library(silicate)
#> 
#> Attaching package: 'silicate'
#> The following object is masked from 'package:stats':
#> 
#>     filter
geo_plot(triangle_TRI(TRI0(silicate::minimal_mesh)))

^{Created on 2020-04-20 by the reprex package (v0.3.0)}

A fast mesh_plot for vectorized coloured triangles for geovctrs_triang

mesh_plot.triangle <- function(x, col = NA, border = NA, ..., add = FALSE) {
  x0 <- field(x, "x0")
  x1 <- field(x, "x1")
  x2 <- field(x, "x2")
  xx  <- list(x = c(t(cbind(field(x0, "x"), field(x1, "x"), field(x2, "x")))), 
              y = c(t(cbind(field(x0, "y"), field(x1, "y"), field(x2, "y")))), 
              id = rep(seq_along(x), each = 3L))
if (!add) {
  graphics::plot.new()
  graphics::plot.window(xlim = range(xx$x, finite = TRUE), ylim = range(xx$y, finite = TRUE),
                        ...)
}
vps <- gridBase::baseViewports()

grid::pushViewport(vps$inner, vps$figure, vps$plot)

if (!is.na(col[1])) {
  xx$col <- rep(col, length.out = length(x))
}
grid::grid.polygon(xx$x, xx$y, xx$id, gp = grid::gpar(col = border, fill = xx$col),
                   default.units = "native")

grid::popViewport(3)

}


x <- triangle_TRI(DEL0(minimal_mesh, max_area = 0.0001))
system.time(geo_plot(x[1:100]))
   user  system elapsed 
   0.76    0.04    0.86 
system.time(mesh_plot(x[order(x)], col = viridis::viridis(length(x))))
   user  system elapsed 
   0.13    0.11    0.24