Merge branch 'devel' of https://github.com/r-gris/rangl into devel

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: rangl
 Type: Package
 Title: Tidy Flexible Spatial
-Version: 0.3.0
+Version: 0.3.0.9001
 Authors@R: person("Michael D.","Sumner", role = c("aut", "cre"), email =
     "[email protected]")
 Description: Create tidy tables for topological, triangulated, and other spatial data structures.
@@ -10,6 +10,7 @@ Imports:
     dplyr,
     methods,
     proj4,
+    raster,
     rgl,
     RTriangle,
     sp,
@@ -26,4 +27,4 @@ SystemRequirements: PROJ.4
 Encoding: UTF-8
 LazyData: true
 RoxygenNote: 5.0.1
-Remotes: mdsumner/spbabel
+Remotes: mdsumner/spbabel, mdsumner/trip

NAMESPACE

@@ -3,7 +3,9 @@
 S3method(globe,default)
 S3method(plot,linemesh)
 S3method(plot,pointmesh)
+S3method(plot,quad_mesh)
 S3method(plot,trimesh)
+S3method(rangl,RasterLayer)
 S3method(rangl,SpatialLines)
 S3method(rangl,SpatialMultiPoints)
 S3method(rangl,SpatialPoints)
@@ -19,6 +21,12 @@ importFrom(RTriangle,triangulate)
 importFrom(dplyr,inner_join)
 importFrom(dplyr,semi_join)
 importFrom(methods,slotNames)
+importFrom(raster,projection)
+importFrom(raster,values)
+importFrom(raster,xmax)
+importFrom(raster,xmin)
+importFrom(raster,ymax)
+importFrom(raster,ymin)
 importFrom(rgl,shade3d)
 importFrom(sp,CRS)
 importFrom(sp,SpatialPoints)
@@ -31,3 +39,4 @@ importFrom(stats,setNames)
 importFrom(tibble,as_tibble)
 importFrom(tibble,tibble)
 importFrom(utils,head)
+importFrom(utils,tail)

NEWS.md

@@ -1,3 +1,15 @@
+# rangl dev
+
+* raster package is now an Import
+
+* added support for RasterLayer
+
+* fixed globe to keep PROJ.4
+
+* quashed a major bug introduced by use of dplyr::distinct, best to use factor unique classifier on character versions of coords
+
+* several cleanup fixes
+
 # rangl 0.3.0
 
 * rename again, main function is called 'rangl', the term 'mesh' is too often used across R 

R/globe.r

@@ -37,5 +37,8 @@ globe.default <- function(x, gproj = "+proj=geocent +ellps=WGS84", ...) {
   x$v$x_ <- xyz[,1]
   x$v$y_ <- xyz[,2]
   x$v$z_ <- xyz[,3]
+  x$meta <- rbind(x$meta[1, ], x$meta)
+  x$meta$proj[1] <- gproj
+  x$meta$ctime[1] <- format(Sys.time(), tz = "UTC")
   x
 }

R/rangl-lines.r

@@ -37,15 +37,15 @@ line_mesh_map_table1 <- function(tabs) {
   tabs$v$countingIndex <- seq(nrow(tabs$v))
   nonuq <- dplyr::inner_join(tabs$bXv, tabs$v, "vertex_")
 
-  ps <- RTriangle::pslg(P = as.matrix(tabs$v[, c("x_", "y_")]),
+  pl <- list(P = as.matrix(tabs$v[, c("x_", "y_")]),
                         S = do.call(rbind, lapply(split(nonuq, nonuq$branch_),
                                                   function(x) path2seg(x$countingIndex))))
 
-  tabs$v <- tibble::tibble(x_ = ps$P[,1], y_ = ps$P[,2], vertex_ = spbabel:::id_n(nrow(ps$P)))
+  tabs$v <- tibble::tibble(x_ = pl$P[,1], y_ = pl$P[,2], vertex_ = spbabel:::id_n(nrow(pl$P)))
   tabs$b <- tabs$bXv <- NULL
-  tabs$l <- tibble::tibble(segment_ = spbabel:::id_n(nrow(ps$S)), object_ = tabs$o$object_[1])
+  tabs$l <- tibble::tibble(segment_ = spbabel:::id_n(nrow(pl$S)), object_ = tabs$o$object_[1])
   tabs$lXv <- tibble::tibble(segment_ = rep(tabs$l$segment_, each = 2), 
-                             vertex_ = tabs$v$vertex_[as.vector(t(ps$S))])
+                             vertex_ = tabs$v$vertex_[as.vector(t(pl$S))])
 
   tabs
 }
@@ -63,9 +63,6 @@ rangl.SpatialLines <- function(x, ...) {
     tabs_i <- tabs; tabs_i$o <- tabs_i$o[i_obj, ]
     tabs_i <- semi_cascade(tabs_i)
     tt_i <- line_mesh_map_table1(tabs_i)
-    # plot.trimesh(tt_i)
-    # scan("", 1L)
-    # rgl::rgl.clear()
     ll[[i_obj]] <- tt_i
   }
 
@@ -83,7 +80,7 @@ rangl.SpatialLines <- function(x, ...) {
   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_")
+  outlist$meta <- tibble::tibble(proj = pr4, x = "x_", y = "y_", ctime = format(Sys.time(), tz = "UTC"))
   class(outlist) <- "linemesh"
   outlist
 }

R/rangl-polygons.r

@@ -1,5 +1,6 @@
 #' @importFrom utils head
 path2seg <- function(x) {
+  ## this is a trick of array logic to generate paired indexes from a sequence
   head(suppressWarnings(matrix(x, nrow = length(x) + 1, ncol = 2, byrow = FALSE)), -2L)
 }
 
@@ -18,34 +19,52 @@ mesh <- function(x, ...) {
   .Deprecated("rangl", package= "rangl", old = "mesh")
 }
 
+## this internal function does the decomposition to primitives of a 
+##  single Spatial object, i.e. a "multipolygon"
+## we need to do it one object at a time otherwise keeping track
+## of the input versus add vertices is harder (but maybe possible later)
 tri_mesh_map_table1 <- function(tabs, max_area = NULL) {
+  ## the row index of the vertices
+  ## we need this in the triangulation
   tabs$v$countingIndex <- seq(nrow(tabs$v))
+  ## 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) path2seg(x$countingIndex))))
 
-  ## FIXME: need to pick sensible behaviour for a
+  ## build the triangulation, with input max_area (defaults to NULL)
   tr <- RTriangle::triangulate(ps, a = max_area)
 
-  ## process the holes if needed
-  ## may be quicker than testing entire object
+  ## process the holes if present
   if (any(!tabs$b$island_)) {
+    ## filter out all the hole geometry and build an sp polygon object with it
+    ## this 
+    ##   filters all the branches that are holes
+    ##   joins on the vertex instance index
+    ##   joins on the vertex values
+    ##   recomposes a SpatialPolygonsDataFrame using the spbabel::sp convention
     holes <- spbabel::sp(dplyr::inner_join(dplyr::inner_join(dplyr::filter_(tabs$b, quote(!island_)), tabs$bXv, "branch_"), 
                                            tabs$v, "vertex_"))
+    ## 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)
-
+    ## sp::over() is very efficient, but has to use high-level objects as input
     badtris <- !is.na(over(SpatialPoints(centroids), sp::geometry(holes)))
+    ## presumably this will always be true inside this block (but should check some time)
     if (any(badtris)) tr$T <- tr$T[!badtris, ]
   }
 
   ## trace and remove any unused triangles
-
+  ## the raw vertices with a unique vertex_ id
   tabs$v <- tibble::tibble(x_ = tr$P[,1], y_ = tr$P[,2], vertex_ = spbabel:::id_n(nrow(tr$P)))
+  ## drop the path topology
   tabs$b <- tabs$bXv <- NULL
-  #tabs$o <- tabs$o[1,]  ## FIX ME
+  ## add triangle topology
   tabs$t <- tibble::tibble(triangle_ = spbabel:::id_n(nrow(tr$T)), object_ = tabs$o$object_[1])
   tabs$tXv <- tibble::tibble(triangle_ = rep(tabs$t$triangle_, each = 3), 
                              vertex_ = tabs$v$vertex_[as.vector(t(tr$T))])
@@ -96,7 +115,7 @@ rangl.SpatialPolygons <- function(x, max_area = NULL, ...) {
   outlist$tXv <- allverts[, c("triangle_", "vertex_")]
   outlist$v <- dplyr::distinct_(allverts,  "vertex_", .keep_all = TRUE)[, c("x_", "y_", "vertex_")]
   ## finally add longitude and latitude
-  outlist$meta <- tibble::tibble(proj = pr4, x = "x_", y = "y_")
+  outlist$meta <- tibble::tibble(proj = pr4, x = "x_", y = "y_", ctime = format(Sys.time(), tz = "UTC"))
   class(outlist) <- "trimesh"
   outlist
 }

R/rangl-raster.r

@@ -0,0 +1,79 @@
+#' Raster rangl
+#' 
+#' Colours not supported, this just gives the virids palette sequentially. 
+#' @param x \code{\link[raster]{raster}}
+#' @param z \code{\link[raster]{raster}}, by default \code{x} is used
+#' @param na.rm remove missing values
+#' @param ... unused
+#' @return quad_mesh
+#' @export
+#' @importFrom raster projection values xmin xmax ymin ymax
+#' @examples
+#' library(raster)
+#' w <- raster(volcano)
+#' plot(rangl(w/300))
+#' 
+rangl.RasterLayer <-  function(x, z = x, na.rm = FALSE, ...) {
+  x <- x[[1]]  ## just the oneth raster for now
+  pr4 <- projection(x)
+  exy <- edges0(x)
+  ind <- apply(prs0(seq(ncol(x) + 1)), 1, p_4, nc = ncol(x) + 1)
+  ## all face indexes
+  ind0 <- as.vector(ind) +
+    rep(seq(0, length = nrow(x), by = ncol(x) + 1), each = 4 * ncol(x))
+  ind1 <- matrix(ind0, nrow = 4)
+  if (na.rm) {
+    ind1 <- ind1[,!is.na(values(x))]
+  }
+  o <- tibble(object_ = 1, xmin = xmin(x), xmax = xmax(x), ymin = ymin(x), ymax = ymax(x), nrow = nrow(x), ncol = ncol(x), proj = projection(x))
+  qXv <- tibble(vertex_ = as.vector(ind1), quad_ = rep(seq(ncol(ind1)), each = 4))
+  if (!is.null(z)) z <- raster::extract(z, exy, method = "bilinear") else z <- 0
+  v <- tibble(x_ = exy[,1], y_ = exy[,2], z_ = z)
+  l <- list(o = o, qd = tibble(quad = seq(ncol(ind1)), object_ = 1), qXv = qXv, v = v)
+  l$meta <- tibble::tibble(proj = pr4, x = "x_", y = "y_", ctime = format(Sys.time(), tz = "UTC"))
+  class(l) <- "quad_mesh"
+  l
+}
+
+#' Title
+#'
+#' @param x quad_mesh
+#' @param ... args passed to rgl plot
+#'
+#' @return qmesh
+#' @export
+#'
+#' @examples
+#' example(rangl.RasterLayer)
+#' plot(w)
+plot.quad_mesh <- function(x, ...) {
+  ## etc blah
+  ob <- mkq_3d()
+  ob$vb <- t(cbind(as.matrix(x$v[, c("x_", "y_", "z_")]), 1))
+  ob$ib <- matrix(x$qXv$vertex_, nrow = 4)
+  invisible(rgl::shade3d(ob, col = trimesh_cols(nrow(x$qd))[ob$ib], ...))
+}
+
+mkq_3d <- function() {
+  structure(list(vb = NULL, ib = NULL, primitivetype = "quad",
+                 material = list(), normals = NULL, texcoords = NULL), .Names = c("vb",
+                                                                                  "ib", "primitivetype", "material", "normals", "texcoords"), class = c("mesh3d",
+                                                                                                                                                        "shape3d"))
+
+}
+p_4 <- function(xp, nc) {
+  (xp + c(0, 0, rep(nc, 2)))[c(1, 2, 4, 3)]
+}
+#' @importFrom utils tail head
+prs0 <- function(x) {
+  cbind(head(x, -1), tail(x, -1))
+}
+edges0 <- function(x) {
+  #eps <- sqrt(.Machine$double.eps)
+  #as.matrix(expand.grid(seq(xmin(x), xmax(x) -eps, length = ncol(x) + 1),
+  #                      seq(ymax(x), ymin(x) + eps, length = nrow(x) + 1)
+  #))
+  as.matrix(expand.grid(seq(xmin(x), xmax(x), length = ncol(x) + 1),
+                        seq(ymax(x), ymin(x), length = nrow(x) + 1)
+  ))
+}

R/rangl-trip.r

@@ -31,7 +31,7 @@ rangl.trip <-
     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_")
+    outlist$meta <- tibble::tibble(proj = pr4, x = "x_", y = "y_", ctime = format(Sys.time(), tz = "UTC"))
     class(outlist) <- "linemesh"
     outlist
   }

README.Rmd

@@ -36,7 +36,12 @@ Plot methods take those tables and generate the "indexed array" structures neede
 
 The core work for translating "Spatial" classes is done by the unspecialized 'spbabel::map_table' function.
 
-This is likely to be replaced by a 'primitives()' function that takes any lines or polygons data and returns just the linked edges. 
+This is likely to be replaced by a 'primitives()' function that takes any lines or polygons data and returns just the linked edges. Crucially, polygons and lines are described by the same 1D primitives, and this is easy to do. Harder is to generate 2D primitives and for that we rely on [Jonathan Richard Shewchuk's Triangle](https://www.cs.cmu.edu/~quake/triangle.html). 
+
+Triangulation is with `RTriangle` package using "constrained mostly-Delaunay Triangulation" from the Triangle library, but could alternatively use `rgl` with its ear clipping algorithm. 
+
+(With RTriangle we can set a max area for the triangles, so it can wrap around curves like globes and hills.)
+
 ## Installation
 
 This package is in active development and will see a number of breaking changes before release. 
@@ -120,7 +125,7 @@ rgl.snapshot("readme-figure/README-sids-globe.png"); rgl.clear()
 
 ## Holes are trivially supported. 
 
-It's trivial to have "holes", because there are no "holes" because we have a true surface, composed of 2D primitives (triangles). 
+It's trivial to have "holes", because there are no holes, because we have a true surface, composed of 2D primitives (triangles). 
 
 ```{r,eval=TRUE}
 library(spbabel)
@@ -182,7 +187,7 @@ rgl::rgl.snapshot("readme-figure/README-Platonic.png"); rgl.clear()
 
 ![Platonic](readme-figure/README-Platonic.png?raw=true "Platonic")
 
-
+To complete the support for these rgl objects we need quads, and to allows a mix of quads and triangles in one data set (that's what `extrude3d` uses).
 
 
 ## Points
@@ -211,15 +216,59 @@ rgl::rgl.snapshot("readme-figure/README-MultiPoints.png"); rgl::rgl.clear()
 
 ## Trips
 
-Trust me. 
+The soon to be released update to trip includes a 'walrus818' data set courtesy of Anthony Fischbach.  Zoom around and see if you can find them. 
 
-```{r,eval=FALSE}
+```{r,eval=TRUE}
 library(trip)
-example(trip)
-plot(tr)
+library(rangl)
+data(walrus818)
+
+library(graticule)
+prj <-"+proj=laea +lon_0=0 +lat_0=90 +ellps=WGS84"
+gr <- graticule(lats = seq(40, 85, by = 5), ylim = c(35, 89.5), proj = prj)
+library(maptools)
+data(wrld_simpl)
+
+w <- spTransform(subset(wrld_simpl, coordinates(wrld_simpl)[,2] > -70),  prj)
+library(graticule)
+walrus <- spTransform(walrus818, prj)
+
+gr$color_ <- "black"
+rgl.close()
+rgl::par3d(windowRect = c(100, 100, 912 + 100, 912 +100))
+plot(rangl(gr))
+w$color_ <- sample(viridis::inferno(nrow(w)))
+plot(rangl(w), specular = "black")
+plot(rangl(walrus))
+um <- structure(c(0.934230506420135, 0.343760699033737, 0.0950899347662926, 
+0, -0.302941381931305, 0.905495941638947, -0.297159105539322, 
+0, -0.188255190849304, 0.24880850315094, 0.950081348419189, 0, 
+0, 0, 0, 1), .Dim = c(4L, 4L))
+par3d(userMatrix = um)
+rgl::rgl.snapshot("readme-figure/README-walrus.png"); rgl.clear()
+
+```
+
+![Walrus](readme-figure/README-walrus.png?raw=true "Walrus")
+
+
+## Rasters
+
+Single layer rasters are supported. 
+
+```{r}
+wrld <- rasterize(wrld_simpl, raster())
+plot(rangl(wrld/10))
+rgl::rgl.snapshot("readme-figure/README-raster.png"); rgl.clear()
+
+plot(globe(rangl(wrld*100000)), specular = "black")
+rgl::rgl.snapshot("readme-figure/README-rasterglobe.png"); rgl.clear()
+
 ```
 
+![Raster](readme-figure/README-raster.png?raw=true "Raster")
 
+![RasterGlobe](readme-figure/README-rasterglobe.png?raw=true "RasterGlobe")
 
 ## Open topics