---
title: "High Quality Constrained Triangulation for Simple Features"
author: "Michael D. Sumner"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    fig_width: 7
    fig_height: 7
vignette: >
  %\VignetteIndexEntry{"High Quality Constrained Triangulation for Simple Features}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

## Example

This is a basic example which shows you how to decompose a MULTIPOLYGON `sf` data frame object into a GEOMETRYCOLLECTION `sf` data frame object made of triangles:

```{r example}


library(sf)
library(sfdct)
nc <- read_sf(system.file("shape/nc.shp", package="sf"), quiet = TRUE)
nc_triangles <- ct_triangulate(nc)

plot(st_geometry(nc_triangles),   col = viridisLite::viridis(nrow(nc_triangles)))
```

We can use the underlying `RTriangle::triangulate` arguments to hone the triangles we get. 

```{r hone}
i_feature <- 25
nc1 <- nc[c(i_feature, unlist(st_touches(nc[i_feature, ], nc))), ]
plot(st_geometry(nc1),col = viridisLite::viridis(nrow(nc1)))

## subvert st_area because we really don't want m^2
st_crs(nc1) <- NA
areas <- st_area(nc1)
st_crs(nc1) <- st_crs(nc)
nc1_triangles <- ct_triangulate(nc1, a = min(areas)/5)
bcol <- viridisLite::viridis(nrow(nc1_triangles))
plot(st_geometry(nc1_triangles), col = NA, border = bcol)

nc2_triangles <- ct_triangulate(nc1, a = min(st_area(st_set_crs(nc1, NA)))/25)
plot(st_geometry(nc2_triangles), col = NA, border = bcol)

```

Get a grouped triangulated set from a MULTIPOINT. Note how these aren't constrained by 
the edges of the input polygons (because we threw those away!) but these are controlled to have a smaller maximum area. 

Area is calculated in the native coordinates, assuming "planar coordinates", with no respect to the real world. 

```{r MULTIPOINT}
## manual cast to MULTIPOINT originally required
#st_geometry(nc1) <- st_sfc(lapply(unlist(unlist(st_geometry(nc1), recursive = FALSE), recursive = FALSE), st_multipoint), crs = st_crs(nc1))
mp_nc1 <- st_cast(nc1, "MULTIPOINT")
mtriangs <- ct_triangulate(nc1, a = 0.0005)
plot(st_geometry(mtriangs), col = viridisLite::viridis(nrow(mtriangs)), border = "#00000033")
```

```{r}
plot(nc[4, ]$geometry)
## q, minimum angle
## D, Delaunay criterion is met
plot(ct_triangulate(nc[4, ]$geometry, q = 35, D = TRUE), add = TRUE, col = "transparent")

```

## Geometry collection is in development

POLYGON triangles in GEOMETRYCOLLECTION will be re-triangulated. All vertices in the GC will be included, as well as all 
edges of all component geometries, but each component is triangulated individually, not with reference to the entire set. 



## Chaining together operations 

We can use piping to chain things together. 

```{r}
data("map_world", package= "sfdct")
library(dplyr)
g <-  map_world %>% dplyr::filter(startsWith(ID, "Indonesia")) %>% ct_triangulate() %>% st_geometry() 
plot(g, col = "aliceblue", main = "")

nc_triangles[1:2, c(1, 5)] %>% st_transform("+proj=laea") %>% ct_triangulate(a = 2e6) %>% plot()
```

## Holes better work (or else)

```{r}
data("antarctica")
plot(antarctica[0])
a <- ct_triangulate(st_difference(antarctica[1], antarctica[2, ]), a = 5e10)
plot(a[0], col = "firebrick")
```


The output of `ct_triangulate` can be the input to another call to it. 

```{r}
m <- map_world %>% dplyr::filter(startsWith(ID, "Papua New Guinea"))
plotme <- function(x) {plot(st_geometry(x), col = "aliceblue"); x}
m %>% ct_triangulate(a = 0.2, D = TRUE) %>% plotme() 
```