Module12 (6/2/2026)

Author

Ethan Schatz

This module explores spatial framing packages for R.

Part 1

Load the Data

library(dplyr)


Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(sf)

Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

library(terra)

terra 1.9.27

library(ggplot2)

countries <- st_read("E:/Summer 2026/data/countries/countries.shp", quiet = TRUE)
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"

Plot Median GDP

countries$gdp_plot <- ifelse(
  countries$gdp_md_est > 0,
  countries$gdp_md_est,
  NA
)

gdp_breaks <- c(1e3, 1e4, 1e5, 1e6, 1e7)

ggplot() +
  geom_sf(data = countries, aes(fill = gdp_plot)) +
  scale_fill_gradient(
    low = "white",
    high = "darkblue",
    trans = "log10",
    breaks = gdp_breaks,
    labels = c("$1B", "$10B", "$100B", "$1T", "$10T"),
    na.value = "gray80"
  ) +
  theme_bw()

Plot Population

countries$pop_plot <- ifelse(
  countries$pop_est > 0,
  countries$pop_est,
  NA
)

pop_breaks <- c(1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9)

ggplot() +
  geom_sf(data = countries, aes(fill = pop_plot)) +
  scale_fill_gradient(
    low = "white",
    high = "darkgreen",
    trans = "log2",
    breaks = pop_breaks,
    labels = c("1,000", "10,000", "100,000", "1,000,000", "10,000,000", "100,000,000", "1,000,000,000"),
    na.value = "gray80"
  ) +
  theme_bw()

Part 2

Load the data

ca_places <- st_read("E:/Summer 2026/data/ca_places/ca_places.shp", quiet = TRUE)

bethel <- ca_places %>% 
  dplyr::filter(NAME %in% c("Bethel Island", "Oakley"))

# Blue
b2 <- rast("E:/Summer 2026/data/rs/LC08_044034_20170614_B2.tif")

# Green
b3 <- rast("E:/Summer 2026/data/rs/LC08_044034_20170614_B3.tif")

# Red
b4 <- rast("E:/Summer 2026/data/rs/LC08_044034_20170614_B4.tif")

# Near Infrared
b5 <- rast("E:/Summer 2026/data/rs/LC08_044034_20170614_B5.tif")

Calculate NDVI

ndvi <- (b5 - b4) / (b5 + b4)

bethel_vect <- vect(bethel)
bethel_vect <- project(bethel_vect, crs(ndvi))

median_ndvi <- extract(
  ndvi,
  bethel_vect,
  fun = median,
  na.rm = TRUE
)

median_ndvi$NAME <- bethel$NAME

median_ndvi <- median_ndvi |>
  dplyr::rename(Median = 2) |>
  dplyr::select(NAME, Median)

median_ndvi

           NAME    Median
1        Oakley 0.3005387
2 Bethel Island 0.4180438

DISCLAIMER:

ChatGPT was used during the process of writing the code for the purpose of debugging and fixing errors in the code.