Skip to content

Commit

Permalink
Updated class 16 slides
Browse files Browse the repository at this point in the history
  • Loading branch information
ldemaz committed Mar 20, 2024
1 parent ac8acb3 commit 8f2814a
Show file tree
Hide file tree
Showing 2 changed files with 66 additions and 347 deletions.
187 changes: 9 additions & 178 deletions docs/class16.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,15 @@ districts %>%
summarize(mu = mean(areadiff)) #%>%
```

---
## Calculate length (exercise)
- Use `st_length()` to add column 'road_length' to `roads`
- 'road_length' should be in km
- Use ggplot to plot districts with black outline (color), and empty fill (fill = NA)
- Then plot roads based on length (use `scale_color_continuous()` )
- Plot roads with line width 1 (`lwd = 1`)
- Use `theme_bw()`

---

```{r, eval=FALSE}
Expand All @@ -132,181 +141,3 @@ ggplot() +
- Add aesthetics in `aes()`. These are columns used for plotting.
- For `geom_point`, `geom_line` you need 'x' and 'y' aesthetics.
- For `geom_sf`, spatial context used for 'x' and 'y'. But you might want 'fill' or 'color'

---
## Calculate length (exercise)
- Use `st_length()` to add column 'road_length' to `roads`
- 'road_length' should be in km
- Use ggplot to plot districts with black outline (color), and empty fill (fill = NA)
- Then plot roads based on length (use `scale_color_continuous()` )
- Plot roads with line width 1 (`lwd = 1`)
- Use `theme_bw()`

---

## Spatial Operations
Data setup
```{r, eval = FALSE, message = FALSE}
coords <- cbind("x" = c(25, 25.1, 26.4, 26.4, 26.4, 25),
"y" = c(-15.5, -14.5, -14.6, -14.8, -15.6, -15.5))
pol <- st_polygon(x = list(coords)) %>%
st_sfc %>%
st_sf(ID = 1, crs = 4326)
roads_rpj <- st_transform(roads, crs = st_crs(districts))
```

---
## Plotting
```{r, eval = FALSE, message = FALSE}
par(mar = rep(0, 4))
plot(st_geometry(districts), col = "grey")
plot(pol, col = "purple", add = TRUE)
plot(roads_rpj, col = 'red', add = TRUE)
```


---
## `st_intersects`
- Similar to "Select by Spatial Location"
- Finds indices of second object that intersect each element of first object
```{r, eval = FALSE, message = FALSE}
st_intersects(x = roads_rpj, y = districts)
```

---
## `st_intersects`
```{r, eval = FALSE, message = FALSE}
highway <- roads_rpj[401, ]
dists_int <- districts %>%
slice(st_intersects(x = highway, y = districts)[[1]])
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
dists_int
plot(st_geometry(dists_int), col = 'grey')
plot(highway, col = 'red', add = TRUE)
```


---
## `st_intersection()`
- finds intersecting pairs between two objects
- Each row in result is one "intersection object" (i.e. an object from x and y that intersect)
```{r, eval = FALSE, message = FALSE}
pol_int_dists <- st_intersection(x = pol, y = districts)
par(mar = rep(0, 4))
plot(st_geometry(districts), col = "grey")
plot(st_geometry(pol_int_dists), col = rainbow(n = nrow(pol_int_dists)),
add = TRUE)
```

---
## `st_difference`
- Similar to "erase" in ArcGIS
- Returns first object without second object.
```{r, eval = FALSE, message = FALSE}
d1 <- st_difference(x = districts, y = pol)
d2 <- st_difference(x = pol, y = districts)
par(mfrow = c(1, 2), mar = c(0, 0, 1, 0))
d1 %>%
st_geometry %>%
plot(col = "grey")
d2 %>%
st_geometry %>%
plot(col = "grey")
```


---
## `st_union`
- Combines objects without resolving borders.
- Really takes all possible combinations of first and 2nd object
- Switching order of inputs changes column order
```{r, eval = FALSE, message = FALSE}
uni1 <- st_union(x = pol, y = districts)
uni1
uni2 <- st_union(x = districts, y = pol)
uni2
par(mfrow = c(1, 2), mar = c(0, 0, 1, 0))
plot(uni1[2], reset = FALSE, main = "Union of pol with districts")
plot(uni2[1], reset = FALSE, main = "Union of districts with pol")
```

---
## `st_buffer()`
- Best to use projected coordinates (e.g. Albers)
```{r, eval = FALSE, message = FALSE}
## project everything to Albers
farmers_alb <- farmers_sf %>% st_transform(st_crs(roads))
long_roads <- roads %>% filter(as.numeric(st_length(.)) > 500000)
districts_alb <- districts %>% st_transform(st_crs(roads))
roads_buff_20km = st_buffer(long_roads, 20000)
plot(st_geometry(districts_alb), col = "grey")
plot(st_geometry(roads_buff_20km), col = "blue", add = TRUE)
plot(st_geometry(long_roads), col = "red", add = TRUE)
```

---
## `st_sample`

```{r, eval = FALSE, message = FALSE}
par(mar = rep(0, 4))
districts %>% st_union %>% plot(col = "grey")
set.seed(1)
districts %>%
st_union %>%
st_sample(size = 100, exact = TRUE) %>%
plot(pch = 20, add = TRUE)
```


---
## Exercise
Find all districts within 50 km of sampled points.

```{r, eval = FALSE, message = FALSE}
```

---
## Exercise
Of the sampled farmers, how many are within 100km of one of the "long roads"?
```{r, eval = FALSE, message = FALSE}
```


---
## Further Exercises

### Round 1
- Calculate the length of roads in kilometers
- Buffer roads > 100 km by 30 kilometers, save as new object
- Buffer roads > 100 km 10 kilometers, save as new object
- Take the difference between the 30 km and 10 km buffers, i.e. erase/difference the two
- Plot the resulting difference using `ggplot` with the difference in red over Zambia's districts

```{r, eval=FALSE}
# code to adapt
roads %>% mutate(km = as.numeric(st_length(.)) / 1000) %>%
filter(km > 200) %>%
st_buffer(dist = 20000) %>%
st_transform(crs = 4326)
st_difference(larger_object, smaller_object)
ggplot() + geom_sf
```

---

### Round 2
- Randomly sample 100 points within the 30 km road buffer (use a seed of 1)
- Plot the 30 km buffer (red) the points in it (blue) over the districts

```{r, eval=FALSE}
buffered_object %>% st_sample()
```
Loading

0 comments on commit 8f2814a

Please sign in to comment.