Spatial Data Analysis (450:320)

library(geodata)
library(dplyr)
library(sf)

tza <- rnaturalearth::ne_countries(country = "United Republic of Tanzania")

car_plants <- read.csv("https://sselp.github.io/sdawr//Riched_Carnivorous_Plant.csv") %>% 
    st_as_sf(coords = c("longitude", "latitude"), crs = 4326)

pop_global <- population(
  year = 2020, res = 5,
  path = tempdir())

pop_global <- log(pop_global + 1)

# Overlay: crop and mask, compare the difference
pop_tza1 <- mask(pop_global, tza)
pop_tza2 <- mask(crop(pop_global, tza), tza)

# Extract, compare the difference
pop_plants <- extract(pop_global, car_plants)
pop_plants <- extract(pop_global, car_plants, cells = TRUE, xy = TRUE)
pop_plants <- extract(pop_global, car_plants, bind = TRUE)

fname <- 'your_class_note_folder/S2A_MSIL2A_20201014T154231_N0214_R011_T19TBG_20201014T200156_20m.tif'
s2 <- rast(fname)

names(s2) <- paste0('B', c(2:7, 11:12, 8))

ndvi <- (s2$B8 - s2$B4) / (s2$B8 + s2$B4)
water_mask <- ndvi < 0
plot(water_mask)

# A clear way to define the reclassify matrix
rcl <- matrix(
  c(-1, 0, 1, # class 1
    0, 0.2, 2, # class 2
    0.2, 0.5, 3, # class 3
    0.5, 1, 4), # class 4
  ncol = 3, byrow = T)

lc <- classify(ndvi, rcl, include.lowest = TRUE)

plot(lc)
lc

# Convert to categorical map
lc <- as.factor(lc)
levels(lc) <- data.frame(
  value = 1:4,
  landcover = c(
    "water",
    "Urban",
    "sparse vegetation",
    "dense vegetation"
  )
)

plot(lc)
lc
levels(lc)

# Aggregate, result is a raster with coarse resolution
pop_tza_coarse <- aggregate(pop_tza, fact = 8, fun = mean)

# Zonal, result is a matrix! 
cnts <- rnaturalearth::ne_countries()
cnts_rst <- rasterize(cnts, pop_global, "name")
pop_zonal <- zonal(pop_global, cnts_rst, fun = mean)

# Observe values in pop_zonal, what might be wrong? 

# Focal, result is a raster as well.
## Use s2 NIR band as an example
fname <- 'your_class_note_folder/S2A_MSIL2A_20201014T154231_N0214_R011_T19TBG_20201014T200156_20m.tif'
nir_band <- rast(fname, lyrs = 9)

# weight matrix could not the square, 
# but ncol and nrow must be uneven numbers

## Reproduce Sobel edge detection
wx <- matrix(c(-1, -2, -1, 0, 0, 0, 1, 2, 1), nrow = 3, byrow = TRUE)
wy <- matrix(c(1, 0, -1, 2, 0, -2, 1, 0, -1), nrow = 3, byrow = TRUE)

# Apply focal calculation
nir_x <- focal(nir_band, w = wx)
nir_y <- focal(nir_band, w = wy)

# Simply band calculation
nir_edge <- sqrt(nir_x^2 + nir_y^2)

writeRaster(nir_edge, 'your_class_note_folder/edge.tif')

# Calculate the distance to the carnivore plant observations
car_plants_rst <- rasterize(car_plants, field = 1, pop_tza2)
dis_plants <- terra::distance(car_plants_rst)