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.
trails <- macleish_layers %>%
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) %>%
units::set_units("miles")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?
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.proj4_aea <- "+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"boundaries_aea <- macleish_layers %>%
pluck("boundary") %>%
st_transform(proj4_aea)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.8564224set_units() function to do the conversion.# SAMPLE SOLUTION
boundaries_aea %>%
st_area() %>%
units::set_units("acres")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"))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.
trail_stream_xings <- macleish_layers %>%
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 only way to do this. This is just one way that makes sense to me.
stream_buffer <- macleish_layers %>%
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:
my_points with columns for lat and lon.my_points <- tribble(
~point, ~lat, ~lon,
"A", 42.449285, -72.679370,
"B", 42.448988, -72.678093,
"C", 42.449455, -72.677567,
)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_sf <- my_points %>%
st_as_sf(coords = c("lon", "lat"))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?