shapes.Rmd

---
title: "Shapes"
author: "Michael Sumner"
date: "20/09/2017"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width=7, fig.height = 7)
```

# Shapes

This is not a shapefile. 

It's a `SpatialPolygonsDataFrame`, rendered in a particular way using the R programming language. 

```{r}
library(aceecostats)
library(sp)
plot(aes_region_simple, col = aes_region_simple$colour)
```


```{r}
library(mapview)
mapview(aes_region_ll, col.regions = aes_region_ll$colour) # viridis::viridis(nrow(aes_region_ll)))
```

It can be rendered in various ways. 

```{r}
library(anglr)

mesh <- anglr(aes_region_simple)
plot(mesh)
rgl::rglwidget()
```

```{r}
library(anglr)
rgl::rgl.clear()
mesh <- anglr(aes_region_simple)
plot(globe(mesh))
rgl::rglwidget()
```

Finally, we can nowadays actually edit or create our own drawings interactively, but this is still 
an early stage. 

```{r}
library(mapedit)
```

There's no easy way into this, it's a DEEP DIVE. 


```{r}
library(sfdct)
aes_triangles <- ct_triangulate(st_as_sf(aes_region_simple))

plot(aes_triangles$geometry, col = NA)
```
# Simple shape

A map of Antarctica is pretty complicated, so let's take it back a little. 

```{r}
library(silicate)
data("minimal_mesh")
library(sf)
print(minimal_mesh)
plot(minimal_mesh)
```


This simple situation hides a lot of problems. 

## First problem

It's pretty complicated. 

```{r}
str(unclass(minimal_mesh$geom))
```

## Second problem

There are three shapes here, but only two values.  (One of the shapes is a hole). 

```{r}
minimal_mesh$a
```
## Third problem

There are shared-elements between these two objects, two coordinates and one edge. But how can we know this? 

```{r}
op <- par(xpd = NA)
plot(minimal_mesh, col = NA)
ring1 <- minimal_mesh$geom[[1]][[1]][[1]]
points(ring1)
text(ring1, label = sprintf("%s,%s", ring1[, 1], ring1[, 2]), pos = 2, adj = 3, col = "dodgerblue")
ring2 <- minimal_mesh$geom[[2]][[1]][[1]]
points(ring2)
text(ring2, label = sprintf("%s,%s", ring2[, 1], ring2[, 2]), pos = 4, adj = 3, col = "firebrick")
ring3 <- minimal_mesh$geom[[1]][[1]][[2]]
points(ring3)
text(ring3, label = sprintf("%s,%s", ring3[, 1], ring3[, 2]), pos = 1, adj = 3, col = "black")

par(op)
```


## Standard algorithms don't know about the shapes

```{r}
## two POLYGONs, composed of three rings p1+p2 and p3
p1 <- cbind(x = c(0, 0, 0.75, 1,   0.5, 0.8, 0.69, 0), 
            y = c(0, 1, 1,    0.8, 0.7, 0.6, 0,    0))
p2 <- cbind(x = c(0.2, 0.2, 0.3, 0.5, 0.5, 0.2), 
            y = c(0.2, 0.4, 0.6, 0.4, 0.2, 0.2))
p3 <- cbind(x = c(0.69, 0.8, 1.1, 1.23, 0.69), 
            y = c(0, 0.6, 0.63, 0.3, 0))
library(sf)
g <- data.frame(two_polygons = 1:2)
g[["geometry"]] <- st_sfc(st_polygon(list(p1, p2[nrow(p2):1, ])), st_polygon(list(p3)))
g <- st_as_sf(g)

## Delaunay triangulation of the polygon vertices, grouped by feature
gt <- g
st_geometry(gt)  <- st_triangulate(st_geometry(g))

op <- par(mfrow = c(1, 2))
plot(g) 
plot(gt, col = scales::alpha(c("blue", "yellow"), 0.3), main = "convex triangulation")
plot(st_geometry(g), add = TRUE, lwd = 4)
```



# This is simple features

- shapes are represented as paths so only planar polygonal shapes are possible
- the standard allows for XY[Z[M]] geometry,but this is not extensible  - no capacity to store data against component geometry elements
- no internal topology (no vertex-, edge-, or path-sharing).



## Return to the blocky triangles in original

- resize the triangles for a smoother curve

- elevate with bathymetry

- show the North Carolina example



# Software