Update answer to c4 exercises

_04-ex.Rmd

@@ -1,8 +1,7 @@
 ```{r 04-ex-e0, include=TRUE, message=FALSE}
 library(sf)
 library(dplyr)
-data(nz, package = "spData")
-data(nz_height, package = "spData")
+library(spData)
 ```
 
 E1. It was established in Section \@ref(spatial-vec) that Canterbury was the region of New Zealand containing most of the 100 highest points in the country.
@@ -43,19 +42,93 @@ E3. Generalizing the question to all regions: how many of New Zealand's 16 regio
 - Bonus: create a table listing these regions in order of the number of points and their name.
 
 ```{r 04-ex-e3}
+# Base R way:
 nz_height_count = aggregate(nz_height, nz, length)
 nz_height_combined = cbind(nz, count = nz_height_count$elevation)
+plot(nz_height_combined)
+
+# Tidyverse way:
+nz_height_joined = st_join(nz_height, nz %>% select(Name))
+# Calculate n. points in each region - this contains the result
+nz_height_counts = nz_height_joined %>% 
+  group_by(Name) %>% 
+  summarise(count = n())
+
+# Optionally join results with nz geometries:
+nz_height_combined = left_join(nz, nz_height_counts %>% sf::st_drop_geometry())
+# plot(nz_height_combined) # Check: results identical to base R result
+
+# Generate a summary table
 nz_height_combined %>% 
   st_drop_geometry() %>% 
   dplyr::select(Name, count) %>% 
   arrange(desc(count)) %>% 
   na.omit()
 ```
 
-E4. Use `dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))`, and reclassify the elevation in three classes: low (<300), medium and high (>500).
+E4. Test your knowledge of spatial predicates by finding out and plotting how US states relate to each other and other spatial objects.
+
+The starting point of this exercise is to create an object representing Colorado state in the USA. Do this with the command 
+`colorado = us_states[us_states$NAME == "Colorado",]` (base R) or with with the  `filter()` function (tidyverse) and plot the resulting object in the context of US states.
+
+- Create a new object representing all the states that geographically intersect with Colorado and plot the result (hint: the most concise way to do this is with the subsetting method `[`).
+- Create another object representing all the objects that touch (have a shared boundary with) Colorado and plot the result (hint: remember you can use the argument `op = st_intersects` and other spatial relations during spatial subsetting operations in base R).
+- Bonus: create a straight line from the centroid of the District of Columbia near the East coast to the centroid of California near the West coast of the USA (hint: functions `st_centroid()`, `st_union()` and `st_cast()` described in Chapter 5 may help) and identify which states this long East-West line crosses.
+
+```{r 04-ex-4-1}
+colorado = us_states[us_states$NAME == "Colorado", ]
+plot(us_states$geometry)
+plot(colorado$geometry, col = "grey", add = TRUE)
+```
+
+```{r 04-ex-4-2}
+intersects_with_colorado = us_states[colorado, , op = st_intersects]
+plot(us_states$geometry, main = "States that intersect with Colorado")
+plot(intersects_with_colorado$geometry, col = "grey", add = TRUE)
+```
+
+```{r 04-ex-4-3}
+# Alternative but more verbose solutions
+# 2: With intermediate object, one list for each state
+sel_intersects_colorado = st_intersects(us_states, colorado)
+sel_intersects_colorado_list = lengths(sel_intersects_colorado) > 0
+intersects_with_colorado = us_states[sel_intersects_colorado_list, ]
+
+# 3: With intermediate object, one index for each state
+sel_intersects_colorado2 = st_intersects(colorado, us_states)
+sel_intersects_colorado2
+us_states$NAME[unlist(sel_intersects_colorado2)]
+
+# 4: With tidyverse
+us_states %>% 
+  st_filter(y = colorado, .predicate = st_intersects)
+```
+
+```{r 04-ex-4-4}
+touches_colorado = us_states[colorado, , op = st_touches]
+plot(us_states$geometry, main = "States that touch Colorado")
+plot(touches_colorado$geometry, col = "grey", add = TRUE)
+```
+
+
+```{r 04-ex-4-5}
+washington_to_cali = us_states %>% 
+  filter(grepl(pattern = "Columbia|Cali", x = NAME)) %>% 
+  st_centroid() %>% 
+  st_union() %>% 
+  st_cast("LINESTRING")
+states_crossed = us_states[washington_to_cali, , op = st_crosses]
+states_crossed$NAME
+plot(us_states$geometry, main = "States crossed by a straight line\n from the District of Columbia to central California")
+plot(states_crossed$geometry, col = "grey", add = TRUE)
+plot(washington_to_cali, add = TRUE)
+```
+
+
+E5. Use `dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))`, and reclassify the elevation in three classes: low (<300), medium and high (>500).
 Secondly, read the NDVI raster (`ndvi = rast(system.file("raster/ndvi.tif", package = "spDataLarge"))`) and compute the mean NDVI and the mean elevation for each altitudinal class.
 
-```{r 04-ex-e4}
+```{r 04-ex-e5}
 library(terra)
 dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))
 ndvi = rast(system.file("raster/ndvi.tif", package = "spDataLarge"))
@@ -70,11 +143,11 @@ plot(dem_reclass)
 zonal(c(dem, ndvi), dem_reclass, fun = "mean")
 ```
 
-E5. Apply a line detection filter to `rast(system.file("ex/logo.tif", package = "terra"))`.
+E6. Apply a line detection filter to `rast(system.file("ex/logo.tif", package = "terra"))`.
 Plot the result.
 Hint: Read `?terra::focal()`.
 
-```{r 04-ex-e5}
+```{r 04-ex-e6}
 # from the focal help page (?terra::focal()):
 # Laplacian filter: filter=matrix(c(0,1,0,1,-4,1,0,1,0), nrow=3)
 # Sobel filters (for edge detection): 
@@ -93,11 +166,11 @@ sobel_y = focal(r, w = filter_y)
 plot(sobel_y, col = c("black", "white"))
 ```
 
-E6. Calculate the Normalized Difference Water Index	(NDWI; `(green - nir)/(green + nir)`) of a Landsat image. 
+E7. Calculate the Normalized Difference Water Index	(NDWI; `(green - nir)/(green + nir)`) of a Landsat image. 
 Use the Landsat image provided by the **spDataLarge** package (`system.file("raster/landsat.tif", package = "spDataLarge")`).
 Also, calculate a correlation between NDVI and NDWI for this area.
 
-```{r 04-ex-e6}
+```{r 04-ex-e7}
 file = system.file("raster/landsat.tif", package = "spDataLarge")
 multi_rast = rast(file)
 
@@ -120,12 +193,12 @@ two_rasts_df = as.data.frame(two_rasts)
 cor(two_rasts_df$ndvi, two_rasts_df$ndwi)
 ```
 
-E7. A StackOverflow [post](https://stackoverflow.com/questions/35555709/global-raster-of-geographic-distances) shows how to compute distances to the nearest coastline using `raster::distance()`.
+E8. A StackOverflow [post](https://stackoverflow.com/questions/35555709/global-raster-of-geographic-distances) shows how to compute distances to the nearest coastline using `raster::distance()`.
 Try to do something similar but with `terra::distance()`: retrieve a digital elevation model of Spain, and compute a raster which represents distances to the coast across the country (hint: use `geodata::elevation_30s()`).
 Convert the resulting distances from meters to kilometers.
 Note: it may be wise to increase the cell size of the input raster to reduce compute time during this operation.
 
-```{r 04-ex-e7}
+```{r 04-ex-e8}
 # Fetch the DEM data for Spain
 spain_dem = geodata::elevation_30s(country = "Spain", path = ".", mask = FALSE)
 
@@ -148,10 +221,10 @@ distance_to_coast_km = distance_to_coast / 1000
 plot(distance_to_coast_km, main = "Distance to the coast (km)")
 ```
 
-E8. Try to modify the approach used in the above exercise by weighting the distance raster with the elevation raster; every 100 altitudinal meters should increase the distance to the coast by 10 km.
+E9. Try to modify the approach used in the above exercise by weighting the distance raster with the elevation raster; every 100 altitudinal meters should increase the distance to the coast by 10 km.
 Next, compute and visualize the difference between the raster created using the Euclidean distance (E7) and the raster weighted by elevation.
 
-```{r 04-ex-e8}
+```{r 04-ex-e9}
 # now let's weight each 100 altitudinal meters by an additional distance of 10 km
 distance_to_coast_km2 = distance_to_coast_km + ((spain_dem / 100) * 10)
 # plot the result