07-topology-and-geocoding.Rmd

```{r echo=FALSE, warning=FALSE}
yml_content <- yaml::read_yaml("chapterauthors.yml")
author <- yml_content[["topologyGeocoding"]][["author"]]
```

# Topology and Geocoding {#topology-and-geocoding}

Written by
```{r results='asis', echo=FALSE}
cat(author)
```

``` {r 7-load-packages-and-setup, echo = FALSE, include = FALSE}
baseurl <- "https://ubc-geomatics-textbook.github.io/geomatics-textbook/"
canadalicense <- "https://open.canada.ca/en/open-government-licence-canada/"
covlicense <- "https://opendata.vancouver.ca/pages/licence/"
library(leaflet)
library(dplyr)
library(deldir)
library(tripack)
library(animation)
library(alphashape3d)
library(alphahull)
library(lidR)
options(rgl.useNULL=TRUE)
library(rgl)
library(mgcv)
library(sf)
knitr::knit_hooks$set(webgl = hook_webgl)
```

Frequently, we need spatial data to behave and relate in specific and predictable ways. Many types of analyses may expect spatial data to be represented and interact in a standard form. In this chapter, we will extend our knowledge of data models using topology, which unlocks many advanced spatial analyses. We will look at a specific example of an analysis that requires topology, geocoding, which will be a convenient segue into network analysis discussed in the following chapter.

:::: {.box-content .learning-objectives-content}

::: {.box-title .learning-objectives-top}
## Learning Objectives {-}
::: 

1. Understand the role of topology in governing data behaviour and data organization
2. Recognize some examples and uses of 2D and 3D topologies
3. Understand the role of bounding a set of points from triangulation and convex hulls
4. Synthesize the process of geocoding
5. Practice geocoding addresses and reverse geocoding addresses to other coordinate systems

::::

## Key Terms {-}

Vertex, Node, Pseudonode, Dangle, Planar Topology, Non-Planar Topology, Geocoding, Adjacency, Overlap, Connect, Inside, Reverse Geocoding, Singlepart, Multipart, Holes, Delaunay Triangulation, Thiessen Polygons, Voronoi Diagram, Centroid, Convex Hull, Convex Alpha Hull, Multipatch

## Topology

**Topology** describes the relationships of spatial data. This is a very broad definition that encompasses the wide range of possible arrangements of spatial data in practice. If we drill down into this concept, topology is really what allows us undertake specific types of analysis that requires or expects spatial data to behave in a certain way. If you think about the feature geometries that we have at our disposal, then there are no fewer than nine combinations of how these geometries can interact as illustrated in Figure \@ref(fig:7-topology-grid) below.

```{r 7-topology-grid, echo=FALSE, out.width = '75%', fig.cap=fig_cap}
fig_cap <- paste0("Grid showing all the combinations of how point, line and polygon geometries can interact. Pickell, CC-BY-SA-4.0.")
knitr::include_graphics("images/07-topology-grid.png")
```

It is important to recognize that there may be cases where we "expect" that a given combination of features will conform to a specific interaction. For example, the provinces and territories of Canada are typically represented as polygons that share adjacent boundaries. That is, the adjacent boundaries shared between any two provinces or territories cannot logically overlap as this representation (model) would contravene the legal definitions of the provinces and territories. In another case, a human technician may erroneously digitize a road that crosses another road without indicating that the two roads share an intersection, which could have consequences for how traffic flow can be modeled between the two roads (i.e., intersection with traffic light versus overpass). These are both examples of situations where topology is needed. Topology applies logic to define how features are expected to relate to other features in order to conform to knowledge systems like legal definitions of land and traffic flow. In short, topology ensures data integrity for other types of analysis.

## Planar vs. Non-Planar Topology

In the context of topology, **planar** refers to the concept that all vertices of feature vector geometry are mapped onto the same plane. So in a planar world-view, all lines and polygons share coincident vertices. For example, if two polygons overlap, then the overlapping area forms a new polygon with a boundary of vertices defined by the union of the two other polygons. Also, if two lines overlap, then the two lines are divided into four new segments and a new vertex is formed at the intersection. In other words, planar topology does not allow polygons or lines to lay "underneath" or "on top" of another line or polygon and feature geometries must always be distinct.

On the other hand, **non-planar** topology is the concept that vertices of feature vector geometry can be mapped to different planes. It is important to emphasize here that when we are talking about planes that we are not referring to projected coordinate systems. It is generally assumed that any two spatial data layers containing feature geometries are interacting within the same projected coordinate system. Non-planar topology allows for other knowledge systems to be represented in spatial data. The case where a pipeline runs underneath a river or a territory that was traditionally used by several Indigenous peoples (Figure \@ref(fig:7-native-land-leaflet)) are examples of valid non-planar topology.

``` {r, 7-native-land-data-load, echo=F, message=F, warning=F, cache=T}
native_land <- st_as_sfc("Data/07/07-native_land_vancouver_island.geojson", crs = st_crs("Data/07/07-native_land_vancouver_island.geojson"), GeoJSON = TRUE)
labels <- native_land$Name
```

``` {r 7-native-land-leaflet-setup, cached=T, echo = FALSE, include = FALSE}
getTilesURL <- "https://tile.thunderforest.com/landscape/{z}/{x}/{y}.png?apikey=06d706f23fc4461db5db742f930174fc"
m <- leaflet(native_land) %>%
  addTiles(getTilesURL) %>%
  addScaleBar(position = c("topleft")) %>%
  setView(lat = 49.25, lng = -125, zoom = 7) %>%
  addPolygons(data=native_land,
              fillColor = native_land$color,
              weight = 0.6,
              opacity = 0.5,
              color = "#ffffff",
              dashArray = NULL,
              fillOpacity = 0.25,
              highlight = highlightOptions(
                weight = 1.8, 
                color = "#7d7d7d",
                dashArray = NULL,
                fillOpacity = 0.7,
                bringToFront = TRUE),
              label = labels,
              labelOptions = labelOptions(
                style = list("font-weight" = "normal", padding = "3px 8px"),
                textsize = "12px", direction = "auto"))
```

(ref:7-native-land-leaflet-caption1) Non-planar topology of 36 indigenous territories overlapping Vancouver Island, British Columbia. Data from [@native_land_native_nodate], CC0. Animated figure can be viewed in the web browser version of the textbook: https://www.opengeomatics.ca/topology-and-geocoding.html.

(ref:7-native-land-leaflet-caption2) Non-planar topology of 36 indigenous territories overlapping Vancouver Island, British Columbia. Data from [@native_land_native_nodate], CC0.

```{r 7-native-land-leaflet, fig.cap = fig_cap, out.width= "100%", echo = FALSE}
if (knitr:::is_latex_output()) { fig_cap ="(ref:7-native-land-leaflet-caption1)"
knitr::include_graphics("images/07-native-land-leaflet.png")
} else { fig_cap = "(ref:7-native-land-leaflet-caption2)"
m
}
```

## Implementing Planar Topology

Implementing planar topology involves defining specific rules for how features should relate to one another given some analytical context. This process also requires that the spatial data are housed a relational database or data model that supports topology. In other words, topology is enforced only by data models that support topological rules. When a topological rule is violated, the relational database identifies the contravening features and displays them on the map and in the attribute table. Then, it is up to an analyst to decide how the error should be corrected. For example, some errors like intersecting lines can automatically be split at the intersection while overlapping polygons might need to be manually edited to reflect the correct adjacency. Thus, the process of applying topology is first to work within a data model that supports topology, then choose the topological rules that reinforce a particular knowledge system, and finally to inspect and decide how to deal with any contraventions. Since planar topology is only supported by certain data models, and some data models are proprietary to certain software, the exact topological rules that can be implemented in a GIS are mostly dependent on the software that you are using. Instead of examining a specific GIS software package, we will discuss the "fundamental" planar topological relationships that are common across nearly all implementations of topology. (If you want to know more about how topology is implemented within specific data models, skip ahead to the "Data models supporting planar topology" section.)

So far, we have seen that there are six ways to combine feature geometries (points, lines, and polygons). We can extend this understanding to include at least six different ways that they can relate to one another: adjacent; overlap; intersect; connect; cover; and inside. Some of these relationships can be modeled *between* two different spatial layers (e.g., two point layers) or *within* a single spatial layer. In the following sections, we will look at different planar topological rules that apply both between and within feature geometries.

## Adjacency and Overlap

There are times when we need to ensure that two polygons are **adjacent** to one another by sharing a common edge. If two polygons are not adjacent to one another, then a gap, known as a **sliver**, exists between them or they must **overlap**. Consider the case where we are mapping land covers. If we have a formal scheme that describes all possible land covers, then we expect that a map of land covers will have perfect adjacency between all polygons so that there are no areas that are not mapped (i.e., slivers) and that no area has multiple, overlapping land covers. Since lines are also 2-dimensional, lines can overlap other lines. Depending on the context, a topological rule may be needed to promote or prevent this relationship. For example, if you are modeling bus routes, then one road might support several different routes.

``` {r echo=FALSE, warning=FALSE, message=FALSE}
countries <- st_read("data/02/world-countries.shp", quiet=T)
uscanada <- countries[countries$CNTRY_NAME == "Canada" | countries$CNTRY_NAME == "United States",]
uscanada_sf <- st_as_sf(uscanada)
uscanada_simp <- st_simplify(uscanada_sf, preserveTopology = FALSE, dTolerance = 1000)
```

``` {r echo=FALSE, warning=FALSE, message=FALSE}
pal <- colorFactor(palette = c("#a50026","#313695"), levels = c("Canada","United States"))
m = leaflet(uscanada_simp) %>%
  addTiles() %>%
  setView(lng = -94.0637905778631, lat = 48.64363272747508, zoom = 9) %>%
  addScaleBar(position = c("topleft")) %>%
  addPolygons(stroke = T, color=~pal(CNTRY_NAME), weight=1, smoothFactor = 0.3, fillOpacity = 0.8)
```

``` {r echo=FALSE, out.width="100%", fig.cap=fig_cap}
if (knitr:::is_latex_output()) { fig_cap <- paste0("Simplified boundary of Canada (red) and the United States (blue), centred at the Rainy River boundary between Minnesota and Ontario. Without preserving topology, the output results in illogical overlaps in some places and slivers in other places. Pickell, CC-BY-SA-4.0. Data from Natural Resources Canada and licensed under the Open Government Licence - Canada. Animated figure can be viewed in the web browser version of the textbook: https://www.opengeomatics.ca/topology-and-geocoding.html")
  knitr::include_graphics("images/07-Canada-US-border-simplified-leaflet.PNG")
} else { fig_cap <- paste0("Simplified boundary of Canada (red) and the United States (blue), centred at the Rainy River boundary between Minnesota and Ontario. Without preserving topology, the output results in illogical overlaps in some places and slivers in other places. Pickell, CC-BY-SA-4.0. Data from Natural Resources Canada and licensed under the Open Government Licence - Canada.")
m
}
```

Some examples of adjacency and overlap topological rules:

- Polygons within the same layer must not have gaps
- Polygons within the same layer must not overlap
- Polygons must not overlap other polygons
- Lines must not overlap other lines
- Lines must not self-overlap

## Intersect and Connect

As we have seen from [Chapter 3](https://ubc-geomatics-textbook.github.io/geomatics-textbook/types-of-data.html), lines are often used to represent phenomena that flow, so intersection and connection are important concepts for these representations. Important to understanding how connection and intersection work in planar topology, we need to understand that lines are comprised of a set of vertices and nodes. A **node** is simply the terminating vertex in a set of vertices for a line. For example, suppose the line segment $A$ has a set of vertices, $[[1,0],[1,3],[1,5]]$. Then the nodes for $A$ are $[1,0]$ and $[1,5]$ (Figure \@ref(fig:7-node)). Since nodes define the end points of a line segment, they are key to enforcing connection rules. We will look at network analysis in more detail in the next chapter. For now, let us consider two different networks that can help us conceptualize some fundamental line topology using nodes.

```{r, 7-node, fig.align = 'center', out.width = '75%', fig.cap=fig_cap, cached=T, echo = FALSE}
fig_cap <- paste0("Lines are always comprised of two nodes. Line A shown here has nodes at [1,0] and [1,5]. Pickell, CC-BY-SA-4.0.")
Ax <- c(1,1,1)
Ay <- c(0,3,5)
plot(1,1,ylim=c(-1,6),xlim=c(-1,6),type="n",xlab="",ylab="",xaxt="n",yaxt="n")
axis(1,at=-1:6)
axis(2,at=-1:6)
abline(h=-1:6,v=-1:6,col="grey")
lines(Ax,Ay,type="o",lwd=2,pch=19,col="#7b3294")
labs <- c("(1,0)","(1,3)","(1,5)")
text(Ax,Ay-0.2,labels=labs, cex=0.9, font=2, pos=4,col="#7b3294")
text(0.75,2.5,"A", cex=1, font=2,col="#7b3294")
```

A network of streams and rivers is based on the hydrological knowledge system that explains how water moves over a terrain surface. In both theory and practice, we know that water flows from higher elevations to lower elevations with limited exceptions. Thus, we expect that streams will connect with other streams and continue to flow towards some outlet such as an ocean. **Connection** refers to the fact that the endpoint node of one stream will fall somewhere on another stream segment. Where two line segments come together, it is possible for one segment $A$ to "undershoot" the other segment $B$, resulting in the end node of segment $A$ appropriately named a **dangle** (Figure \@ref(fig:7-dangle)) and a loss of connection.

```{r, 7-dangle, fig.align = 'center', out.width = '75%', fig.cap=fig_cap, cached=T, echo = FALSE}
fig_cap <- paste("A dangle forms when a line (B) does not connect to another line (A). Pickell, CC-BY-SA-4.0.")
Ax <- c(1,1,1)
Ay <- c(0,3,5)
Bx <- c(2,4)
By <- c(2,2)
plot(1,1,ylim=c(-1,6),xlim=c(-1,6),type="n",xlab="",ylab="",xaxt="n",yaxt="n")
axis(1,at=-1:6)
axis(2,at=-1:6)
abline(h=-1:6,v=-1:6,col="grey")
lines(Ax,Ay,type="o",lwd=2,pch=19,col="#7b3294")
lines(Bx,By,type="o",lwd=2,pch=19,col="#008837")
labs <- c("(1,0)","(1,3)","(1,5)")
text(Ax,Ay-0.2,labels=labs, cex=0.9, font=2, pos=4,col="#7b3294")
labs <- c("(2,2)","(4,2)")
text(Bx,By-0.2,labels=labs, cex=0.9, font=2, pos=4,col="#008837")
text(0.75,2.5,"A", cex=1, font=2,col="#7b3294")
text(3.5,2.5,"B", cex=1, font=2,col="#008837")
```

Dangles are the opposite case to **intersections**, which occur when two line segments cross each other. With planar topology, intersections must be modeled with a shared node representing the intersection location. For example, suppose line segment $B$ has a set of vertices, $[[0,1],[2,1],[4,1]]$. If line segments $A$ (defined above) and $B$ are mapped together with non-planar topology, then they will intersect at $[1,2]$, which is not a vertex represented in either segment (Figure \@ref(fig:7-intersection-A-B)). 

```{r, 7-intersection-A-B, fig.align = 'center', out.width = '75%', fig.cap=fig_cap, cached=T, echo = FALSE}
fig_cap <- paste0("Line A mapped with Line B in non-planar topology. Pickell, CC-BY-SA-4.0.")
Ax <- c(1,1,1)
Ay <- c(0,3,5)
Bx <- c(0,2,4)
By <- c(2,2,2)
plot(1,1,ylim=c(-1,6),xlim=c(-1,6),type="n",xlab="",ylab="",xaxt="n",yaxt="n")
axis(1,at=-1:6)
axis(2,at=-1:6)
abline(h=-1:6,v=-1:6,col="grey")
lines(Ax,Ay,type="o",lwd=2,pch=19,col="#7b3294")
lines(Bx,By,type="o",lwd=2,pch=19,col="#008837")
labs <- c("(1,0)","(1,3)","(1,5)")
text(Ax,Ay-0.2,labels=labs, cex=0.9, font=2, pos=4,col="#7b3294")
labs <- c("(0,2)","(2,2)","(4,2)")
text(Bx,By-0.2,labels=labs, cex=0.9, font=2, pos=4,col="#008837")
text(0.75,2.5,"A", cex=1, font=2,col="#7b3294")
text(3.5,2.5,"B", cex=1, font=2,col="#008837")
```

Thus, the intersection of $A$ and $B$ with planar topology would yield four new segments: $C=[[0,2],[1,2]]$, $D=[[1,2],[1,3],[1,5]]$, $E=[[1,2],[2,2],[4,2]]$, and $F=[[1,0],[1,2]]$. Figure \@ref(fig:7-intersection-C-D-E-F) illustrates how all four of these new segments share the same node $[1,2]$ at the intersection of $A$ and $B$.

```{r, 7-intersection-C-D-E-F, fig.align = 'center', out.width = '75%', fig.cap=fig_cap, cached=T, echo = FALSE}
fig_cap <- paste("Line A mapped with Line B in planar topology yields segments C, D, E, and F. All segments share (1,2) as a node. Pickell, CC-BY-SA-4.0.")
Cx <- c(0,1)
Cy <- c(2,2)
Dx <- c(1,1,1)
Dy <- c(2,3,5)
Ex <- c(1,2,4)
Ey <- c(2,2,2)
Fx <- c(1,1)
Fy <- c(0,2)
plot(1,1,ylim=c(-1,6),xlim=c(-1,6),type="n",xlab="",ylab="",xaxt="n",yaxt="n")
axis(1,at=-1:6)
axis(2,at=-1:6)
abline(h=-1:6,v=-1:6,col="grey")
lines(Cx,Cy,type="o",lwd=2,pch=19,col="#e41a1c")
lines(Dx,Dy,type="o",lwd=2,pch=19,col="#377eb8")
lines(Ex,Ey,type="o",lwd=2,pch=19,col="#f781bf")
lines(Fx,Fy,type="o",lwd=2,pch=19,col="#ff7f00")
labs <- c("(0,2)","")
text(Cx,Cy-0.2,labels=labs, cex=0.9, font=2, pos=4,col="#e41a1c")
labs <- c("","(1,3)","(1,5)")
text(Dx,Dy-0.2,labels=labs, cex=0.9, font=2, pos=4,col="#377eb8")
labs <- c("","(2,2)","(4,2)")
text(Ex,Ey-0.2,labels=labs, cex=0.9, font=2, pos=4,col="#f781bf")
labs <- c("(1,0)","")
text(Fx,Fy-0.2,labels=labs, cex=0.9, font=2, pos=4,col="#ff7f00")
text(0.5,2.5,"C", cex=1, font=2,col="#e41a1c")
text(3.5,2.5,"E", cex=1, font=2,col="#f781bf")
text(0.75,0.5,"F", cex=1, font=2,col="#ff7f00")
points(1,2,pch=19)
text(1,2-0.2,"(1,2)",pos=4,font=2,cex=0.9)
```

As well, **pseudonodes** can occur when a node does not actually terminate a line segment at a junction, for example, between two streams or roads. In other words, a pseudonode is a node that is shared by two lines. Figure \@ref(fig:7-pseudonode) illustrates a pseudonode occurring at $[3,5]$.

```{r, 7-pseudonode, fig.align = 'center', out.width = '75%', fig.cap=fig_cap, cached=T, echo = FALSE}
fig_cap <- paste0("Lines A and B share a pseudonode at [3,5], indicated in red. Pickell, CC-BY-SA-4.0.")
Ax <- c(0,1,3)
Ay <- c(0,3,5)
Bx <- c(3,4,3)
By <- c(5,3,1)
plot(1,1,ylim=c(-1,6),xlim=c(-1,6),type="n",xlab="",ylab="",xaxt="n",yaxt="n")
axis(1,at=-1:6)
axis(2,at=-1:6)
abline(h=-1:6,v=-1:6,col="grey")
lines(Ax,Ay,type="o",lwd=2,pch=19,col="#7b3294")
labs <- c("(0,0)","(1,3)","(3,5)")
text(Ax,Ay-0.2,labels=labs, cex=0.9, font=2, pos=4,col="#7b3294")
text(0.5,3.5,"A", cex=1, font=2,col="#7b3294")
lines(Bx,By,type="o",lwd=2,pch=19,col="#008837")
labs <- c("(3,5)","(4,3)","(3,1)")
text(Bx,By-0.2,labels=labs, cex=0.9, font=2, pos=4,col="#008837")
text(3.5,3,"B", cex=1, font=2,col="#008837")
points(3,5,col="red",pch=19)
```

Some examples of intersection and connection topological rules:

- Lines must not intersect other lines
- Lines must intersect other lines
- Lines must not self-intersect
- Lines within a same layer must not self-intersect
- Lines must not have dangles

## Coincident and Disjoint

Point features can be either **coincident** or **disjoint** with other point features. Point features that need to be disjoint may be representing trees, mountain peaks, or any similar type of feature that would be expected to be discrete in geographic space. There are also instances where we might need one set of point features to be coincident with another such as field plots that are centered using a tree or other spatially-discrete feature on the landscape.

Some examples of coincident and disjoint topological rules:

- Points must be disjoint with other points
- Points must by coincident with other points

## Cover

**Cover** refers to planar topology where a feature lays on or within another feature. For example, dams represented as point features must be covered by a line representing a river (Figure \@ref(fig:7-dam-river-watershed)). Similarly, lines representing rivers must be covered by polygons representing watersheds. As well, property parcel polygons must be covered by the municipal or regional tax authority polygon.

```{r 7-dam-river-watershed, echo=FALSE, out.width = '50%', fig.cap=fig_cap}
fig_cap <- paste("Topological relationship between dam (point) covered by a river (line), which is covered by a watershed (polygon). Pickell, CC-BY-SA-4.0.")
knitr::include_graphics("images/07-dam-river-watershed.png")
```

Some examples of cover topological rules:

- Point must be covered by a line
- Point must be covered by a polygon
- Line must be covered by a polygon
- Polygon must be covered by a polygon

### Multipart geometry

Sometimes we need to represent several points, lines or polygons as a collection, which is known as **multipart geometry**. Multipart geometry allow us to represent several disjoint and non-adjacent geometries as a single feature. In this way, we can assign attribute values to the collection of features rather than each geometry individually. The territorial boundary of Canada is a good example of an instance where a multipart geometry can be useful because all of the contiguous land and non-contiguous land (i.e., islands) can be represented and associated with a single feature in the attribute table. However, if the distinction of features is important, such as identifying the names of islands in the Haida Gwaii archipelago, then a singlepart geometry should be used (Figure \@ref(fig:7-haida-gwaii-leaflet)).

``` {r 7-haida-leaflet-setup, echo=FALSE, warning=FALSE, message=FALSE}
haida <- st_as_sfc("Data/07/07-haida-gwaii-islands-names.geojson", crs = st_crs("Data/07/07-haida-gwaii-islands-names.geojson"), GeoJSON = TRUE)
labels2 <- haida$Island_Name
```

``` {r 7-haida-leaflet-setup2, echo=FALSE, warning=FALSE, message=FALSE, cached=T, echo = FALSE}
getTilesURL2 <- "https://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}{r}.png"
m2 = leaflet(haida) %>%
  addTiles(getTilesURL2) %>%
  addScaleBar(position = c("topleft")) %>%
  setView(lat = 53, lng = -132, zoom = 7) %>%
  addPolygons(data=haida,
              fillColor = "#993300",
              weight = 0.6,
              opacity = 1,
              color = "#000000",
              dashArray = NULL,
              fillOpacity = 1,
              highlight = highlightOptions(
                weight = 1.8, 
                color = "#7d7d7d",
                dashArray = NULL,
                fillOpacity = 0.7,
                bringToFront = TRUE),
              label = labels2,
              labelOptions = labelOptions(
                style = list("font-weight" = "normal", padding = "3px 8px"),
                textsize = "12px", direction = "auto"))
```

(ref:7-haida-gwaii-leaflet-caption) Singlepart geometry of the Haida Gwaii archipelago off the west coast of British Columbia, Canada. Hover over the islands to see the names. Polygon data from [@statistics_canada_statistics_nodate] and island placenames from [@natural_resources_canada_natural_nodate]. Open Government License - Canada: https://open.canada.ca/en/open-government-licence-canada/. Interactive web map can be viewed in the web browser version of the textbook: https://www.opengeomatics.ca/topology.html#fig:7-haida-gwaii-leaflet.

```{r 7-haida-gwaii-leaflet, fig.cap = fig_cap, out.width= "50%", echo = FALSE}
if (knitr:::is_latex_output()) { 
  fig_cap <- "(ref:7-haida-gwaii-leaflet-caption)"
  knitr::include_graphics("images/07-Singlepart-geometry-of-the-Haida-Gwaii-archipelago-off-the-west-coast-of-British-Columbia-Canada-static.png")
} else { 
  fig_cap <- paste0("Singlepart geometry of the Haida Gwaii archipelago off the west coast of British Columbia, Canada. Hover over the islands to see the names. Polygon data from [@statistics_canada_statistics_nodate] and island placenames from [@natural_resources_canada_natural_nodate]. <a href='",canadalicense,"'>Open Government License - Canada</a>.")
m2
}
```
<br/>

Although it is possible to convert from a multipart geometry into singlepart geometry, you need to carefully consider how your features should be represented in the attribute table. For example, if you will be undertaking calculations using area or perimeter of the constituent polygons that comprise a multipart geometry of Canada, then you will return a single value for all of Canada while singlepart geometry would return values for each individual polygon. As well, area calculations can vary between multipart and singlepart geometry. For example, approximately 27% of Canada's land area (including freshwater), is comprised of more than 52,000 islands, which is a statistic you could only calculate with singlepart geometry. Thus, the choice of representing a feature using singlepart or multipart geometry should be based on how the features will be used in your analysis (i.e., aggregated versus disaggregated).

### Holes

When dealing with polygon features, **holes** may occur, which represent discontinuity of the interior polygon space. Imagine the case of a forested land cover that surrounds a lake. If we consider the forested land cover polygon on its own, then the polygon will have a hole where the lake exists (Figure \@ref(fig:7-hole)).

``` {r 7-hole, echo=FALSE, out.width = '50%', fig.cap=fig_cap}
fig_cap <- paste("Conceptual forest land cover polygon that contains a lake causing a hole. Pickell, CC-BY-SA-4.0.")
knitr::include_graphics("images/07-hole.png")
```

Topologically, holes in polygons imply that another polygon shares an adjacent boundary where the hole exists, for example, from the union of two layers (see [Chapter 6](https://ubc-geomatics-textbook.github.io/geomatics-textbook/overlay-and-proximity-analysis.html)). In our example, the lake would comprise its own polygon that would completely fill the hole.

### Delaunay triangulation

**Delaunay triangulation** is method for forming a triangle mesh over a set of points. The Delaunay triangulation method connects all points in a set such that no point in the set lays _within_ a circumcircle formed by any of the triangles in the mesh [@delaunay_sur_1934]. A circumcircle is a circle that passes through all the vertices of a cyclic polygon such as a triangle. In other words, the circumcircles are empty. To illustrate this, consider the four points in Figure \@ref(fig:7-delaunay-triangulation-circumcircles). There are only two circumcircles that can be formed from this set of points that ensures that no point lays within a circumcircle. The triangulation is then simply the lines connecting the three points that fall on any given circumcircle. One important property of the Delaunay triangulation is that the smallest angle in the resulting triangles is maximized from the circumcircle fitting, which minimizes sliver triangles that might from with very shallow angles. 

``` {r 7-delaunay-triangulation-circumcircles-setup, echo=FALSE, cached=T}
## Not run: creates image
## set.seed(7)
## x <- runif(4)
## y <- runif(4)
## tesselation <- deldir(x,y)
## png("/images/07-delaunay-triangulation-circumcircles.png",height=600,width=600)
## par(mar=c(0,0,0,0))
## plot(x,y,xlim=c(-0.25,1.25),ylim=c(-0.25,1.25),pch=19,xaxt="n",yaxt="n",ylab="",xlab="",bty="n")
## c1 <- circum(x[2:4],y[2:4])
## circles(c1$x,c1$y,c1$radius,col="grey")
## points(c1$x,c1$y,pch=19,col="red")
## c2 <- circum(x[c(1:3)],y[c(1:3)])
## circles(c2$x,c2$y,c2$radius,col="grey")
## points(c2$x,c2$y,pch=19,col="red")
## plot.deldir(tesselation,add=T,xlim=c(-0.25,1.25),ylim=c(-0.25,1.25),cmpnt_col=c("black","red","black","black","black"),pch=19,bty="n")
## dev.off()
```

``` {r 7-delaunay-triangulation-circumcircles, echo=FALSE, out.width = '75%', fig.cap=fig_cap}
fig_cap <- paste("Delaunay triangulation of four points. Black lines show the triangulation, grey lines represent the circumcircles connecting the three points of each triangle, red points represent the centres of the circumcircles, and the red dotted lines show that connecting the centres of the circumcircles forms the Voronoi diagram. Pickell, CC-BY-SA-4.0.")
knitr::include_graphics("images/07-delaunay-triangulation-circumcircles.png")
```

Figure \@ref(fig:7-delaunay-triangulation) shows a Delaunay triangulation for a set of 50 points. We can see that sliver triangles mostly occur on the edge of the extent of the points. Delaunay triangulations can be performed both in 2- and 3-dimensional Euclidean space and are therefore important for representing 3D surfaces as well as performing spatial estimation over 2D areas from a set of points.

``` {r 7-delaunay-triangulation-setup2, echo=FALSE, cached=T}
## Not run: creates image
## set.seed(14)
## x <- runif(50)
## y <- runif(50)
## tesselation <- deldir(x,y)
## tiles <- tile.list(tesselation)
## png("/images/07-delaunay-triangulation.png",height=600,width=600)
## plot(tiles,xlim=c(-0.25,1.25),ylim=c(-0.25,1.25),border="white")
## plot.deldir(tesselation,add=T,wlines="triang",xlim=c(-0.25,1.25),ylim=c(-0.25,1.25),cmpnt_col=c("black",NA,NA,"black","black"),pch=19)
## dev.off()
```

``` {r 7-delaunay-triangulation, echo=FALSE, out.width = '75%', fig.cap=fig_cap}
fig_cap <- paste("Delaunay triangulation of 50 random points. Pickell, CC-BY-SA-4.0.")
knitr::include_graphics("images/07-delaunay-triangulation.png")
```

### Thiessen polygons

**Thiessen polygons** are an implementation of a nearest neighbour algorithm in Euclidean space: given some set of input point features mapped on a plane, partition the plane into polygon areas that represent the nearest locations on the plane to those points. These resulting polygons are also sometimes referred to as proximal polygons, representing the proximal areas given some set of points. When Thiessen polygons are created for geographic data, the resulting diagrams are called **Voronoi diagrams** and sometimes referred to as Voronoi maps (Figure \@ref(fig:7-2d-thiessen-polygons)). Voronoi maps have many uses such as partitioning geographic space into areas that are nearest to weather stations, airports, or cellular towers. Thiessen polygons can be intersected with other geographic data layers in a GIS using map algebra to efficiently solve proximal questions like, "what is the nearest X?" without having to search or calculate the exact distances of all nearby features, which can be computationally time-consuming [@okabe_nearest_1994]. 

``` {r 7-create-thiessen-polygons, echo=FALSE, cached=T}
## Not run: creates image
## set.seed(14)
## x <- runif(50)
## y <- runif(50)
## tesselation <- deldir(x,y)
## tiles <- tile.list(tesselation)
## png("/images/07-2d-thiessen-polygons.png",height=600,width=600)
## plot(tiles, pch=19, fillcol=hcl.colors(50, "Purple-Yellow"),xlim=c(-0.25,1.25),ylim=c(-0.25,1.25))
## dev.off()
```

``` {r 7-2d-thiessen-polygons, echo=FALSE, out.width = '75%', fig.cap=fig_cap}
fig_cap <- paste("Thiessen polygons of 50 random points. Pickell, CC-BY-SA-4.0.")
knitr::include_graphics("images/07-2d-thiessen-polygons.png")
```

Thiessen polygons are a product of Delaunay triangulation described in the previous section. Figure \@ref(fig:7-delaunay-triangulation-thiessen-polygons) shows the relationship between the points, triangulation, circumcircles, and the Thiessen polygons. Connecting the circumcentres of the circumcircles produces the Voronoi diagram (Figure \@ref(fig:7-delaunay-triangulation-thiessen-polygons)).

``` {r 7-delaunay-triangulation-thiessen-polygons-setup, echo=FALSE, cached=T}
## Not run: creates image
## set.seed(14)
## x <- runif(50)
## y <- runif(50)
## tesselation <- deldir(x,y)
## tiles <- tile.list(tesselation)
## png("/images/07-delaunay-triangulation-thiessen-polygons.png",height=600,width=600)
## plot(tiles, pch=19, fillcol=hcl.colors(50, "Purple-Yellow"),xlim=c(-0.25,1.25),ylim=c(-0.25,1.25))
## plot.deldir(tesselation,add=T,wlines="both",xlim=c(-0.25,1.25),ylim=c(-0.25,1.25),cmpnt_col=c("red","black","black","black","black"),pch=19,bty="n")
## dev.off()
```

``` {r 7-delaunay-triangulation-thiessen-polygons, echo=FALSE, out.width = '75%', fig.cap=fig_cap}
fig_cap <- paste("Delaunay triangulation (red lines) overlaid onto the Thiessen polygons. Pickell, CC-BY-SA-4.0.")
knitr::include_graphics("images/07-delaunay-triangulation-thiessen-polygons.png")
```

### Centroids

A **centroid** is a point that represents the geometric centre of a polygon. For convex polygons, the centroid will always lay within the polygon, but for concave polygons, the centroid may lay outside the polygon (Figure \@ref(fig:7-concave-polygon-centroid)). Circular polygons always have centroids that are equidistant to the boundary of the polygon (Figure \@ref(fig:7-circle-polygon-centroid)).

``` {r 7-concave-polygon-centroid, echo=FALSE, out.width = '50%', fig.cap=fig_cap}
fig_cap <- paste("Concave polygon with the centroid (red dot) laying outside its boundary. Pickell, CC-BY-SA-4.0.")
knitr::include_graphics("images/07-concave-polygon-centroid.png")
```

``` {r 7-circle-polygon-centroid, echo=FALSE, out.width = '50%', fig.cap=fig_cap}
fig_cap <- paste("Circle polygon with the centroid (red dot) laying equidistant from the boundary of the polygon. Pickell, CC-BY-SA-4.0.")
knitr::include_graphics("images/07-circle-polygon-centroid.png")
```

### Convex hull

A **convex hull** is the smallest polygon that contains some set of points. It is sometimes also referred to as a "convex envelope" or "convex closure" because the perimeter of the polygon is formed by connecting the outermost points and closing or enveloping the remaining points. The convex hull is therefore the mathematical implementation of **topological closure**, where closure refers to the smallest closed set of points that contain the set of points. In practice, the convex hull is a bit like applying a rubber band around the outermost points so that the tension of the rubber band forms straight lines between the pairs of points in the closed set (Figure \@ref(fig:7-2d-convex-hull)). There are several algorithms for computing the convex hull, including Jarvis' March [@jarvis_identification_1973], Graham's Scan [@graham_efficient_1972], quickhull [@barber_quickhull_1996], and CudaHull [@stein_cudahull_2012].

``` {r 7-create-convex-hull, echo=FALSE}
## Not run: creates image
## set.seed(14)
## x <- runif(50)
## y <- runif(50)
## xy <- cbind(x,y)
## hpts <- chull(x,y)
## hpts <- c(hpts, hpts[1])
## png("/images/07-2d-convex-hull.png",height=600,width=600)
## plot(x,y,pch = 19,xlab="x-coordinate",ylab="y-coordinate",xlim=c(-0.1,1.1),ylim=c(-0.1,1.1))
## lines(xy[hpts,],lwd=2)
## dev.off()
```

``` {r 7-2d-convex-hull, echo=FALSE, out.width = '75%', fig.cap=fig_cap}
fig_cap <- paste("Convex hull formed by topological closure of the smallest closed set of points around the entire set of points. Arrangement of the points are the same as in the Thiessen polygons figure above. Pickell, CC-BY-SA-4.0.")
knitr::include_graphics("images/07-2d-convex-hull.png")
```

### Convex alpha hulls and alpha shapes

Convex hulls can be generalized to the concave case, called **convex alpha hulls** or **α-shapes (alpha shapes)**, by adjusting the maximum radius of the circumcircles through a parameter, alpha $α$. The objective of a convex alpha hull is to minimize the α-shape formed by circumcircles of radius less than or equal to $α$. Similar to the Delaunay triangulation, the circumcircles must be _open_, meaning they contain no other points in the set. Defined in this way, the final α-shape may not result in closure of the full set of points and can result in holes where the distance between points exceeds $2α$. Surprisingly, α-shapes are prone to not existing at all. For the case of $α=0$, applying circumcircles of radius 0 results in an empty α-shape and only the input set of points are returned without any boundaries. For $α=∞$, the α-shape is equivalent to the convex hull because if we use circumcircles with an infinitely large radius, then all points in the set are bound to be enclosed by the resulting α-shape, which like the convex hull must also minimize the bounding area. Figure \@ref(fig:7-2d-alpha-hull) shows an animation of the α-shapes for $α=1$ to $α=0$ for our set of 50 points.

``` {r 7-create-alpha-hull, echo=F}
## Not run; creates images for gif
## set.seed(14)
## x <- runif(50)
## y <- runif(50)
## xy <- cbind(x,y)
## alphas <- seq(1, 0, -0.05)
## alphas.str <- as.character(format(round(alphas,digits=2)),nsmall=2)
## for (a in 1:length(alphas)){
##   hpts <- ahull(xy,alpha=alphas[a])
##   png(paste0("c:/gif/alpha-",a,".png"), width=600,height=600)
##   plot(hpts,pch=19,xlab="x-coordinate",ylab="y-coordinate",xlim=c(-0.25,1.25),ylim=c(-0.25,1.25),main=paste0("Alpha = ",alphas.str[a]))
##   box(which = "plot", lty = "solid")
##   dev.off()
## }
```

``` {r 7-2d-alpha-hull, echo=FALSE, out.width = '75%', fig.cap=fig_cap}
fig_cap <- paste0("Concave alpha hull generates an alpha shape around a set of points. Online version of the figure is animated by alpha values from 0 to 1 by increments of 0.05. <a href='",baseurl,"#7-2d-alpha-hull'>Animated figure can be viewed in the web browser version of the textbook.</a> Pickell, CC-BY-SA-4.0.")
if (knitr:::is_latex_output()) { 
  knitr::include_graphics("images/07-2d-alpha-hull-static.png")
} else { 
  knitr::include_graphics("images/07-2d-alpha-hull.gif")
}
```

## 3D topologies

In this next section, a number of topologies are described that are important for 3D modeling and several examples are given using LiDAR (Light Detection and Ranging), which is the topic of [Chapter 15](https://ubc-geomatics-textbook.github.io/geomatics-textbook/LiDAR-acquisition-and-analysis.html). It is beyond the scope of this chapter to discuss the technology of LiDAR, so the reader is referred to [Chapter 15](https://ubc-geomatics-textbook.github.io/geomatics-textbook/LiDAR-acquisition-and-analysis.html) for a more in-depth discussion of LiDAR.

### Multipatch geometries

Similar to multipart geometries, **multipatch geometries** associate several *faces* or *facets* to a single 3D feature such as a building or tree. In order for multipatch geometries to be topologically valid, they must form a closed set of faces, known as a **polyhedron**. Polyhedrons are comprised of *flat* faces that connect 3 or more vertices. Figure \@ref(fig:7-platonic-polyhedrons) illustrates the 5 Platonic polyhedrons, so-named after Plaot who initially wrote about them. The Platonic polyhedrons are a special type of *regular* polyhedron because they are the only polyhedrons that are highly symmetrical and have special transitive properties on the edges, faces, and vertices. As well, the Platonic polyhedrons are all examples of the 3-dimensional case of a convex hull (more on that in the next section). Most polyhedrons that we come across in environmental management like trees, lakes, glaciers, and buildings are very irregular and not Platonic.

(ref:7-platonic-polyhedrons-caption) The five Platonic solids are examples of regular, convex polyhedrons and multipatch geometries. From left to right: tetrahedron (4 faces); hexahedron (6 faces); octahedron (8 faces); dodecahedron (12 faces); and icosahedron (20 faces). [@cyp_polyhedra_2005], CC-BY-SA-3.0. Animated figure can be viewed in the web browser version of the textbook: https://www.opengeomatics.ca/topology.html#7-platonic-polyhedrons.

```{r 7-platonic-polyhedrons, fig.cap = fig_cap, out.width= "75%", echo = FALSE}
if (knitr:::is_latex_output()) { 
  fig_cap <- "(ref:7-platonic-polyhedrons-caption)"
knitr::include_graphics("images/07-platonic-polyhedrons-static.png")
} else { 
  fig_cap <- paste0("The five Platonic solids are examples of regular, convex polyhedrons and multipatch geometries. From left to right: tetrahedron (4 faces); hexahedron (6 faces); octahedron (8 faces); dodecahedron (12 faces); and icosahedron (20 faces). [@cyp_polyhedra_2005], CC-BY-SA-3.0.")
  knitr::include_graphics("images/07-platonic-polyhedrons.gif")
}
```
<br/>

### 3D Convex hull

The 3D convex hull is the smallest convex polyhedron that contains a set of 3D points. The 3D problem is not unlike the 2D problem for finding the 2D convex hull, except instead of using 2D circumcircles we use 3D circumscribed spheres. Otherwise, the overall objective is the same, minimize the 3D polyherdron that encloses the 3D set of points. The 3D convex hull is frequently produced in order to generate a 3D object from a laser scan. Since points are 1-dimensional they have limited use beyond their enumeration within a volume or on a plane or on a line. By contrast, a 3D convex hull produces a polyhedron, which can be used to visualize the object that was initially scanned into a 3D point cloud. The 2D polygon faces that comprise a polyhedron can provide _shape_ and interact with a simulated light source to improve perception of _depth_, which are qualities that are not provided by 3D points alone.  

### 3D Convex alpha hull

The 3D convex alpha hull is Figure \@ref(fig:7-3d-alpha-hull-deciduous) below shows a 3D convex alpha hull for a deciduous tree near the Forest Sciences Centre at the University of British Columbia.

```{r 7-setup-3d-alpha-hull-deciduous, echo=FALSE, cached=T, warning=F, include=F}
las <- lidR::readLAS("Data/07/07-deciduous-tree.las")
z <- las@data$Z
x <- las@data$X
y <- las@data$Y
x <- (x-min(x))/(max(x)-min(x))
y <- (y-min(y))/(max(y)-min(y))
z <- (z-min(z))/(max(z)-min(z))
xyz <- cbind(x,y,z)
ahull.d <- ashape3d(xyz,alpha=0.05)
```

(ref:7-3d-alpha-hull-deciduous-caption1) 3D concave alpha hull for a deciduous tree. The alpha hull was generated using alpha = 0.05. Data collected by Spencer Dakin Kuiper with a GeoSlam terrestrial laser scanner in Vancouver, Canada. Pickell, CC-BY-SA-4.0. Animated figure can be viewed in the web browser version of the textbook: https://www.opengeomatics.ca/topology-and-geocoding.html#fig:7-3d-alpha-hull-deciduous.

(ref:7-3d-alpha-hull-deciduous-caption2) 3D concave alpha hull for a deciduous tree. The alpha hull was generated using alpha = 0.05. Data collected by Spencer Dakin Kuiper with a GeoSlam terrestrial laser scanner in Vancouver, Canada. Pickell, CC-BY-SA-4.0.

```{r 7-3d-alpha-hull-deciduous, echo=FALSE, out.width = '75%', fig.cap=fig_cap, webgl = TRUE}
if (knitr:::is_latex_output()) { fig_cap = "(ref:7-3d-alpha-hull-deciduous-caption1)"
knitr::include_graphics("images/07-3d-alpha-hull-deciduous-static.png")
} else { fig_cap = "(ref:7-3d-alpha-hull-deciduous-caption2)"
plot3d(ahull.d,col="#2D5A27")
rglwidget(width = 600, height = 600, reuse=FALSE, elementId="7-3d-alpha-hull-deciduous")
}
```

``` {r 7-setup-3d-alpha-hull-conifer, echo=F, cached=T, warning=F, include=F}
las <- readLAS("Data/07/07-conifer-tree.las")
z <- las@data$Z-1156
x <- las@data$X-640600
y <- las@data$Y-5674900
xyz <- cbind(x,y,z)
ahull.c <- ashape3d(xyz,alpha=0.1)
```

(ref:7-3d-alpha-hull-conifer-caption1) 3D concave alpha hull for a coniferous tree. The alpha hull was generated using alpha = 0.05. Data collected by Spencer Dakin Kuiper with a GeoSlam terrestrial laser scanner in Vancouver, Canada. Pickell, CC-BY-SA-4.0. Animated figure can be viewed in the web browser version of the textbook: https://www.opengeomatics.ca/topology-and-geocoding.html#fig:7-3d-alpha-hull-conifer.

(ref:7-3d-alpha-hull-conifer-caption2) 3D concave alpha hull for a coniferous tree. The alpha hull was generated using alpha = 0.05. Data collected by Spencer Dakin Kuiper with a GeoSlam terrestrial laser scanner in Vancouver, Canada. Pickell, CC-BY-SA-4.0.

``` {r 7-3d-alpha-hull-conifer, echo=FALSE, out.width = '75%', fig.cap=fig_cap, webgl = TRUE}
if (knitr:::is_latex_output()) { fig_cap = "(ref:7-3d-alpha-hull-conifer-caption1)"
knitr::include_graphics("images/07-3d-alpha-hull-conifer-static.png")
} else { fig_cap = "(ref:7-3d-alpha-hull-deciduous-caption2)"
plot3d(ahull.c,xlim=c(40,64),ylim=c(21,45),col="#2D5A27")
rglwidget(width = 600, height = 600, reuse=FALSE, elementId="7-3d-alpha-hull-conifer")
}
```

We can derive some useful information from these 3D convex alpha hulls. One obvious use of a polyhedron such as this is to calculate the volume, which for a tree might inform the amount of merchantable wood in the stem or the size of the canopy. Also, the polyhedron has replaced the points with polygon faces, which can be used to model shadowing and shading by simulating the position of the Sun in the sky. 

### 3D Voronoi tessellation

Similar to Thiessen polygons, Voronoi tessellations can be undertaken for 3D point clouds. 3D Voronoi tessellations produce nearest neighbour polyhedrons around the space of each 3D point. Just as in the 2D case, these 3D Voronoi tessellations can be used to partition the 3D space. Figure \@ref(fig:7-neper) shows the use of a 3D Voronoi tessellation for use in mapping and visualizing polycrystals.

(ref:7-neper-caption) 3D Voronoi tesselation of simulated polycrystals. [@quey_neper_nodate], GPL 3.0.

```{r 7-neper, fig.cap = fig_cap, out.width= "75%", echo = FALSE}
if (knitr:::is_latex_output()) { 
  fig_cap <- "(ref:7-neper-caption)"
} else { 
  fig_cap <- paste("3D Voronoi tesselation of simulated polycrystals. @quey_neper_nodate, GPL 3.0.")
}
knitr::include_graphics("images/07-neper.png")
```
<br/>

## Reflection Questions {-}

1. Give some examples of situations where you might use planar and non-planar topology.
2. Define topology.
3. How are Thiessen Polygons related to the Delaunay Triangulation?
4. How do you interpret a Voronoi diagram? Give an example for how you could use this layer to solve a spatial problem.
5. What are some uses of a multipart geometry?

## Practice Questions {-}

1. Draw a convex hull for the points below. Estimate where the centroid is for the polygon you drew.
2. Draw the Delaunay Triangulation for the points below.
3. Draw two lines from the points below that share a pseudonode.
4. Draw a polygon with a hole from the points below.

```{r 7-scatterplot, echo=FALSE, out.width = '75%', fig.cap=fig_cap}
fig_cap <- paste("Ten random points. Pickell, CC-BY-SA-4.0.")
knitr::include_graphics("images/07-scatterplot.png")
```

## Recommended Readings {-}

Barber, C.B., Dobkin, D.P., and Huhdanpaa, H. 1996. The quickhull algorithm for convex hulls. ACM Transactions on Mathematical Software, 22(4): 469-483. https://doi.org/10.1145/235815.235821

Delaunay, B.N. 1934. Sur la Sph`ere Vide. Izvestia Akademia Nauk SSSR, VII Seria, Otdelenie
Matematicheskii i Estestvennyka Nauk, 7:793–800.

Graham, R.L. 1972. An Efficient Algorithm for Determining the Convex Hull of a Finite Planar Set. Information Processing Letters, 1(4): 132–133. doi:10.1016/0020-0190(72)90045-2.

Jarvis, R.A. 1974. On the identification of the convex hull of a finite set of points in the plane. Information Processing Letters, 2(1):18-21. https://doi.org/10.1016/0020-0190(73)90020-3

Okabe, A., Boots, B., and Sugihara, K. 2007. Nearest neighbourhood operations with generalized Voronoi diagrams: a review. International Journal of Geographical Information Systems, 8(1): 43-71. https://doi.org/10.1080/02693799408901986

Stein, A., Geva, E., and El-Sana, J. 2012. CudaHull: Fast parallel 3D convex hull on the GPU. Computers and Graphics, 36(4): 265-271. https://doi.org/10.1016/j.cag.2012.02.012