MODELS
Our capitalized model forms are all built from workflows that use
these verbs. Here we illustrate with a very simple simple
features model minimal_mesh, two polygons that share a
single edge.
library(silicate)
plot(matrix(attr(minimal_mesh$geom, "bbox"), 2, byrow = TRUE), type = "n", xlab = "", ylab = "")
rbind_na <- function(x) head(do.call(rbind, unlist(lapply(x, function(a) rbind(a, NA)), recursive = F)), -1)
cols <- sc_colours(nrow(minimal_mesh))
junk <- lapply(seq_along(minimal_mesh$geom), function(y) polypath(rbind_na(minimal_mesh$geom[[y]]), col = cols[y] ))
SC
This is the universal model, the one model that can be used to represent anything at all. The key entity in SC is the edge, a single binary relationship identifying a line segment joining two vertices. Note how each edge and each vertex has a unique label. Since each primitive has exactly the same structure this model presents a number of optimizations, and all structures can be expressed in these terms. We can even abandon any record of sequential paths in these models because they can be re-derived by tracing through the labels at a later time.
x <- SC(minimal_mesh)
names(x)## [1] "object" "object_link_edge" "edge" "vertex"
## [5] "meta"x$edge## # A tibble: 15 × 4
## .vx0 .vx1 path_ edge_
## <chr> <chr> <int> <chr>
## 1 o9jchu phJiJJ 1 kpBBjm
## 2 phJiJJ TINXMR 1 gYd7QP
## 3 TINXMR cpsLif 1 wA9gmq
## 4 gtuG1S cpsLif 1 UYEaEX
## 5 gtuG1S aSzy6I 1 HyW4tM
## 6 u3yvlS aSzy6I 1 Dbtzw1
## 7 o9jchu u3yvlS 1 hb6yYN
## 8 SaLOSt vzDobV 2 eb7dvh
## 9 vzDobV 8yO6N3 2 Bw4uGG
## 10 x1bPEZ 8yO6N3 2 TpEhHN
## 11 fgHLvc x1bPEZ 2 rjfqR8
## 12 SaLOSt fgHLvc 2 evVPcA
## 13 aSzy6I zxh9ud 3 pFK3IG
## 14 zxh9ud vTfrp9 3 9G3MzP
## 15 u3yvlS vTfrp9 3 OmBgldprint(x)## class : SC
## type : Primitive
## vertices : 14 (2-space)
## primitives : 15 (1-space)
## crs : NAplot(x)
text(x$vertex[c("x_", "y_")], label = x$vertex$vertex_)
PATH
At the highest level silicate provides a normalized form of a complicated structure containing paths. This data structure contains three paths, each of which is a sequentially joined set of coordinates.
The key aspect here is that each component entity is formally named, with a unique ID that we persist for subsequent usage. This persistence is required as it records relationships between existing entitites implicitly rather than baking data into structure that are both explicit and also discard a lot of information about relationships.
Even though the path scheme requires that two vertices are visited separately by two polygon rings, we don’t expand out the vertices but simply record a reference to it in the path.
x <- PATH(minimal_mesh)
names(x)## [1] "object" "path" "path_link_vertex" "vertex"
## [5] "meta"At a basic level, the tables in this PATH are
essentially the same as the kinds of structures we normally use, but we
need to de-normalize them to see this exactly.
purrr::reduce(x[c("object", "path", "path_link_vertex", "vertex")], dplyr::inner_join)## Joining, by = "object_"
## Joining, by = "path_"
## Joining, by = "vertex_"## # A tibble: 19 × 10
## a object_ ncol type subobject path_ ncoor…¹ vertex_ x_ y_
## <int> <chr> <int> <chr> <int> <chr> <int> <chr> <dbl> <dbl>
## 1 1 zGYdgL 2 MULTIPOLYGON 1 UBt593 8 nyKRKa 0 0
## 2 1 zGYdgL 2 MULTIPOLYGON 1 UBt593 8 SMTDgZ 0 1
## 3 1 zGYdgL 2 MULTIPOLYGON 1 UBt593 8 cSx5S0 0.75 1
## 4 1 zGYdgL 2 MULTIPOLYGON 1 UBt593 8 SEqPLv 1 0.8
## 5 1 zGYdgL 2 MULTIPOLYGON 1 UBt593 8 KNG9qn 0.5 0.7
## 6 1 zGYdgL 2 MULTIPOLYGON 1 UBt593 8 VxB7KB 0.8 0.6
## 7 1 zGYdgL 2 MULTIPOLYGON 1 UBt593 8 8J5Qau 0.69 0
## 8 1 zGYdgL 2 MULTIPOLYGON 1 UBt593 8 nyKRKa 0 0
## 9 1 zGYdgL 2 MULTIPOLYGON 1 oPl1uN 6 jAqzPc 0.2 0.2
## 10 1 zGYdgL 2 MULTIPOLYGON 1 oPl1uN 6 vMI8t4 0.5 0.2
## 11 1 zGYdgL 2 MULTIPOLYGON 1 oPl1uN 6 cbgT0w 0.5 0.4
## 12 1 zGYdgL 2 MULTIPOLYGON 1 oPl1uN 6 C6Fxpl 0.3 0.6
## 13 1 zGYdgL 2 MULTIPOLYGON 1 oPl1uN 6 9AIulm 0.2 0.4
## 14 1 zGYdgL 2 MULTIPOLYGON 1 oPl1uN 6 jAqzPc 0.2 0.2
## 15 2 VdjyDC 2 MULTIPOLYGON 1 dAu4HA 5 8J5Qau 0.69 0
## 16 2 VdjyDC 2 MULTIPOLYGON 1 dAu4HA 5 VxB7KB 0.8 0.6
## 17 2 VdjyDC 2 MULTIPOLYGON 1 dAu4HA 5 o3ov7m 1.1 0.63
## 18 2 VdjyDC 2 MULTIPOLYGON 1 dAu4HA 5 rEZvqE 1.23 0.3
## 19 2 VdjyDC 2 MULTIPOLYGON 1 dAu4HA 5 8J5Qau 0.69 0
## # … with abbreviated variable name ¹ncoords_Paths are easily re-expressed as a set of edges, treating each pair of coordinates as a new kind of entity. Some edges are shared, in that two objects might have a particular edge in their own path. If the shared neighbours are two polygons, then it’s likely that each polygon encounters that edge in a different direction, but this is not a stable pattern so it’s best to not assume it to be the case. It’s important here because we need to differentiate between an instance of an edge, which is a particular line segment that is part of one particular polygon. If that polygon segment is a shared edge with a neighbour polygon, then we distinguish the instances as a particular segment of a unique edge. This allows us the opportunity to treat the edge as completely abstract, and also to decide or record what it’s orientation is, which records which vertex is the first one, and which is the second. Not all applications need to care about this distinction, though.
In this way we have an analogy for edges and segments compared to vertices and coordinates, and I think this terminology makes sense!
ARC
There’s only one more entity we need to describe an alternative view of this polygon data, the actual component paths that describe these structures topologically. These are paths that occur between any two nodes, or that have no nodes at all. This will occur for any isolated “island”, or for any hole within a polygon. These arcs are exactly analogous to a LINESTRING, but are used within a context where the nodes are important information about the topology.
Nodes are any coordinate where three or more edges meet each other, so they are inherently about relationships between neighbouring shapes. This is exactly the same model as was used in the olden days, called arc-node topology, a one-dimensional topology that pre-dates the more modern use of closed planar paths to describe polygons.
library(dplyr)
arc <- ARC(minimal_mesh)
nodes <- sc_node(arc)
plot(arc)
inner_join(nodes, arc$vertex) %>% dplyr::select(x_, y_) %>% points(pch = "N")## Joining, by = "vertex_"
We haven’t yet done any coordinate precision checking, so there will be mis-identified arcs where coordinates are intended to be but are not exactly the same. (spdep glosses over this in a way we can’t control)
Example from wrld_simpl …
TRI
Triangulations. Here we strictly want constrained triangulations, where there is no divergence from the planar coverage of input shape and its triangulated decomposition. This is not easy to do, and there are two main methods. One is high-quality mostly-Delaunay triangulations where the Delaunay condition is relaxed in order to always align to input boundaries. Holes in a near-Delaunay triangulation are either removed by flood-fill from a seed-vertex, or by post-filtering the triangles whose centroid does not intersect the input boundaries. The second method is ear-clippping, a relatively inexpensive algorithm that produces low-quality triangle meshes.
Currently silicate can only provide the second form, the ear-clipping
method provided by decido.
tri <- TRI(minimal_mesh)
plot(tri)
Performance is pretty good for this method of triangulation, despite also being constrained to the input shapes.
#system.time(sf::st_triangulate(inlandwaters))
# user system elapsed
# 4.699 0.125 4.823
#system.time(sfdct::ct_triangulate(inlandwaters))
# user system elapsed
# 15.476 0.225 15.460
system.time(tri <- TRI(inlandwaters))## user system elapsed
## 5.744 0.051 5.806The triangles are low-quality (but triangles are very easy to subdivide).
plot(tri)
plot(NA, xlim = c(625000, 1060000), ylim = c(-1350000, -550000))
plot(tri, add = TRUE)
The CRAN package
sfdct will use the first method (package RTriangle
wraps the Triangle library) for simple features.
#plot(sfdct::ct_triangulate(minimal_mesh))