Getting home

0.1 R Markdown

The three rings to rule the world in R and spatial analysis are sf, tigris, and ggplot or leaflet. OK, it is also easier to take over places in the US if you also have the tidycensus package. But first, let’s look at some bird data to go through each of these super powers.

https://gist.github.com/edzer/f461a3a95570c4ab7edf3125c2f19d20

Say you want to measure how far a bird must fly from their nest in St. Paul, MN down to the the center of Smoky Mountain National Park in Tennesse. Let’s set up 2 points, and then measure how far apart they are. Next we will develop a polygon that is a straight path between these two nesting locations and join with the states to develop a final list of states that the birds will fly through and make some maps. We will of course follow up this exercise by setting up strict conservation laws to protect these flight paths.

A tree in the center of St. Paul.

library(sf)

lat_y <- 44.9283966 #Always remember that latitude is a ladder.
long_x <- -93.1734448

coords <- c(long_x, lat_y)

nest_mn <- st_sfc(st_point(coords))

nest_mn <- st_as_sf(nest_mn, crs = 4326)

A tree in the center of Smoky Mountain National Park

lat_y <- 35.5812 # Remember that latitude is a ladder going up and down
long_x <- -83.8609

coords <- c(long_x, lat_y)

nest_smnp <- st_sfc(st_point(coords)) # setting up the datum

nest_smnp <- st_as_sf(nest_smnp, crs = 4326) # providing the coordinate reference system

How far part are these two nests? We need to convert the coordinate reference systems to make sure that the two locations are comparable and our final answer is in meters.

nest_mn_meters <- st_transform(nest_mn, crs = 9822)

nest_smnp_meters <- st_transform(nest_smnp, crs = 9822)

tree_distances <- st_distance(nest_mn_meters, nest_smnp_meters)

tree_distances <- tree_distances/1000

Each nest’s coordinates are currently in latitude and longitude, so we will first convert to UTMs that have units that are meters.

My fave CRS info sheet

We will also make a data frame with our two polygons. This isn’t necessary, but you may have to do this some day. Then, let’s say that the birds fly in roughly a straight line. So, let’s cast these two points into a line between them.

tree_points <- rbind(nest_mn, nest_smnp)

tree_route <-  tree_points %>%
    summarize(geometry = st_union(.)) %>%
    st_cast("LINESTRING") %>%
  mutate(route = "Minnesota to Smoky Mountain National Park")
    
plot(tree_route$geometry)

But we are fairly sure the birds won’t fly in an exactly straight line from Minnesota to Tennessee. So, let’s put a buffer around this line and make it a polygon. Let’s assume the birds will fly within 1000 meters on either side of our line. So, REMEMBER about units? Let’s check what our units currently are (hint latitude and longitude have what units?), and then change the units and draw our new polygon. It will be glorious.

st_crs(tree_route)

tree_route <- st_transform(tree_route, crs = 9822)

st_crs(tree_route)

tree_route <- st_buffer(tree_route, dist = 1000)

plot(tree_route)

tree_route <- st_transform(tree_route, crs = 4326)

plot(tree_route)

Looking at that shape in our Plots window is satisfying, but doesn’t tell us anything about what states the polygon goes through, and we’d like a pretty map in the end. So, let’s get started. First, we are going to pull in our tigris package and grab the U.S. states’ polygons.

https://github.com/walkerke/tigris

library(tigris)

all_states <- states()

Now let’s see which states intersect with our polygon bird path. We bring back some sf magic for this. But remember, CRS must be the same!! SO, we will need to do some transformation.

st_crs(all_states)
st_crs(tree_route)

tree_route <- st_transform(tree_route, crs = st_crs(all_states))

bird_route_states <- st_filter(all_states, tree_route)

Now we have two types of polygons, a set of states that intersect with our bird path and the bird flight path itself. We will make a map of these two shapes in two ways. First with leaflet then with ggplot2.

First, leaflet

library(leaflet)

leaflet() %>% 
  addTiles() %>% 
  addPolygons(data = tree_route, weight = 8, fillColor = 'darkorange', color = 'darkorange') %>%
  addPolygons(data = bird_route_states)

Now, ggplot2

library(ggthemes)

ggplot() +
  geom_sf(data = bird_route_states, aes(fill = NAME, label = NAME)) +
  geom_sf_text(data = bird_route_states, aes(label = NAME)) +
  geom_sf(data = tree_route, aes(fill = route), size = 10) +
  theme(legend.position = "none")

# I get an error if I skip the step of including data =

Now let’s protect some birds!