Module12_Exercises: Spatial Data Analysis

Setup

library(sf)
library(terra)
library(ggplot2)
library(dplyr)
library(viridis)

Exercise 1: Countries Map

# Load the countries shapefile
countries <- st_read("./data/countries/countries.shp", quiet = TRUE)

# Check the variables
names(countries)

 [1] "SP_ID"      "scalerank"  "featurecla" "labelrank"  "sovereignt"
 [6] "sov_a3"     "adm0_dif"   "level"      "type"       "admin"     
[11] "adm0_a3"    "geou_dif"   "geounit"    "gu_a3"      "su_dif"    
[16] "subunit"    "su_a3"      "brk_diff"   "name"       "name_long" 
[21] "brk_a3"     "brk_name"   "brk_group"  "abbrev"     "postal"    
[26] "formal_en"  "formal_fr"  "note_adm0"  "note_brk"   "name_sort" 
[31] "name_alt"   "mapcolor7"  "mapcolor8"  "mapcolor9"  "mapcolor13"
[36] "pop_est"    "gdp_md_est" "pop_year"   "lastcensus" "gdp_year"  
[41] "economy"    "income_grp" "wikipedia"  "fips_10"    "iso_a2"    
[46] "iso_a3"     "iso_n3"     "un_a3"      "wb_a2"      "wb_a3"     
[51] "woe_id"     "adm0_a3_is" "adm0_a3_us" "adm0_a3_un" "adm0_a3_wb"
[56] "continent"  "region_un"  "subregion"  "region_wb"  "name_len"  
[61] "long_len"   "abbrev_len" "tiny"       "homepart"   "geometry"  

# Create cleaned variables for mapping
countries_map <- countries %>%
  mutate(
    gdp_clean = ifelse(gdp_md_est <= 0, NA, gdp_md_est),
    pop_clean = ifelse(pop_est <= 0, NA, pop_est),
    log_gdp = log10(gdp_clean),
    log_pop = log10(pop_clean)
  )

# Map 1: Median GDP
ggplot(countries_map) +
  geom_sf(aes(fill = log_gdp), color = "gray70", linewidth = 0.1) +
  scale_fill_viridis_c(option = "magma", na.value = "gray90") +
  labs(
    title = "Global Median GDP by Country",
    fill = "Log10 GDP"
  ) +
  theme_void()

# Map 2: Population
ggplot(countries_map) +
  geom_sf(aes(fill = log_pop), color = "gray70", linewidth = 0.1) +
  scale_fill_viridis_c(option = "viridis", na.value = "gray90") +
  labs(
    title = "Global Population by Country",
    fill = "Log10 Population"
  ) +
  theme_void()

# Optional categorical map for income group
ggplot(countries) +
  geom_sf(aes(fill = factor(income_grp)), color = "gray70", linewidth = 0.1) +
  scale_fill_brewer(palette = "Set3", na.value = "gray90") +
  labs(
    title = "Income Group by Country",
    fill = "Income Group"
  ) +
  theme_void()

Exercise 2: Median NDVI for Bethel Island and Oakley

# Load the raster bands used to calculate NDVI
b4 <- rast("./data/rs/LC08_044034_20170614_B4.tif")
b5 <- rast("./data/rs/LC08_044034_20170614_B5.tif")

# Calculate NDVI
ndvi <- (b5 - b4) / (b5 + b4)
names(ndvi) <- "ndvi"

# Plot NDVI if needed
plot(ndvi, col = rev(terrain.colors(10)), main = "NDVI")

# Load the California places shapefile
ca_places <- st_read("./data/ca_places/ca_places.shp", quiet = TRUE)

# Match the vector CRS to the NDVI raster CRS
ca_places <- st_transform(ca_places, crs(ndvi))

# Select Bethel Island and Oakley
places_two <- ca_places %>%
  filter(NAME %in% c("Bethel Island", "Oakley"))

# Check the selected place names
places_two$NAME

[1] "Oakley"        "Bethel Island"

# Optionally plot the two places on top of NDVI
plot(ndvi, col = rev(terrain.colors(10)), main = "Bethel Island and Oakley on NDVI")
plot(st_geometry(places_two), add = TRUE, border = "red", lwd = 2)

# Extract median NDVI for each place
median_ndvi <- terra::extract(
  ndvi,
  terra::vect(places_two),
  fun = median,
  na.rm = TRUE
)

# Attach place names and clean the output
median_ndvi <- median_ndvi %>%
  mutate(place = places_two$NAME[ID]) %>%
  select(place, ndvi)

# 11) Show final result
median_ndvi

          place      ndvi
1        Oakley 0.3005387
2 Bethel Island 0.4180438