-## Spring 2023
+## Spring 2024
Updated on: `r Sys.Date()`
Unit 1 Practice answers
-GEOG246-346
- - - - - - - - - -1 Module 1 practice answers
-N/A
-2 Module 2 practice answers
-N/A
-3 Module 3 practice answers
-3.1 Practice 1
-3.1.1 Questions
--
-
Class of
ais integer.
-Because the vector is contained in a list, and using a single
[]to pull out an element returns a 1-element list, so you have to use list notation to get the vector.
-ais vector of integers,lis a list of containing 2 integer vectors and 1 character vector.
-We applied a function that is the fourth element of list
l2to the integer vectorfthat is the second element ofl2.
-
3.1.2 Code
--
-
-
a <- 20:30
-b <- letters-
-
-
names(a) <- b[1:length(a)]
-a-
-
-
l <- list(a = a, b = b)-
-
-
a[a >= 26]
-a[c(1, 7)]
-a[c(length(a) - 1, length(a))]-
-
-
b[b %in% c("a", "c", "g")]-
-
-
l[[1]]
-l[1]
-l$a[l$a < 25]
-l[["a"]][l[["a"]] == 25]
-l[[2]][l[[2]] %in% c("d", "e", "f")]3.2 Practice 2
-3.2.1 Questions
--
-
The vector is coerced to character.
-The whole
matrixis coerced to character. Ifdata.frame, onlym$bis coerced to character.
--
-
- Have to use both column and row indices to attract row, col subset; -
- Can use list notation for
data.frame(e.g.m$b), but not formatrix.
-
-
3.2.2 Code
--
-
-
m <- cbind(1:10, 11:20, 21:30)-
-
-
m[4:5, 2:3]-
-
-
colnames(m) <- letters[1:3]-
-
-
m[m[, "b"] > 14 & m[, "b"] <= 18, ]-
-
-
d <- as.data.frame(m)-
-
-
d[d$a > 4, "a"] <- -1 -
-
-
d$c <- letters[1:10]-
-
-
l <- list(m, d)
-l[[2]][2:3, ]-
-
-
d %>% filter(b >= 14 & b <= 18)3.3 Practice 3
-3.3.1 Questions
--
-
1 and 5 pulled from
band multiplied in sequence onmby row then column.
-Would multiply 1st columns of each, 2nd columns of each, and 3rd columns of each.
-drow 1 columns 1 to 3 multiplied by the value from the last row, first column ofm.
-
3.3.2 Code
--
-
Not shown here due to length
-
-
sin(m)
-cos(m)-
-
-
rowSums(d[, 1:3])
-colSums(d[, 1:3])
-rowMeans(d[, 1:3])
-colMeans(d[, 1:3])3.4 Practice 4
-3.4.1 Code
--
-
-
It adds an NA to print statement for the 5th letter of the alphabet. Expand sscript by adding another “th” to the vector, so it has 5 elements.
sscript <- c("st", "nd", "rd", "th", "th") # vector of superscripts
-for(i in 1:5) { # for loop with iterator i over vector 1:4
- stmnt <- paste0(letters[i], " is the ", i, sscript[i],
- " letter in the alphabet")
- print(stmnt) # print statement
-}-
-
-
Turns 4 of the letters red and adds axes to plots
-sscript <- c("st", "nd", "rd", "th") # vector of superscripts
-par(mfrow = c(1, 4), mar = c(0, 0, 1, 0.5))
-for(i in 1:4) {
- stmnt <- paste0(letters[i], " is the ", i, sscript[i],
- " letter in the alphabet")
- plot(1:4, rep(3, 4), ylim = c(1, 5), pch = letters[1:4], #axes = FALSE,
- xlab = "", ylab = "", main = stmnt, cex = 2)
- # points(i, 3, pch = letters[i], col = "red", cex = 2)
-}-
-
-
# Chunk 26
-for(i in 1:10) {
- if(i < 5) { # condition 1
- print(paste(i, "is less than", i + 1))
- } else if(i >= 3 & i <= 8) { # condition 2
- print(paste(i, "is between", i - 1, "and", i + 1))
- } else { # remaining conditions
- print(paste(i, "is greater than", i - 1))
- }
-}-
-
-
Create a for loop that iterates over a vector 1:20. Insert a condition into it such that it only prints out a result when the iterator’s value is 11.
for(i in 1:20) {
- if(i == 11) print(i)
-}3.5 Practice 5
-3.5.1 Code
--
-
-
Turns the summary of “a” variables to NA.
-dat_list <- list(data.frame(a = 1:10, b = 21:30),
- data.frame(a = 31:40, b = 41:50),
- data.frame(a = 51:60, b = 61:70))
-dat_modify <- function(x) {
- x[1:3, 1] <- 999
- return(x)
-}
-dat_list <- lapply(dat_list, dat_modify)
-dat_list
-
-dat_list2 <- c(dat_list, mean) # add another element to dat_list
-lapply(1:length(dat_list2), function(x) { # x <- 1
- d <- dat_list2[[x]] # extract element of list
- if(is.data.frame(d)) { # check if it is a data.frame
- d[d == 999] <- NA # convert any 999 values to NA
- o <- c(colSums(d, na.rm = FALSE), # column sums, dropping NAs
- "total" = sum(d, na.rm = TRUE)) # sum dropping NAs
- } else { # if it is not a data.frame, make an error statement
- o <- paste("Operation not valid for a", class(d))
- }
- return(o) # return result
-})-
-
-
Can’t coerce list to double.
-flist <- list(mean, sd, range)
-lapply(1:3, function(x) flist[[x]](dat_list[[1]]))-
-
-
lapply(dat_list, function(x) mean(unlist(x)))
-lapply(dat_list, function(x) sd(unlist(x)))-
-
-
sapply(dat_list, function(x) mean(unlist(x)))
-sapply(dat_list, function(x) sd(unlist(x)))-
-
-
lapply(dat_list, function(x) sum(x[1:2, 2]))4 Module 4 practice answers
-4.1 Practice 1
-4.1.1 Questions
--
-
A
tibbleis an enhanceddata.frame. Among other features, it provides more information on data types it contains when printing. It can be treated exactly like adata.framethough, in terms of indexing and other operations.
-Base methods:
tb_a$a;tb_a[["a"]];tb_a[, "a"]. tidyverse, withtb_a %>% dplyr::select(a). Note that that gives back atibble(ordata.frame, if you have adata.frame). To get a vector back, usetb_a %>% dplyr::pull(a)
-The data are messy, because the “column headers are values, not variable names” (see here). We would gather the data in the month rows, setting the month as the key and the mean flu cases as the value.
-inner_joinpreserves just the values in x and drops non-matching rows from y.left_joinfills non-matching values in y with NAs.right_joinpreserves values in y dropping non-matching values in x.full_joinpreserves all non-matching rows in both x and y, filling NAs into non-matched rows.
-
4.1.2 Code
--
-
-
set.seed(1)
-dat <- tibble(a = sample(1:10), b = rnorm(10))
-td <- "/path/where/you/want/to/write" # REPLACE THIS WITH YOUR OWN!!!!
-readr::write_csv(dat, path = file.path(td, "dummy.csv"))-
-
-
td <- "/path/where/you/want/to/write" # REPLACE THIS WITH YOUR OWN!!!!
-readr::read_csv(file.path(td, "dummy.csv"))-
-
-
Recreate tibbles first
# Chunk 13
-set.seed(1)
-t1 <- tibble(v1 = paste0("N", 1:5), v2 = rnorm(5))
-t2 <- tibble(v1 = paste0("N", 1:5), v3 = runif(5))
-t3 <- tibble(v1 = paste0("N", 1:7), v4 = sample(1:100, 7))
- # v5 = letters[sample(1:26, 7)])
-t4 <- tibble(v1 = paste0("N", c(1:2, 4:7, 11)),
- v5 = letters[sample(1:26, 7)])Then do joins:
-left_join(t1, t2) %>% left_join(., t3) %>% left_join(., t4)
-right_join(t1, t2) %>% right_join(., t3) %>% right_join(., t4)-
-
-
left_join(t1, t2) %>% left_join(., t3) %>% left_join(., t4) %>% arrange(v5)
-right_join(t1, t2) %>% right_join(., t3) %>% right_join(., t4) %>%
- arrange(desc(v5))-
-
-
fs <- dir(system.file("extdata/", package = "geospaar"), pattern = "FAOSTAT",
- full.names = TRUE)
-crops <- lapply(fs, readr::read_csv)
-crops_df <- do.call(rbind, lapply(crops, function(x) {
- x %>% dplyr::select(Item, Area, Element, Year, Value) %>%
- pivot_wider(names_from = Element, values_from = Value) %>%
- rename(crop = Item, country = Area, year = Year,
- harv_area = `Area harvested`, prod = Production)
-}))
-crop_ylds <- crops_df %>% mutate(yield = prod / harv_area)
-crop_ylds <- crop_ylds %>%
- mutate(country = ifelse(country == "South Africa", "ZAF", country)) %>%
- mutate(country = ifelse(country == "Zambia", "ZMB", country)) %>%
- mutate(harv_km2 = harv_area / 100)-
-
-
crop_ylds %>% rename(harv_area_km2 = harv_km2)-
-
-
my_tb <- tibble(v1 = 1:10, v2 = 11:20) %>%
- rbind(., tibble(v1 = 11:20, v2 = 21:30)) %>% mutate(v3 = v2^2) %>%
- arrange(-v3)-
-
-
my_tb %>% slice(1, 10, 17) %>% dplyr::select(v2, v3)4.2 Practice 2: Analysis
-4.2.1 Questions
--
-
dplyr::filter
-group_by(crop, country, y2k)is doing the splitting, on crop type, then country, and then year.summarize(...)is doing the apply using amean. There is no combine line, as it is implicit.
-Chunk 30 adds a
filterfor crop type (selecting out maize), and then simply groups on the y2k variable.
-It doesn’t work when the output of the analysis is not tabular/a list, as with
cor.testandlm. We can overcome this by 1) creating individual functions that reproduce the component outputs of the analysis (e.g. Chunk 33) and these as a list of functions usingfunstosummarise_all, 2) doing the splits outside of the pipeline (e.g. Chunk 36), or 3) using functions such asdoandbroom::tidywithin the pipeline (e.g. Chunk 39).
-
4.2.2 Code
--
-
-
crop_ylds %>% filter(crop == "sorghum" & country == "ZAF" & year >= 2000)
-crop_ylds[crop_ylds$crop == "sorghum" & crop_ylds$country == "ZAF" &
- crop_ylds$year >= 2000, ]-
-
-
crop_ylds %>% filter(crop == "sorghum" & country == "ZAF" & year >= 2000) %>%
- select(prod, harv_area, yield) %>% summarise_all(funs(mean, sd))-
-
-
crop_ylds %>% group_by(crop, country) %>% select(prod, harv_area) %>%
- summarise_all(funs(mean, sd))-
-
-
crop_ylds %>% filter(crop == "maize" & country == "ZMB") %>%
- select(yield, harv_area) %>% cor()
-
-dat <- crop_ylds %>% filter(crop == "maize" & country == "ZMB")
-cor.test(dat$harv_area, dat$yield)-
-
-
crop_ylds %>% filter(crop == "maize") %>% summarize(mu_yld = mean(yield))
-# 2.07 t/ha-
-
-
South Africa shows larger yields gains (0.061 t/ha/yr versus 0.03 t/ha/yr)
-summary(lm(yield ~ year,
- data = crop_ylds %>% filter(crop == "maize" & country == "ZMB")))
-summary(lm(yield ~ year,
- data = crop_ylds %>% filter(crop == "maize" & country == "ZAF")))
-
-crop_ylds %>% filter(crop == "maize" & country == "ZMB") %>%
- lm(yield ~ year, data = .) %>% summary()
-crop_ylds %>% filter(crop == "maize" & country == "ZAF") %>%
- lm(yield ~ year, data = .) %>% summary()-
-
-
crop_ylds %>% filter(crop != "sorghum") %>% group_by(crop, country) %>%
- do(prod_ha_lm = lm(yield ~ year, data = .)) %>%
- broom::tidy(., prod_ha_lm)4.3 Practice 3: Visualization
-4.3.1 Questions
--
-
ggplot2is built fromgridgraphics, and is based on an underlying visualization philosophy. It builds up graphics objects using the+operator, and easily does splits within the data using the “color” argument withinaesand/or thefacet_gridfunction.graphicsplots can be faster to implement for exploratory analysis,ggplot2has more attractive, presentation-grade defaults.
-Because the syntax used in
graphicsplots is used in many of the plotting functions developed for spatial packages, including newer ones such assfandstars.
-Each of the three plots takes the axis labels exactly as they are specified to the axis arguments (either as they were specified in the formula in the case of Chunk 40 or to the “x” and “y” arguments in Chunk 41). You can change the names using the “xlab” and “ylab” arguments.
-Using “col”, “pch”, and “cex” arguments.
-You have to add the
.to the “data” argument ofplot, e.g.dat %>% plot(y ~ x, data = .)
-You get just an empty grey background–
ggplotwon’t plot anything without ageom_*function added to theggplotobject.
-
4.3.2 Code
--
-
-
crop_ylds %>% filter(crop == "sorghum") %>%
- ggplot() + geom_histogram(aes(x = yield), bins = 15) +
- ggtitle("Distribution of sorghum yields")
-crop_ylds %>% filter(crop == "sorghum") %>%
- ggplot() + geom_histogram(aes(x = yield), bins = 15, fill = "red") +
- ggtitle("Distribution of sorghum yields")-
-
-
crop_ylds %>% filter(crop == "wheat" & country == "ZAF") %>%
- plot(harv_area ~ year, data = ., pch = 16, col = "blue",
- xlab = "", ylab = "Harvested area (ha)",
- main = "South Africa wheat (1961-2017)")
-crop_ylds %>% filter(crop == "wheat" & country == "ZAF") %>%
- plot(harv_area ~ year, data = ., pch = 16, type = "l", col = "blue",
- xlab = "", ylab = "Harvested area (ha)",
- main = "South Africa wheat (1961-2017)")-
-
-
crop_ylds %>% filter(crop == "wheat" & country == "ZAF") %>%
- ggplot() + geom_point(aes(year, harv_area), col = "blue") +
- xlab("") + ylab("Harvested area (ha)") +
- ggtitle("South Africa wheat (1961-2017)")
-crop_ylds %>% filter(crop == "wheat" & country == "ZAF") %>%
- ggplot() + geom_line(aes(year, harv_area), col = "blue") +
- xlab("") + ylab("Harvested area (ha)") +
- ggtitle("South Africa wheat (1961-2017)")-
-
-
crop_ylds %>% filter(crop == "wheat") %>%
- ggplot() + geom_line(aes(year, harv_area, color = country)) +
- scale_color_manual(values = c("red", "blue")) +
- xlab("") + ylab("Harvested area (ha)") +
- ggtitle("Wheat (1961-2017)")
-
-# extra
-crop_ylds %>% filter(crop == "wheat") %>%
- ggplot() + geom_line(aes(year, log10(harv_area), color = country)) +
- scale_color_manual(values = c("red", "blue")) +
- xlab("") + ylab("Harvested area (ha)") +
- ggtitle("Wheat (1961-2017)")-
-
-
crop_ylds %>% filter(crop == "wheat" & country == "ZAF") %>%
- ggplot() + geom_point(aes(year, harv_area)) +
- geom_smooth(aes(year, harv_area)) +
- xlab("") + ylab("Harvested area (ha)") +
- ggtitle("South African wheat (1961-2017)")-
-
-
# ggplot2
-crop_ylds %>% filter(crop == "wheat" & country == "ZMB") %>%
- ggplot() +
- geom_histogram(aes(x = harv_area), bins = 10, col = "black", fill = "blue") +
- xlab("Harvested area (ha)") + ggtitle("Zambian Wheat (1961-2017)")
-
-# hist
-# with dplyr
-crop_ylds %>% filter(crop == "wheat" & country == "ZMB") %>% pull(harv_area) %>%
- hist(., main = "Zambian Wheat (1961-2017)", xlab = "Harvested area (ha)",
- col = "blue")
-# with base subsetting
-hist(crop_ylds$harv_area[crop_ylds$crop == "wheat" &
- crop_ylds$country == "ZMB"],
- main = "Zambian Wheat (1961-2017)", xlab = "Harvested area (ha)",
- col = "blue")-
-
-
crop_ylds %>% filter(country == "ZAF") %>%
- ggplot() + geom_point(aes(x = year, y = harv_area)) +
- geom_smooth(aes(x = year, y = harv_area)) +
- facet_grid(cols = vars(crop)) +
- scale_color_manual(values = c("red", "blue")) +
- ylab("Yield (tons/ha)") + xlab("")-
-
Unit 2 Practice answers
-GEOG246-346
- - - - - - - - - -1 Module 1 practice answers
-1.1 Practice 1
-1.1.1 Questions
--
-
Answers here.
-Assuming the
tibblehas x and y or lat/long coordinates, you apply the functionst_as_sfwith the “coords” argument set to specify which columns contain the x and y coordinates.
-sf::plotby default plots one panel per variable. You can create a single panel by specifying the variable you want, or by using thest_geometryargument to strip out the geometry from object. It also prevent distortions that sometimes occur when overlaying subsequent features on top of a base map.
-For a single point you provide an x and y coordinate, otherwise you give an input matrix containing x and y coordinates. A polygon requires that the last point pair in the matrix is the same as the first point pair, to close the polygon.
-
1.1.2 Code
--
-
-
farmers %>% st_geometry()
-st_geometry(farmers)-
-
-
st_crs(farmers) <- st_crs(districts)
-# p <- "path/to/your/project/notebooks/data"
-p <- "~/Desktop"
-st_write(farmers, dsn = file.path(p, "farmers.sqlite"))
-rm(farmers)
-st_read(file.path(p, "farmers.sqlite"))-
-
-
class(roads)
-str(roads)
-class(districts)
-str(districts)-
-
-
plot(roads %>% st_geometry(), col = "blue")-
-
-
plot(districts %>% select(distName), main = "Zambia Districts")-
-
-
pts <- st_multipoint(x = cbind(x = c(27, 28, 29), y = c(-13, -14, -15)))
-plot(districts %>% st_geometry(), col = "grey")
-plot(pts, add = TRUE, col = "orange", pch = 16)1.2 Practice 2
-1.2.1 Questions
--
-
At least with this example, pretty negligible–well less than 1% mean absolute error. It might matter more in other places and with other scales though. The reason
st_areaknows how to estimate areas is because it invokeslwgeom::st_geod_area, which calculates a geodetic surface area.
-Because for the time being
sf::plotis a fair bit faster. However, this recent twitter thread suggests that may change soon.
-By using
mutatewithcutthat has break values based on those properties. In our example, we found the breaks usingquantileand different probabilities/percentile levels that creating tertiles of area.
-
1.2.2 Code
--
-
-
set.seed(1)
-districts %>% sample_n(20) %>% st_area() %>% units::set_units("ha") %>% mean()-
-
-
set.seed(1)
-roads %>% sample_n(100) %>% st_length() %>% units::set_units("km") %>% mean()-
-
-
plot(st_geometry(districts), col = "lightgrey")
-set.seed(1)
-farmers %>% filter(season == 2) %>% sample_n(200) %>% st_geometry() %>%
- plot(col = "red", pch = 20, add = TRUE)-
-
-
districts %>% st_transform(st_crs(roads)) %>% st_geometry() %>%
- plot(col = "lightgrey")
-roads %>%
- mutate(length = as.numeric(st_length(.) / 1000)) %>%
- filter(length > 50 & length < 100) %>% st_geometry() %>%
- plot(col = "red", pch = 20, add = TRUE)-
-
-
deciles <- function(x) quantile(x, probs = seq(0, 1, 0.1))
-dist_deciles <- districts %>% mutate(area = as.numeric(st_area(.)) / 10^6) %>%
- mutate(acls = cut(area, breaks = deciles(area), include.lowest = TRUE)) %>%
- group_by(acls) %>% summarize(sum_area = sum(area))
-dist_deciles
-#
-# #2
-cols <- heat.colors(10)
-par(mar = rep(0, 4))
-plot(st_geometry(dist_deciles), col = cols)
-legend(x = "bottomright", pch = 15, col = cols, bty = "n",
- legend = paste0(1:10, c("st", "nd", "rd", rep("th", 7))))1.2.3 Practice 3
-1.2.3.1 Questions
--
-
It changes features from one type to another (e.g. POLYGON to MULTIPOLYGON), either one specified by the user or the simplest possible common feature, if left unspecified. Casting is sometimes necessary to avoid mixed feature types that cause failures for subsequent operations.
-st_unionruns under the hood ofsummarise.sf, so asummarizeoperation on ansfwill result in a merged/dissolved set of spatial features.
-It affects the order of the fields in the resulting unioned
sfobject–the fields from the object passed to the “x” argument appear first.
-
1.2.3.2 Code
--
-
-
coords <- cbind("x" = c(27, 27.5, 27.5, 27, 27),
- "y" = c(-13, -13, -13.5, -13.5, -13))
-pol2 <- st_polygon(x = list(coords)) %>% st_sfc %>% st_sf(ID = 1, crs = 4326)
-
-par(mar = rep(0, 4))
-plot(st_geometry(districts), col = "grey")
-plot(pol2, col = "blue", add = TRUE)-
-
-
pol2_int_dists <- st_intersection(pol2, districts)
-districts[1] %>% plot(col = "grey", main = NULL)
-pol2_int_dists[2] %>% plot(col = rainbow(nrow(pol2_int_dists)), add = TRUE)-
-
-
pol2 is nearly 26 hectares larger
pol2_int_dists %>% st_area() %>% as.numeric() %>% sum() / 10000 -
- pol2 %>% st_area() %>% as.numeric() %>% sum() / 10000-
-
-
st_difference(districts, pol2)[2] %>% plot(col = "grey", main = NULL)-
-
-
# Try compare the
-set.seed(1)
-farmers_alb %>% filter(season == 2) %>%
- sample_n(size = 5) %>% st_buffer(dist = 30000) %>% st_geometry %>% plot()
-
-# With
-set.seed(1)
-farmers_alb %>% filter(season == 2) %>%
- st_sample(size = 5) %>% st_buffer(dist = 30000) %>% st_geometry %>% plot()-
-
-
roads_gt400 <- roads %>% filter(as.numeric(st_length(.)) / 1000 > 400)
-par(mar = rep(0, 4))
-st_transform(districts, st_crs(roads)) %>% st_geometry %>%
- plot(col = "grey", main = NULL)
-roads_gt400 %>% st_buffer(25000) %>% plot(col = "green", add = TRUE)-
-
-
par(mar = rep(0, 4))
-st_transform(districts, st_crs(roads)) %>% st_geometry %>%
- plot(col = "grey", main = NULL)
-roads_gt400 %>% st_buffer(25000) %>% plot(col = "tan", add = TRUE)
-set.seed(1)
-roads_gt400 %>% st_buffer(25000) %>%
- st_sample(size = rep(20, nrow(.)), exact = TRUE) %>%
- plot(col = "red", pch = 20, add = TRUE)2 Module 2 practice answers
-2.1 Practice 1
-2.1.1 Questions
--
-
rasteruses the S4 object-oriented system, wheresfuses S4. Slots in S4 objects are accessed using the@operator.
-brick, becauserasteronly allows you to read and write a single layer.
-stackorganizes multiple rasters that might be stored in separation locations into a single multi-layer object, whereasbrickrequires a single file on disk.brickstake longer to create, but are faster to work with once they exist.
-The output of
raster’s vectorization functions arespobjects, so you have to convert them tosfobjects usingst_as_sf.
-
2.1.2 Code
--
-
-
# recreate r, r2, r3
-e <- extent(c("xmin" = 27, "xmax" = 29, "ymin" = -16, "ymax" = -14))
-r <- raster(x = e, res = 0.25, crs = crs(districts))
-set.seed(1)
-values(r) <- sample(1:100, size = ncell(r), replace = TRUE) # 3
-r2 <- r > 50
-r3 <- r
-set.seed(1)
-values(r3) <- rnorm(n = ncell(r3), mean = 10, sd = 2)
-
-#
-r4 <- r3
-set.seed(1)
-values(r4) <- runif(n = ncell(r4), 0, 1)
-r5 <- r4 > 0.5
-
-s2 <- list(r, r2, r3, r4, r5) %>% stack
-names(s2) <- c("r", "r2", "r3", "r4", "r5")
-plot(s2)-
-
- We use
tempdir()here, but you should use yournotebooks/datafolder.
-
b <- brick(s2, file = file.path(tempdir(), "b2.tif"))-
-
-
zamr3 <- raster(x = extent(districts), crs = crs(districts), res = 0.2)
-values(zamr3) <- 1:ncell(zamr3)
-
-farmersr <- farmers %>% distinct(uuid, .keep_all = TRUE) %>% select(x, y) %>%
- mutate(count = 1) %>% st_as_sf(coords = c("x", "y")) %>%
- rasterize(x = ., y = zamr3, field = "count", fun = sum)
-
-par(mar = c(0, 0, 1, 4))
-districts %>% st_union %>% plot(col = "grey", border = "grey")
-plot(farmersr, add = TRUE)-
-
-
zamr_alb <- projectRaster(from = zamr, res = 20000, crs = crs(roads),
- method = "ngb")
-farmersr_alb <- projectRaster(from = farmersr, to = zamr_alb,
- method = "bilinear")
-
-par(mar = c(0, 0, 1, 4), mfrow = c(1, 2))
-districts %>% st_union %>% plot(col = "grey", border = "grey")
-plot(farmersr, add = TRUE)
-districts %>% st_transform(crs(roads)) %>% st_union %>%
- plot(col = "grey", border = "grey")
-plot(farmersr_alb, add = TRUE)-
-
-
par(mar = c(0, 0, 0, 0))
-farmersr %>% rasterToPolygons(dissolve = TRUE) %>% st_as_sf %>%
- plot(main = NULL)2.2 Practice 2
-2.2.1 Questions
--
-
cellStatsprovides summary statistics over an entire layer;zonalcalculates statistics within pre-defined zones;focalcalculates statistics within a moving window.
-Use
method = bilinear.
-You don’t. You have to resample them to a common resolution and extent first. And then you have to
stackthem.
-
2.2.2 Code
--
-
-
as.Date("10-11-2017", "%m-%d-%Y")
-as.Date("10-11-17", "%m-%d-%y")
-as.Date("101117", "%m%d%y")
-as.Date("10112017", "%m%d%Y")
-lubridate::mdy("10-11-2017")
-lubridate::as_date("20171011")-
-
-
farmersr2 <- farmers %>% distinct(uuid, .keep_all = TRUE) %>%
- mutate(count = 1) %>% select(x, y, count) %>%
- st_as_sf(coords = c("x", "y")) %>%
- rasterize(., distsr, field = "count", fun = sum)-
-
-
zonal(farmersr2, distsr, fun = sum) %>% data.frame %>%
- subs(distsr, .) %>% plot_noaxes-
-
-
wmat3 <- matrix(1, nrow = 3, ncol = 3)
-wmat5 <- matrix(1, nrow = 5, ncol = 5)
-fstack <- list(sd3 = focal(x = chirpsz[[20]], w = wmat3, fun = sd),
- sd5 = focal(x = chirpsz[[20]], w = wmat5, fun = sd),
- max3 = focal(x = chirpsz[[20]], w = wmat3, fun = max),
- max5 = focal(x = chirpsz[[20]], w = wmat5, fun = max)) %>% stack
-fstack %>% plot_noaxes-
-
-
chirps_d57 <- chirpsz[[1]] %>% crop(., extent(districts[57, ]))
-s <- list(d1 = disaggregate(chirps_d57, fact = 5),
- d2 = disaggregate(chirps_d57, fact = 5, method = "bilinear")) %>%
- stack
-plot_noaxes(s)-
-
-
s <- lapply(list(mean, cv, median), function(x) {
- calc(chirpsz, fun = x)
-}) %>% stack
-names(s) <- c("Mean", "CV", "Median")
-plot_noaxes(s, nr = 1)2.3 Practice 3
-2.3.1 Questions
--
-
Using
cutwith vector of breakpoints (e.g. quantiles), or usereclassifywith a reclassification matrix.
-sampleRandom,sampleStratified. There is alsosampleRegular.
-
2.3.2 Code
--
-
-
chirps_sd <- calc(chirpsz, fun = sd)-
-
-
(chirps_sd < cellStats(chirps_sd, mean)) %>% plot-
-
-
quantile(raintot, probs = seq(0, 1, 0.2)) %>% cut(raintot, .) %>%
- plot_noaxes-
-
-
set.seed(11)
-randdistsr <- districts %>% sample_n(size = 15) %>% rasterize(., raintot)
-plot_noaxes(randdistsr)
-
-newrandrain <- mask(raintot, randdistsr)
-newrandrain %>% plot_noaxes-
-
-
set.seed(1)
-randsamp <- sampleRandom(x = newrandrain, size = 300)
-
-set.seed(1)
-ind <- sampleStratified(x = randdistsr, size = 300 / 15, cells = TRUE)
-stratsamp <- newrandrain[ind[, 1]]
-
-rand_rain_stats <- bind_rows(
- tibble(rain = randsamp, dat = "Simple"),
- tibble(rain = stratsamp, dat = "Stratified"),
-) %>% drop_na
-
-bp_theme <- theme(legend.title = element_blank(), axis.text.x = element_blank(),
- axis.ticks.x = element_blank(),
- panel.grid.major.x = element_blank(),
- panel.grid.minor.x = element_blank(),
- panel.background = element_rect(fill = "grey95"))
-
-rand_rain_stats %>% ggplot() +
- geom_boxplot(mapping = aes(y = rain, fill = dat), position = "dodge2") +
- scale_fill_manual(values = c("lightblue", "steelblue")) +
- ggtitle("Rainfall distributions") + xlab(NULL) + ylab("mm") + bp_theme2.4 Practice 4
-2.4.1 Questions
--
-
The
areafunction is your friend for this.
-Use
expressiontogether combined withpasteas needed for more complex labels.
-Make
names(predstack)matches the predictor names used by the model.
-
2.4.2 Code
--
-
-
demalb42 <- districts %>% filter(ID == 42) %>% st_transform(st_crs(roads)) %>%
- crop(demalb, .)
-vars <- c("slope", "aspect", "flowdir", "tri")
-terrvars <- stack(lapply(1:length(vars), function(x) {
- tv <- terrain(x = demalb42, opt = vars[x], unit = "degrees")
-}))
-names(terrvars) <- vars
-
-plot_noaxes(terrvars)-
-
-
library(gstat)
-
-# #1
-raintotalb <- projectRaster(from = raintot, res = 5000, crs = crs(roads))
-names(raintotalb) <- "rain"
-r <- raster(extent(raintotalb), res = res(raintotalb), crs = crs(raintotalb), vals = 1)
-
-# lapply approach to interpolation
-idw_list <- lapply(c(250, 500, 1000), function(x) {
- set.seed(1)
- rainsamp <- sampleRandom(raintotalb, size = 1000, xy = TRUE)
- rainsamp <- as.data.frame(rainsamp)
- invdist <- gstat(id = "rain", formula = rain ~ 1, locations = ~x + y,
- data = rainsamp)
- invdistr <- interpolate(object = r, model = invdist)
- invdistrmsk <- mask(x = invdistr, mask = raintotalb)
-})
-
-idws <- stack(c(raintotalb, idw_list))
-
-titles <- c("CHIRPS rainfall", "1000 pt IDW", "500 pt IDW", "250 pt IDW")
-plot_noaxes(idws, main = titles, zlim = c(0, 150))-
-
-
districts %>% filter(ID %in% seq(15, 50, 5)) %>%
- st_transform(st_crs(roads)) %>% st_geometry %>% st_centroid %>%
- as_Spatial(.) %>%
- distanceFromPoints(object = raintotalb, xy = .) %>%
- mask(., raintotalb) %>% plot_noaxes-
-
-
Pretty close to results of highest density sample. Slightly higher mean error.
-# #1
-data(zamprec)
-zamprecalb <- projectRaster(from = zamprec, to = raintotalb)
-names(zamprecalb) <- "rain"
-elev <- resample(aggregate(x = demalb, fact = 5), y = raintotalb)
-
-# #2
-set.seed(1)
-pts <- sampleRandom(x = zamprecalb, size = 25, sp = TRUE) %>% st_as_sf
-pts <- pts %>% mutate(elev = raster::extract(x = elev, y = .))
-pts_dat <- bind_cols(pts %>% data.frame %>% select(-geometry) %>% as_tibble,
- st_coordinates(pts) %>% as_tibble) %>% drop_na
-
-# #3
-p1 <- ggplot(pts_dat) + geom_point(aes(X, rain), col = "steelblue") +
- ylab("Rainfall (mm)")
-p2 <- ggplot(pts_dat) + geom_point(aes(Y, rain), col = "blue2") + ylab("")
-p3 <- ggplot(pts_dat) + geom_point(aes(elev, rain), col = "darkblue") + ylab("")
-cowplot::plot_grid(p1, p2, p3, nrow = 1)
-
-# #4
-rain_lm <- lm(rain ~ X + Y + elev, data = pts_dat)
-summary(rain_lm)
-
-# #5
-xs <- xFromCell(object = raintotalb, cell = 1:ncell(raintotalb))
-ys <- yFromCell(object = raintotalb, cell = 1:ncell(raintotalb))
-X <- Y <- raintotalb
-values(X) <- xs
-values(Y) <- ys
-
-# #6
-predst <- stack(X, Y, elev)
-names(predst) <- c("X", "Y", "elev")
-predrainr <- predict(object = predst, model = rain_lm)
-
-# #7
-s <- stack(zamprecalb, predrainr, (predrainr - zamprecalb) / zamprecalb * 100)
-mae <- round(cellStats(abs(zamprecalb - predrainr), mean), 1)
-
-pnames <- c("'Observed' Rainfall", "Predicted Rainfall", "% Difference")
-par(mfrow = c(1, 3), mar = c(0, 0, 1, 4))
-for(i in 1:3) {
- plot_noaxes(s[[i]], main = pnames[i])
- if(i %in% 1:2) {
- pts %>% st_geometry %>%
- plot(pch = 20, cex = 0.2, col = "grey70", add = TRUE)
- } else {
- mtext(side = 1, line = -3, cex = 0.8,
- text = paste("Mean abs err =", mae, "mm"))
- }
-}-
-