sf: Simple features in R

class: right, middle, my-title, title-slide

.title[
# <code>sf</code>: Simple features in R
]
.author[
### Gavin Fay (acknowledgements Sean Hardison)
]

---

.medium .remark-code { /*Change made here*/
  font-size: 80% !important;
}

.tiny .remark-code { /*Change made here*/
  font-size: 60% !important;
}
</style>

# Today's sf foray

.pull-left[
* Introduction to simple features

* Interacting with simple features in R using `sf`

* `sf` geometry types

* Common `sf` operations and plotting

* Integration with the `tidyverse`

* Example analyses

]

.pull-right[

![r-sf](https://user-images.githubusercontent.com/520851/34887433-ce1d130e-f7c6-11e7-83fc-d60ad4fae6bd.gif)

]

.footnote[
  Pebesma, E., 2018. Simple Features for R: Standardized Support for Spatial Vector
  Data. The R Journal 10 (1), 439-446, https://doi.org/10.32614/RJ-2018-009
]

---

## Simple features: Models for things in space

What are simple features? Simple features are any object or thing that can be represented by a point in space and an attribute

For example, SMAST West can be represented as a simple feature. It's associated with a latitude longitude pair, and we can assign the attribute "building" to it.

---

## Simple features: Models for things in space

---
## Simple features: Models for things in space

We assigned the building a point in space, but simple features can also be described by multiple points, such as this polygon. We might then assign the attribute "green space" to the feature.

---
## Simple features: Models for things in space

A third feature could be this road and pier, which we describe spatially as a line.

---
## Simple features: Models for things in space

We can also describe all of the these features as a single simple feature, which is referred to as a geometrycollection.

---
## Simple features: A standard for spatial information

* A data structure with spatial `geometries` and non-spatial `attributes`
* Used in spatial databases and commercial GIS applications (e.g. PostGIS and ArcGIS)

.footnote[
https://r-spatial.github.io/sf/articles/sf1.html
]

---
## Simple features in R using `sf`

* Classic methods for dealing with spatial data in R are a pain

* `SpatialPointsDataFrame`, `SpatialLinesDataFrame`, `SpatialPolygonsDataFrame`...
  
--
  
* `sf` simplifies this experience while taking advantage of widely used standards

* In `sf-worksheet.Rmd`, run the `view_structure` code chunk

???
**Run view_structure code chunk**

nc is an sf object containing polygons for all counties in the state of north Carolina.

sf: A simple feature
sfc: simple feature geometry list-column

sfg: simple feature geometry

---
## Simple features in R using `sf`

An object of class `sf`:
<img src="lab-08/EDAB_images/sf-data-struc.png" width="70%" style="display: block; margin: auto;" />

* Features in the collection are described by rows in the data.frame

* A column named `geometry` describes the spatial aspect of each feature

* All features in the `sf` object have a coordinate reference system (CRS)
 * Describes a transformation from a 3D surface to a 2D plane (e.g. globe to paper).
 
* CRS is defined by `proj4string` and `epsg`
 
???

* Functions in `sf` all begin with `st_`

---
## Feature geometry types

---
## Geometries and attributes

.medium[
.remark-inline-code[

```r
plot(nc)
```
]
]

---
## Feature geometry types
* All geometry types are defined by groups of point coordinates

* We can make our own special feature geometry (`sfg`) using the `sf` syntax `sf_[geometry]`
  * For example, `sf_point(c(1,1))`

.footnote[
[source](https://r-spatial.github.io/sf/articles/sf1.html)
]

---

## Feature geometry types: `POLYGON` and `MULTIPOLYGON`
* `POLYGON` geometries are allowed to have one external ring and zero or more internal rings
* `MULTIPOLYGON` geometries contain > 1 non-nested polygons

* All individual polygons start and end with the same coordinate pair

---
## A special type of geometry: `GEOMETRYCOLLECTION`

* If a feature (i.e. a row) contains more than one geometry type, it is a `GEOMETRYCOLLECTION`

---
## Feature geometry types

`POINT`: numeric vector  
  * `c(x,y)`
  
`LINESTRING`: A numeric matrix with points in rows
  * `matrix(data = c(x1,y1,...), nrow = 2)`

`POLYGON`: A list of numeric matrices
  * `list(matrix(x1, y1,...,x1,y1), nrow  = 2)`
  
`MULTIPOINT`: A numeric matrix with points in rows
* `matrix(data = c(x1,y1,...), nrow = 2)`
  
`MULTILINESTRING`: A list of numeric matrices
 * `list(matrix(x1, y1,...), nrow  = 2)`

`MULTIPOLYGON`: A list of lists of matrices...
 * `list(list(matrix(x11, y11,...,x11,y11), nrow  = 2), list(matrix(x21, y21,...,x21,y21), nrow  = 2), ...)`

`GEOMETRYCOLLECTION`: A list of `sf` objects
 * `list(sf_object1, sf_object2, ...)`
 
???
Multipolygon class is special case when there are non-nested polygons

---

## Advantages to `sfg` structure

* `sf` objects are also of class `data.frame`, meaning `tidyverse` functions can be applied(!)

.pull-left[

.medium[
.remark-inline-code[

```r
library(dplyr)

wilmington <- nc %>% 
  filter(NAME == "New Hanover")

plot(wilmington)
```
]
]

* Note that all associated attributes are plotted for a given geometry

]

.pull-right[
<img src="lab-08b_files/figure-html/plot-label-out-1.png" style="display: block; margin: auto;" />
]

???
**Turn to Getting tidy with it code chunk**
---

## `ggplot2` ft. `sf`

* `ggplot2` uses the `geom_sf` function to visualize `sf` objects

.medium[
.remark-inline-code[

```r
ggplot() +
  geom_sf(data = nc, 
          aes(fill = AREA))
```
]

]

???
**Turn to ggplot**

---
## Working with `sf`: Reading data

.tiny[
.remark-inline-code[

```r
#Set relative path to data directory
epu.dir <- here::here("data","EPU_shapefile")

#Read in shapefile using sf
 {{ sf_shape <- 
  st_read(file.path(epu.dir,"EPU_extended.shp"),quiet = T) }}

epu <- sf_shape %>% dplyr::select(EPU)

ggplot() +
  geom_sf(data = epu, aes(fill = EPU))
```
]
]

<img src="lab-08b_files/figure-html/epu_shapes_out-1.png" width="40%" style="display: block; margin: auto;" />
---
## Raster data with the `stars` package

.tiny[
.remark-inline-code[

```r
challenger_deep <- raster::raster("../data/challenger_deep.nc") |>
  st_as_stars()

chal_deep <- st_point(c(142.20205,11.332417))

ggplot() + 
  geom_stars(data = challenger_deep) +
  scale_fill_gradientn(colours = terrain.colors(20)) +
  geom_sf(data  = chal_deep) +
  geom_sf_text(data = chal_deep, aes(label = "Challenger Deep (-10,928 m)"), 
               nudge_y = 0.05, nudge_x = 0.05) + 
  ggtitle("Challenger Deep")
```
]
]

```
## trying to read file: /Users/gfay/courses/mar536-biolstats2-s23/data/challenger_deep.nc
```

---
## Working with `sf`: Changing Coordinate Reference Systems

**Identify the existing CRS** with `st_crs`

.medium[
.remark-inline-code[

```r
st_crs(epu)
```
]
]

**Change the CRS** with `st_transform`

.medium[
.remark-inline-code[

```r
#Pass a proj4string
p4s <- "+proj=longlat +datum=NAD27 +no_defs"
st_transform(epu, p4s)

#Pass an EPSG code
epsg <- 4627
st_transform(epu, epsg)
```
]
]

---
## Working with `sf`: Calculating geometric distance to beer

.tiny[
.remark-inline-code[

```r
#Wilmington is here
wilmy <- nc %>% filter(NAME == "New Hanover") %>% st_centroid()

#Beer is here
wicked_weed <- st_geometry(st_point(c(-82.551440,35.591738)))
st_crs(wicked_weed) <- st_crs(nc)

#Draw a line between beer and Wilmington
a_long_ways_to_beer <- st_geometry( sf::st_linestring(
  matrix(rbind(st_coordinates(wilmy), st_coordinates(wicked_weed)), ncol = 2)) )

st_crs(a_long_ways_to_beer) <- st_crs(nc)

#Calculate the total distance with st_distance
* total_dist <- st_distance(wilmy, wicked_weed) 
```
]
]

---
## Working with `sf`: Crop to bounding box
.tiny[
.remark-inline-code[

```r
data.dir <- here::here("data")
load(file.path(data.dir, "topo_4ft.Rdata"))

ymax <-  41.6525; xmin <- -70.5975
ymin <-  41.6425; xmax <- -70.59

#Define a bounding box
{{ bbox <- st_bbox(c(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), 
                crs = st_crs(topography)) }}

#Turn it into an sf object for plotting
grid_box <- st_make_grid(bbox) %>% st_cast("POLYGON")

ggplot() +
  geom_sf(data = topography, aes(color = ELEVATION)) +
  geom_sf(data = grid_box, alpha = 0, color = "red") +
  ggtitle("Topography of Francis-Crane Wildlife Area")
```
]
]

---
## Working with `sf`: Crop to bounding box
.tiny[
.remark-inline-code[

```r
* cropped_topo <- topography %>% st_crop(bbox)

ggplot() +
  geom_sf(data = cropped_topo, aes(color = ELEVATION)) +
  geom_sf(data = grid_box, alpha = 0, color = "red") +
  ggtitle("Topography of Francis-Crane Wildlife Area")
```
]
]

---

## Working with `sf`: Crop to arbitrary polygon

.tiny[
.remark-inline-code[

```r
custom_polygon <- sf::st_read(file.path(data.dir, "Polygon_shapefile", "polygon_for_crop.shp"), 
                              quiet = T) %>%
  st_transform(st_crs(topography))

#Intersection
* topo_intersection <- topography %>% st_intersection(custom_polygon)

#Difference
* topo_difference <- topography %>% st_difference(custom_polygon)

#Plotting
int <- ggplot() +
  geom_sf(data = topo_intersection, aes(color = ELEVATION)) +
  guides(color = F) + ggtitle("Intersection")

difff <- ggplot() +
  geom_sf(data = topo_difference, aes(color = ELEVATION)) +
  guides(color = F) + ggtitle("Difference")

int + difff + plot_layout(ncol = 2)
```
]
]

---
## Working with `sf`: Spatial joins
.tiny[
.remark-inline-code[

```r
sf_shape <- st_read(file.path(epu.dir,"EPU_extended.shp"),quiet = T) 
sf_use_s2(FALSE)
#Get the union of all EPUs
* all_joined <- st_union(sf_shape)

#Create a grouping column to union New England EPU polygons
{{ some_joined <- sf_shape %>% mutate(Region = ifelse(EPU %in% c("GOM","GB","SS"),
                                "New England",
                                "Mid-Atlantic")) %>% 
  group_by(Region) %>% 
  summarise() }}

all_joined <- ggplot() + geom_sf(data = all_joined, color = "blue") +
  guides(color = F)

some_joined <- ggplot() + geom_sf(data = some_joined, aes(color = Region)) +
  guides(color = F)

all_joined + some_joined + plot_layout(ncol = 2)
```
]
]

---
## Putting it together with fake fish

---
## Exercise 3 - Pulling it all together

Create a multi-panel plot with the following 4 elements:
- a map of the 2015 spring NMFS bottom trawl survey catch rates for silver hake off in the Gulf of Maine
- a histogram of silver hake catch rate across tows
- a histogram of the distances of tows from New Bedford (41.636 N, 70.934 W)
- a scatterplot of the silver hake catch rate as a function of depth

_stretch goal: include the population density of states or counties on the land map_. 
_stretch goal: color code the catch rates by EPU_.

---
## Exercise 3 - Pulling it all together

Steps:

1. Read in the survey data. (`hake.csv`)
2. Identify the data for silver hake in the spring, in 2015.
3. Create ggplot objects containing each of the 4 plots
4. For the map:
- map Gulf of Maine
- add bathymetry (raster in `data/gom_bathy.rds`)
- plot points corresponding to positive silver hake tows, make the size of the plotting character proportional to the log(Biomass)
- add point indicating New Bedford
- add appropriate labels and legend
5. Rely on experience from previous labs for the other plots. For the distances, use `st_distance()`

---
## References

#### Packages

[Simple Features for R](https://r-spatial.github.io/sf/index.html)

[Spatiotemporal Arrays, Raster and Vector Datacubes (stars)](https://r-spatial.github.io/stars/)

[gstat](https://r-spatial.github.io/gstat/)

#### Resources, code, and data

[NOAA NEFSC Ecosystem Dynamics and Assessment Github](https://github.com/NOAA-EDAB/)

[Geocomputation with R](https://geocompr.robinlovelace.net/)

[R-spatial](https://www.rspatial.org)

[Mariana Trench Digital Elevation Model](https://catalog.data.gov/dataset/mariana-trench-bathymetric-digital-elevation-model)

[NOAA NEFSC Oceans and Climate Branch](https://www.nefsc.noaa.gov/epd/oceanography/)

[Falmouth, MA GIS Site](http://www.falmouthmass.us/199/Public-Interactive-GIS-Site)