Lab4_PartA

Packages & Working Directory pkgs <- c(“sf”,“dplyr”,“ggplot2”,“tidyr”,“scales”,“RColorBrewer”,“cowplot”,“classInt”,“stringr”,“forcats”) to_install <- setdiff(pkgs, rownames(installed.packages())) if (length(to_install)) install.packages(to_install, dependencies = TRUE) invisible(lapply(pkgs, library, character.only = TRUE))
Set your Lab 4 folder (adjust if needed)

setwd(“C:/Users/ddtmd/Desktop/2025 Fall/GEOG588/Lab4”)

2.Load chapter 5 data

# read them (note the extra "Chapter5/Chapter5")
okcounty      <- sf::st_read("Chapter5/Chapter5/ok_counties.shp")

## Reading layer `ok_counties' from data source 
##   `C:\Users\ddtmd\Desktop\2025 Fall\GEOG588\Lab4\Chapter5\Chapter5\ok_counties.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 77 features and 7 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -103.0025 ymin: 33.62184 xmax: -94.43151 ymax: 37.00163
## Geodetic CRS:  NAD83

tpath  <- sf::st_read("Chapter5/Chapter5/ok_tornado_path.shp")

## Reading layer `ok_tornado_path' from data source 
##   `C:\Users\ddtmd\Desktop\2025 Fall\GEOG588\Lab4\Chapter5\Chapter5\ok_tornado_path.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 4092 features and 22 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -103 ymin: 33.15 xmax: -93.78 ymax: 38.1
## Geodetic CRS:  NAD83

tpoint <- sf::st_read("Chapter5/Chapter5/ok_tornado_point.shp")

## Reading layer `ok_tornado_point' from data source 
##   `C:\Users\ddtmd\Desktop\2025 Fall\GEOG588\Lab4\Chapter5\Chapter5\ok_tornado_point.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 4092 features and 22 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -103 ymin: 33.7535 xmax: -94.4381 ymax: 37
## Geodetic CRS:  NAD83

# Libraries
pkgs <- c("sf","dplyr","ggplot2","tidyr","scales","RColorBrewer","cowplot","classInt")
to_install <- setdiff(pkgs, rownames(installed.packages()))
if (length(to_install)) install.packages(to_install, dependencies = TRUE)
invisible(lapply(pkgs, library, character.only = TRUE))

## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE

## 
## 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

# Find the shapefiles (handles Chapter5/Chapter5 nesting)
shp_files <- list.files("Chapter5", pattern = "\\.shp$", recursive = TRUE, full.names = TRUE)

ok_counties_path <- shp_files[grepl("ok_counties\\.shp$",      shp_files, ignore.case = TRUE)]
ok_path_path     <- shp_files[grepl("ok_tornado_path\\.shp$",  shp_files, ignore.case = TRUE)]
ok_point_path    <- shp_files[grepl("ok_tornado_point\\.shp$", shp_files, ignore.case = TRUE)]

# Safety check
stopifnot(length(ok_counties_path) == 1, length(ok_path_path) == 1, length(ok_point_path) == 1)

# Read data
okcounty          <- sf::st_read(ok_counties_path, quiet = TRUE)
ok_tornado_path   <- sf::st_read(ok_path_path,     quiet = TRUE)
ok_tornado_point  <- sf::st_read(ok_point_path,    quiet = TRUE)

# Rename yr -> year BEFORE plotting
ok_tornado_path  <- ok_tornado_path  |> dplyr::rename(year = yr)
ok_tornado_point <- ok_tornado_point |> dplyr::rename(year = yr)

# (Optional) filter to 2016–2021
ok_tornado_path_f   <- dplyr::filter(ok_tornado_path,  year %in% 2016:2021)
ok_tornado_point_f  <- dplyr::filter(ok_tornado_point, year %in% 2016:2021)

# Quick base map to confirm
ggplot() + geom_sf(data = okcounty, fill = NA) + theme_bw()

# Composite figure (paths + points)
p_paths <- ggplot() +
  geom_sf(data = okcounty, fill = NA, color = "grey60", linewidth = 0.3) +
  geom_sf(data = ok_tornado_path_f, aes(color = factor(year)), linewidth = 0.7, alpha = 0.9) +
  scale_color_brewer(palette = "Set1", name = "Year") +
  coord_sf() + theme_minimal() +
  labs(title = "Oklahoma Tornado Paths (2016–2021)")

p_points <- ggplot() +
  geom_sf(data = okcounty, fill = NA, color = "grey60", linewidth = 0.3) +
  geom_sf(data = ok_tornado_point_f, aes(color = factor(year)), size = 1.3, alpha = 0.9) +
  scale_color_brewer(palette = "Set1", name = "Year") +
  coord_sf() + theme_minimal() +
  labs(title = "Oklahoma Tornado Points (2016–2021)")

cowplot::plot_grid(p_paths, p_points, ncol = 1, labels = "AUTO")

# County join + density (as you did)
tpoint_16_21 <- ok_tornado_point |> dplyr::filter(year %in% 2016:2021) |> dplyr::select(om, year, date)
countypnt <- sf::st_join(tpoint_16_21, okcounty) |> sf::st_drop_geometry()
countysum <- countypnt |> dplyr::group_by(GEOID) |> dplyr::summarize(tcnt = dplyr::n(), .groups = "drop")

countymap <- okcounty |>
  dplyr::left_join(countysum, by = "GEOID") |>
  dplyr::mutate(tcnt = tidyr::replace_na(tcnt, 0)) |>
  dplyr::mutate(area = sf::st_area(okcounty),
                tdens = 1e9 * tcnt / as.numeric(area))  # simple density scaling

ggplot(countymap) + geom_sf(aes(fill = tdens)) + theme_void() +
  labs(fill = "Tornado density")

# --- Prep: ensure we have a 'year' column in both layers ---
if ("yr" %in% names(ok_tornado_path))  ok_tornado_path  <- dplyr::rename(ok_tornado_path,  year = yr)
if ("yr" %in% names(ok_tornado_point)) ok_tornado_point <- dplyr::rename(ok_tornado_point, year = yr)

# OPTIONAL: restrict to 2016–2021 to match the figure request
ok_tornado_path_f  <- dplyr::filter(ok_tornado_path,  year %in% 2016:2021)
ok_tornado_point_f <- dplyr::filter(ok_tornado_point, year %in% 2016:2021)

# use the same palette and legend title in both plots
pal_name <- "Set1"

# Map 1: tornado paths by year (different color per year)
p_paths <- ggplot() +
  geom_sf(data = okcounty, fill = NA, color = "grey60", linewidth = 0.3) +
  geom_sf(data = ok_tornado_path_f, aes(color = factor(year)), linewidth = 0.8, alpha = 0.9) +
  scale_color_brewer(palette = pal_name, name = "Year") +
  coord_sf(datum = NA) +
  theme_minimal(base_size = 12) +
  theme(
    panel.grid.major = element_blank(),
    legend.position  = "none"   # hide here; we'll show legend on the points map
  ) +
  labs(title = "Oklahoma Tornado Paths (2016–2021)")

# Map 2: tornado points by year (same colors)
p_points <- ggplot() +
  geom_sf(data = okcounty, fill = NA, color = "grey60", linewidth = 0.3) +
  geom_sf(data = ok_tornado_point_f, aes(color = factor(year)), size = 1.3, alpha = 0.9) +
  scale_color_brewer(palette = pal_name, name = "Year") +
  coord_sf(datum = NA) +
  theme_minimal(base_size = 12) +
  theme(
    panel.grid.major = element_blank(),
    legend.position  = "right"
  ) +
  labs(title = "Oklahoma Tornado Points (2016–2021)")

# Composite figure (paths on top, points below)
cowplot::plot_grid(
  p_paths, p_points,
  ncol = 1, labels = "AUTO", align = "hv", rel_heights = c(1, 1.05)
)

#Summarize the density of tornado points by both county and year and generate a faceted plot that displays maps of county-level tornado density from 2016-2021.

# 1) Ready
if ("yr" %in% names(tpoint)) tpoint <- dplyr::rename(tpoint, year = yr)
if ("yr" %in% names(tpath))  tpath  <- dplyr::rename(tpath,  year = yr)
tpoint <- dplyr::filter(tpoint, year %in% 2016:2021)
tpath  <- dplyr::filter(tpath,  year %in% 2016:2021)

# 2) Two maps
p_paths  <- ggplot() + geom_sf(data = okcounty, fill = NA) +
  geom_sf(data = tpath,  aes(color = factor(year))) + labs(title="Paths")
p_points <- ggplot() + geom_sf(data = okcounty, fill = NA) +
  geom_sf(data = tpoint, aes(color = factor(year))) + labs(title="Points")

# 3) Combine
cowplot::plot_grid(p_paths, p_points, ncol = 1, labels = "AUTO")

2. Density facet map by year and county

library(sf)
library(dplyr)
library(ggplot2)
library(classInt)
library(cowplot)
library(RColorBrewer)

# 1) point -> county 공간조인 후 연/카운티 집계
pny <- sf::st_join(dplyr::select(tpoint, year), okcounty) |>
  sf::st_drop_geometry() |>
  dplyr::count(GEOID, year, name = "tcnt")

# 2) 카운티 폴리곤에 붙이고 면적/밀도 계산 (연도별)
countymap_year <- okcounty |>
  dplyr::left_join(pny, by = "GEOID") |>
  dplyr::mutate(tcnt = tidyr::replace_na(tcnt, 0)) |>
  dplyr::group_by(year) |>
  dplyr::mutate(
    area_m2 = as.numeric(sf::st_area(geometry)),     # 면적(m^2)
    tdens   = tcnt / area_m2 * 1e9                   # 10^9 m^2 당 개수(해석 쉬운 스케일)
  ) |>
  dplyr::ungroup()

# 3) 2016–2021만 보이도록(원하면)
countymap_year <- dplyr::filter(countymap_year, year %in% 2016:2021)

# 4) 연도별 패싯 지도
p_facet <- ggplot(countymap_year) +
  geom_sf(aes(fill = tdens), color = NA) +
  facet_wrap(~ year, ncol = 3) +
  scale_fill_distiller(palette = "YlOrRd", direction = 1, na.value = "grey90") +
  theme_void() +
  labs(title = "County-level Tornado Point Density (2016–2021)",
       fill = "density\n(points / 1e9 m²)")
p_facet

Quantile classification 4 sheets with 3 to 6 classes

# 한 장의 지도(분위수 nclass 지정) 만들기 함수
make_quantile_map <- function(nclass) {
  # 연도 합산된 밀도(또는 연도 무시하고 총합 밀도 비교용) — 문제 지시에 맞게 연도별이 아닌 '총 밀도' 비교용으로 구성
  sum_by_county <- pny |>
    dplyr::group_by(GEOID) |>
    dplyr::summarise(tcnt = sum(tcnt), .groups = "drop")

  cm_sum <- okcounty |>
    dplyr::left_join(sum_by_county, by = "GEOID") |>
    dplyr::mutate(tcnt   = tidyr::replace_na(tcnt, 0),
                  area_m2 = as.numeric(sf::st_area(geometry)),
                  tdens   = tcnt / area_m2 * 1e9)

  # 분위수 구간 만들기
  brks <- classInt::classIntervals(cm_sum$tdens, n = nclass, style = "quantile")$brks
  cm_sum <- cm_sum |>
    dplyr::mutate(class_q = cut(tdens, breaks = brks, include.lowest = TRUE, dig.lab = 6))

  ggplot(cm_sum) +
    geom_sf(aes(fill = class_q), color = NA) +
    theme_void() +
    scale_fill_brewer(palette = "YlOrRd", na.value = "grey90", drop = FALSE) +
    labs(title = paste0("Quantile breaks (k = ", nclass, ")"),
         fill = "tdens\nquantile")
}

p3 <- make_quantile_map(3)
p4 <- make_quantile_map(4)
p5 <- make_quantile_map(5)
p6 <- make_quantile_map(6)

cowplot::plot_grid(p3, p4, p5, p6, labels = "AUTO", ncol = 2)

Lab4_PartA

Cindy Lee

2025-11-01