-
-
Notifications
You must be signed in to change notification settings - Fork 594
/
Copy path_04-ex.Rmd
152 lines (123 loc) · 5.79 KB
/
_04-ex.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
```{r 04-ex-e0, include=TRUE}
library(sf)
library(dplyr)
data(nz, package = "spData")
data(nz_height, package = "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.
How many of these high points does the Canterbury region contain?
```{r 04-ex-e1}
library(tmap)
tmap_mode("view")
qtm(nz) + qtm(nz_height)
canterbury = nz %>% filter(Name == "Canterbury")
canterbury_height = nz_height[canterbury, ]
nrow(canterbury_height) # answer: 70
```
E2. Which region has the second highest number of `nz_height` points in, and how many does it have?
```{r 04-ex-e2}
nz_height_count = aggregate(nz_height, nz, length)
nz_height_combined = cbind(nz, count = nz_height_count$elevation)
nz_height_combined %>%
st_drop_geometry() %>%
dplyr::select(Name, count) %>%
arrange(desc(count)) %>%
slice(2)
```
E3. Generalizing the question to all regions: how many of New Zealand's 16 regions contain points which belong to the top 100 highest points in the country? Which regions?
- Bonus: create a table listing these regions in order of the number of points and their name.
```{r 04-ex-e3}
nz_height_count = aggregate(nz_height, nz, length)
nz_height_combined = cbind(nz, count = nz_height_count$elevation)
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).
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}
library(terra)
dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))
ndvi = rast(system.file("raster/ndvi.tif", package = "spDataLarge"))
#1
dem_rcl = matrix(c(-Inf, 300, 0, 300, 500, 1, 500, Inf, 2), ncol = 3, byrow = TRUE)
dem_reclass = classify(dem, dem_rcl)
levels(dem_reclass) = c("low", "medium", "high")
plot(dem_reclass)
#2
zonal(c(dem, ndvi), dem_reclass, fun = "mean")
```
E5. 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}
# 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):
# fx=matrix(c(-1,-2,-1,0,0,0,1,2,1), nrow=3)
# fy=matrix(c(1,0,-1,2,0,-2,1,0,-1), nrow=3)
# just retrieve the first channel of the R logo
r = rast(system.file("ex/logo.tif", package = "terra"))
# compute the Sobel filter
filter_x = matrix(c(-1, -2, -1, 0, 0, 0, 1, 2, 1), nrow = 3)
sobel_x = focal(r, w = filter_x)
plot(sobel_x, col = c("white", "black"))
filter_y = matrix(c(1, 0, -1, 2, 0, -2, 1, 0, -1), nrow = 3)
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.
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}
file = system.file("raster/landsat.tif", package = "spDataLarge")
multi_rast = rast(file)
ndvi_fun = function(nir, red){
(nir - red) / (nir + red)
}
ndvi_rast = lapp(multi_rast[[c(4, 3)]], fun = ndvi_fun)
plot(ndvi_rast)
ndwi_fun = function(green, nir){
(green - nir) / (green + nir)
}
ndwi_rast = lapp(multi_rast[[c(2, 4)]], fun = ndwi_fun)
plot(ndwi_rast)
two_rasts = c(ndvi_rast, ndwi_rast)
names(two_rasts) = c("ndvi", "ndwi")
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()`.
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}
# Fetch the DEM data for Spain
spain_dem = geodata::elevation_30s(country = "Spain", path = ".", mask = FALSE)
# Reduce the resolution by a factor of 20 to speed up calculations
spain_dem = aggregate(spain_dem, fact = 20)
# According to the documentation, terra::distance() will calculate distance
# for all cells that are NA to the nearest cell that are not NA. To calculate
# distance to the coast, we need a raster that has NA values over land and any
# other value over water
water_mask = is.na(spain_dem)
water_mask[water_mask == 0] = NA
# Use the distance() function on this mask to get distance to the coast
distance_to_coast = distance(water_mask)
# convert distance into km
distance_to_coast_km = distance_to_coast / 1000
# Plot the result
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.
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}
# 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
plot(distance_to_coast_km2)
# visualize the difference
plot(distance_to_coast_km - distance_to_coast_km2)
```