The `silicate`

package provides a flexible basis for
working with hierarchical data structures. This includes a clean
separation of topology and geometry, allows *naming* of component
entities, and supports *intermediate*forms. Silicate is a
response to a fragmented landscape where many workarounds and
re-implementations of similar patterns are often repeated.

Silicate is composed of *normal-form models* and worker verbs
for extracting underlying entities. By normal-form, we mean the tidy, 3NF,
de-duplicated, *removal of redundancy* sense.

The entities that we identify correspond to the names of the worker verbs. Whether a given model has all or any of these entities explicitly depends on the circumstances, but implicitly these are usually always present and depend only on interpretation.

`sc_object`

- highest level properties, the “features”`sc_coord`

- all instances of coordinates, labelled by vertex if the source model includes them`sc_vertex`

- only unique coordinates (in some geometric space)`sc_path`

- individual paths, sequential traces`sc_edge`

- unique binary relations, unordered segments`sc_segment`

- all instances of edges`sc_arc`

- unique topological paths, arcs either meet two other arcs at a node, or include no nodes`sc_node`

- unique nodes

Finally there is an important extra verb, the `unjoin`

- a
function to *un join* a table, the opposite of the database join
which is used as a key step when building these models, used to remove
duplication at various levels. It’s the primary mechanism for defining
and building-in *topology*, which is (precisely) the
relationships between entities in a model.

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 = "")
<- function(x) head(do.call(rbind, unlist(lapply(x, function(a) rbind(a, NA)), recursive = F)), -1)
rbind_na <- sc_colours(nrow(minimal_mesh))
cols <- lapply(seq_along(minimal_mesh$geom), function(y) polypath(rbind_na(minimal_mesh$geom[[y]]), col = cols[y] )) junk
```

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.

```
<- SC(minimal_mesh)
x names(x)
```

```
## [1] "object" "object_link_edge" "edge" "vertex"
## [5] "meta"
```

`$edge x`

```
## # 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 OmBgld
```

`print(x)`

```
## class : SC
## type : Primitive
## vertices : 14 (2-space)
## primitives : 15 (1-space)
## crs : NA
```

```
plot(x)
text(x$vertex[c("x_", "y_")], label = x$vertex$vertex_)
```

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.

```
<- PATH(minimal_mesh)
x 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.

`::reduce(x[c("object", "path", "path_link_vertex", "vertex")], dplyr::inner_join) purrr`

```
## 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!

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(minimal_mesh)
arc <- sc_node(arc)
nodes
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 …

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(minimal_mesh)
tri 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.806
```

The 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))`