knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.path = "README-" )
suppressMessages({library(gris) library(maptools) library(cgalgris)})
The gris package provides a relational geometry/topology model for spatial data in R. This is inspired by data models used in the commercial products Manifold GIS and Eonfusion. The main aspirations are
Related work
Here we load a commonly used layer of country polygons and convert to the gris format of tables with vertices/branches/objects.
library(gris) library(maptools) data(wrld_simpl) dat <- gris(wrld_simpl) #dat1 <- sbs(dat, filter(dat$o, NAME %in% c("Australia", "Indonesia", "Papua New Guinea"))) dat1 <- dat[which(dat$o$NAME %in% c("Australia", "Indonesia", "Papua New Guinea")), ] library(dplyr) #pl(dat1$v, asp = 1.3) plot(dat1)
The function bld
builds a gris object from a SpatialPolygonsDataFrame, sbs
provides simple subsetting based on dplyr, and pl
plots the individual branches of each object using polypath with a unique colour for objects.
Triangulate with CGAL via cgalgris. The function tri_xy
performs an exact Delaunay triangulation on all vertices, returning a triplet-index for each triangle (zero-based in CGAL). Then process all triangles to detect which fall within the input polygons and plot them separately. (This test is inexact, since some polygon boundaries will not conform to the overall triangulation but CGAL does provide this ability).
library(cgalgris) ## Delaunay triangulation (unconstrained) dat1$vi <- tri_xy(dat1$v$x_, dat1$v$y_) + 1 ## plus 1 for R indexing ## centroids of triangles centr <- data_frame(x = dat1$v$x_[dat1$vi], y = dat1$v$y_[dat1$vi], t = rep(seq(length(dat1$vi)/3), each = 3)) %>% group_by(t) %>% summarize(x = mean(x), y = mean(y)) %>% select(x, y, t) x1 <- dat1$v %>% inner_join(dat1$bXv) %>% mutate(mg = branch_) %>% group_by(mg) %>% do(rbind(., NA_real_)) inside <- which(point.in.polygon(centr$x, centr$y, x1$x_, x1$y_) == 1) ## plot all triangles plot(dat1$v$x_, dat1$v$y_, type = "n", asp = 1.3) apply(matrix(dat1$vi, ncol = 3, byrow = TRUE), 1, function(x) polypath(cbind(dat1$v$x_[x], dat1$v$y_[x]), col = "#FF000080")) ## overplot only those that are internal to the polygons (not exact since triangulation is unconstrained) apply(matrix(dat1$vi, ncol = 3, byrow = TRUE)[inside, ], 1, function(x) polypath(cbind(dat1$v$x_[x], dat1$v$y_[x]), col = "#0066FF99", border = "grey")) points(centr[, c("x", "y")], cex = 0.2)
Future versions will leverage more of CGAL for this functionality.
type
argument to pl
for interpretation of branches when plottingAdd the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.