In this lab, we will learn how to compute various geometric properties of spatial objects.

```
library(tidyverse)
library(macleish)
```

**Goal**: by the end of this lab, you will be able to use `sf`

to perform computations on spatial objects.

The GEOS library is the workhorse for computations on spatial objects. It can do things like compute areas, lengths, intersections, and other geometric properties.

Use the `st_length()`

function to compute the length of lines.

```
library(sf)
%>%
macleish_layers pluck("trails") %>%
st_length()
```

Note that the length computed above is in terms of the units of the current projection, which in this case is meters!

This returns a vector, so we could add the computed length as a new variable to our existing data frame. Note here that `geometry`

is always the name of the column that contains the spatial information in an `sf`

data frame.

```
<- macleish_layers %>%
trails pluck("trails") %>%
mutate(computed_length = st_length(geometry))
```

Note that the `computed_length`

variable in `trails`

is now a numeric vector of class `Units`

.

```
%>%
trails pull(computed_length) %>%
str()
```

You can now use the `set_unts()`

function from the `units`

package to easily convert from meters to, say, miles.

```
%>%
trails pull(computed_length) %>%
::set_units("miles") units
```

Here we compute the area of the MacLeish property.

Note that the area computed below has the units explicitly stated!

```
%>%
macleish_layers pluck("boundary") %>%
st_area()
```

`sf`

is smart enough to pay careful attention to units – as well you should! If we transform our spatial object into a different projection, will the area stay the same? Should it?

- Use the
`st_transform()`

function and the`proj4_aea`

string below to project the MacLeish boundary layers into an Albers Equal Area projection, which is**area-preserving**and in the units of**meters**.

`<- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" proj4_aea `

```
<- macleish_layers %>%
boundaries_aea pluck("boundary") %>%
st_transform(proj4_aea)
```

- Use
`st_area()`

to compute the area of the boundary in the Albers Equal Area projection and convert to**acres**manually using arithmetic. Note that 1 acre is equal to 4046.8564224 square meters.

```
# SAMPLE SOLUTION
%>%
boundaries_aea st_area() / 4046.8564224
```

- Repeat the previous exercise, but this time use the
`set_units()`

function to do the conversion.

```
# SAMPLE SOLUTION
%>%
boundaries_aea st_area() %>%
::set_units("acres") units
```

The purported size of the MacLeish property is 260 acres. Does this accord with your estimate?

Convert the MacLeish boundary into EPSG:3857.

Compute the area in acres of the MacLeish boundary in EPSG:3857. Is it close to 260 acres?

Use

`st_area()`

to compute the area of each of the polygons in the`forests`

layer. Add the corresponding results to the`forests`

layer`sf`

data frame using`mutate()`

. Which one is the largest?

```
# SAMPLE SOLUTION
%>%
macleish_layers pluck("forests") %>%
mutate(computed_area = st_area(geometry)) %>%
arrange(desc(computed_area))
```

You can also detect whether one layer intersects with another. The `trails`

layers contain spatial information about 15 trail segments, while the `streams`

layer contains spatial information about 13 stream segments. Do these layers intersect?

```
%>%
macleish_layers pluck("trails") %>%
st_intersects(pluck(macleish_layers ,"streams"))
```

Do the layers cross each other?

```
%>%
macleish_layers pluck("trails") %>%
st_crosses(pluck(macleish_layers, "streams"))
```

- In this case, there is no difference between the results of
`st_interects()`

and`st_crosses`

applied to the trail and streams layers. Will that**always**be the case? What circumstances would lead to the two functions returning different results?

If you actually want to know where the intersections occur, the `st_intersection()`

function can compute them and return an `sf`

object.

```
<- macleish_layers %>%
trail_stream_xings pluck("trails") %>%
st_intersection(pluck(macleish_layers, "streams"))
```

We can verify this by plotting the intersections on a map. We have a small hoop to jump through here because one of the intersections has type `MULTIPOINT`

, which `leaflet`

cannot handle. To get around this we can use `st_cast()`

to convert the `MULTIPOINT`

feature into `POINT`

.

```
library(leaflet)
leaflet() %>%
addTiles() %>%
addPolylines(data = pluck(macleish_layers, "streams")) %>%
addPolylines(data = pluck(macleish_layers, "trails"), color = "brown") %>%
addMarkers(data = st_cast(trail_stream_xings, "POINT"), popup = ~name)
```

You can use `st_buffer()`

to draw a buffer zone around a spatial object. Note that the `width`

argument to `st_buffer()`

is *in the units of the projection*, which in this case is **degrees**. To specify the buffer depth in **meters**, we have to convert the spatial data into a projection scheme that uses meters. Again, we’ll use the Albers Equal Area projection. Afterwards, we convert back to EPSG:4326 for plotting on Leaflet.

This is not necessarily the

onlyway to do this. This is just one way that makes sense to me.

```
<- macleish_layers %>%
stream_buffer pluck("streams") %>%
st_transform(proj4_aea) %>%
st_buffer(dist = 10) %>%
st_transform(4326)
leaflet() %>%
addTiles() %>%
addPolylines(
data = pluck(macleish_layers, "streams"),
weight = 1, color = "black"
%>%
) addPolygons(data = stream_buffer)
```

You may want to create new spatial objects. At a basic level all you need to do to accomplish this is to create a data frame with a `geometry`

column. [In practice, there may be more hoops to jump through.]

A feasible, if inelegant, way to create a new set of `sf`

object with `POINT`

features is as follows:

- Click on the point you want to use in GoogleMaps
- Right-click to see the coordinates
- Copy-and-paste those into R, and put them into a data frame called
`my_points`

with columns for`lat`

and`lon`

.

```
<- tribble(
my_points ~point, ~lat, ~lon,
"A", 42.449285, -72.679370,
"B", 42.448988, -72.678093,
"C", 42.449455, -72.677567,
)
```

- Use the
`st_as_sf()`

function to convert`my_points`

into an`sf`

object. Note that you will have to specify the variables that contain the spatial coordinates using the`coords`

argument.

```
# SAMPLE SOLUTION
<- my_points %>%
my_sf st_as_sf(coords = c("lon", "lat"))
```

- Put the points on a
`leaflet`

map.

```
leaflet() %>%
addTiles() %>%
addMarkers(data = my_sf, popup = ~point)
```

Please respond to the following prompt on Slack in the `#mod-spatial`

channel.

Prompt: What additional functions or abilities do you think you still need to complete Mini-Project #4?