Geospatial Analysis with R

class: center, middle, inverse, title-slide

.title[
# Geospatial Analysis with R
]
.subtitle[
## Class 22
]

---

## Today

- Modeling review (random forest)
- Project update

---
## Modelling demo
### Random Forests

```r
# Train model
library(caTools)
library(randomForest)
library(geospaar)

# read in train set
train_ref <- system.file("extdata/train_reference.csv", 
                         package = "geospaar") %>% 
  read_csv()  %>% mutate(crop = factor(case_when(
  class == "Maize" ~ 1,
  class == "Rice" ~ 2,
  class == "Other" ~ 3,
  class == "noncrop" ~ 4
))) %>% dplyr::select(id, class, crop, !!names(.))
```

---
### Model

```r
# prepare dataset
dset <- train_ref %>% filter(type == "mean") %>%
  select(-id, -class, -type, -seed)
set.seed(10)
samp <- sample.split(dset$crop, SplitRatio = 0.8)
train <- subset(dset, samp == TRUE)
test <- subset(dset, samp == FALSE)
# train %>% tidyr::drop_na()

# fit model
set.seed(123)
mod <- randomForest(crop ~ ., data = train, ntree = 1000,
                    importance = TRUE)
sum(diag(mod$confusion[1:4, 1:4])) / sum(mod$confusion[1:4, 1:4])
summary(mod)

# Evaluation
pred <- predict(mod, newdata = test %>% dplyr::select(-c(crop)))
cm <- table(test$crop, pred)
# sum(diag(cm)) / sum(cm)

rf_results <- list("train" = train, "test" = test, "mod" = mod, "pred" = pred,
                   "error_mat" = cm)
```

---

```r
imptab <- mod$importance
nms <- rownames(imptab)
imp_tbl <- as_tibble(imptab) %>% 
  mutate(Variable = nms) %>% 
  dplyr::select(Variable, MeanDecreaseAccuracy) %>% 
  arrange(MeanDecreaseAccuracy) %>% 
  mutate(order = 1:n()) %>% 
  mutate(Variable = gsub("\\.", "_", Variable))

p <- ggplot(imp_tbl) + 
  geom_point(aes(MeanDecreaseAccuracy, factor(order))) + 
  scale_y_discrete(labels = toupper(imp_tbl$Variable)) + 
  ylab("Variable") + xlab("Mean Decrease in Accuracy") + 
  theme_linedraw() + 
  theme(axis.text = element_text(size = 7), 
        axis.title = element_text(size = 7))
p
```

---

## Map

Look at Rmd file to see code use to make prediction map for crop types.

<img src="figures/crop_probs.png" width="60%" style="display: block; margin: auto;" />
---
## Advanced - parallelization

```r
library(parallel)
# detectCores()
# cross-platform
cl <- makeCluster(4)
clusterEvalQ(cl, library(raster))
system.time(b <- parLapply(cl, 1:2000, function(x) {  # x <- 1
  r <- raster::raster(extent(30, 31, -14, -13), res = 0.01, 
                      crs = "+proj=longlat +datum=WGS84") 
  set.seed(x)
  values(r) <- sample(1:10, size = ncell(r), replace = TRUE)
  stack(r, r * runif(n = ncell(r), 0.9, 1.2), r * runif(n = ncell(r), 0.8, 1.3))
}))
stopCluster(cl)

# *nix only
library(terra)
system.time(bmc <- mclapply(1:2000, function(x) {
  
  r <- raster(extent(30, 31, -14, -13), res = 0.01, 
              crs = "+proj=longlat +datum=WGS84")
  set.seed(x)
  values(r) <- sample(1:10, size = ncell(r), replace = TRUE)
  stack(r, r * runif(n = ncell(r), 0.9, 1.2), r * runif(n = ncell(r), 0.8, 1.3))
}, mc.cores = 4))

```
---

```r
# serial
system.time(b2 <- lapply(1:2000, function(x) {
  r <- raster(extent(30, 31, -14, -13), res = 0.01, 
              crs = "+proj=longlat +datum=WGS84")
  set.seed(x)
  values(r) <- sample(1:10, size = ncell(r), replace = TRUE)
  stack(r, r * runif(n = ncell(r), 0.9, 1.2), r * runif(n = ncell(r), 0.8, 1.3))
}))
```

## Advanced plotting

---
### mapView

```r
library(mapview)
mapview(districts)
chirps <- terra::rast(system.file("extdata/chirps.tif", package = "geospaar"))

viewRGB(x = terra::as.raster(chirps[[1:3]]))
# mapView(districts, alpha.regions = 1, legend = FALSE)
mapView(districts, alpha.regions = 0, legend = FALSE) + mapview(raintot)

kili_data <- system.file("extdata", "kiliNDVI.tif", package = "mapview") %>%
  raster::stack(.)
mapview(kili_data[[c(1, 10, 23)]])
```

---
### tmap

```r
library(tmap)
data(World)
# tmaptools::palette_explorer()
tm_shape(raintot) + tm_raster(palette = "magma", breaks = seq(0, 200, 25)) +
  tm_shape(districts) + tm_borders() + 
  tm_layout(bg.color = "grey", inner.margins = 0.1) # +
  # tm_shape(World) + tm_borders()
```
---

### leaflet

```r

library(leaflet)
m <- leaflet() %>% addTiles()
# # m  # a map with the default OSM tile layer
# 
m %>% addProviderTiles("Esri.WorldImagery")
# 
# # set bounds
# m %>% fitBounds(0, 40, 10, 50)
# #
# # # move the center to Snedecor Hall
# m <- m %>% setView(-93.65, 42.0285, zoom = 17)
# m
# #
# # popup
# # m %>% addPopups(-93.65, 42.0285,
# #                 "Here is the <b>Department of Statistics</b>, ISU")
rand_lng <- function(n = 10) rnorm(n, -93.65, .01)
rand_lat <- function(n = 10) rnorm(n, 42.0285, .01)

# use automatic bounds derived from lng/lat data
m <- m %>% clearBounds()

# popup
m %>% addPopups(rand_lng(), rand_lat(), "Random popups")
```
---