1 Map 1: Decompose and Recreate a Professional Map

1.1 Original Map

The map I chose to recreate is a severe weather risk forecast published by The New York Times, showing nested storm risk zones across the Midwest using semi-transparent categorical polygons layered over a minimal basemap.

NYT Severe Weather Risk Map Source: The New York Times. “Strong Storms Bring a Threat of Tornadoes on Thursday”

1.2 Decomposition

Map type: Categorical risk zone map using overlapping filled polygons
Projection: Likely Albers Equal Area or Web Mercator (Mapbox default)
Classification: Ordinal — three categories (Some, Moderate, High risk); no numeric classification needed
Color palette: Sequential warm — yellow → orange → dark red; communicates escalating danger intuitively
Legend: Horizontal inline legend in the title area, not a traditional box; very clean
Typography: Sans-serif throughout; city labels are small and subordinate to the risk zones
Annotation: Italic label placed inside the risk blob naming the phenomenon
Contextual layers: State/county boundaries in light gray, city dots + labels, water bodies; no satellite or terrain
White space: Basemap is nearly empty — all visual weight goes to the risk polygons

1.3 Recreation

packages <- c("ggplot2", "mapsf", "tigris", "dplyr", "tidyverse", 
              "rnaturalearth", "sf", "glue", 
              "ggspatial", "cowplot", "patchwork", "scales",
              "readr", "prettydoc")
lapply(packages, library, character.only = TRUE)

paths <- st_read("data/1950-2024-torn-aspath/1950-2024-torn-aspath.shp")

paths_filt <- paths %>%
  filter(mo == 5, yr == 2011, dy == 22, mag == 5)

states <- states(cb = TRUE, resolution = "20m") %>%
  st_transform(5070) %>%
  filter(STUSPS %in% c("MO","KS","OK","AR","IL"))

counties <- counties(cb = TRUE, resolution = "20m") %>%
  st_transform(5070) %>%
  filter(STUSPS %in% c("MO","KS","IL","AR","OK"))

cities <- data.frame(
  name = c("Joplin","Kansas City","St. Louis","Springfield","Tulsa"),
  lon  = c(-94.51, -94.58, -90.20, -93.29, -95.99),
  lat  = c( 37.08,  39.10,  38.63,  37.21,  36.15)
) %>%
  st_as_sf(coords = c("lon","lat"), crs = 4326) %>%
  st_transform(5070)

paths_filt <- paths_filt %>%
  st_transform(5070) %>%
  mutate(mag_f = factor(mag, levels = 0:5,
                        labels = c("EF0","EF1","EF2","EF3","EF4","EF5")))

joplin_bbox <- st_bbox(c(
  xmin = -80000, xmax = 350000,
  ymin = 1380000, ymax = 1750000
), crs = st_crs(5070))

joplin_buff1 <- paths_filt %>% st_buffer(25000) %>% st_union()
joplin_buff2 <- paths_filt %>% st_buffer(60000) %>% st_union()

rivers <- ne_download(scale = 10, type = "rivers_lake_centerlines",
                      category = "physical", returnclass = "sf") %>%
  st_transform(5070) %>%
  st_crop(st_bbox(states))

joplin_length_mi <- paths_filt %>%
  st_length() %>% as.numeric() / 1000 * 0.621371

label_pt <- paths_filt %>%
  st_centroid() %>% st_coordinates() %>% as.data.frame() %>%
  mutate(X = X + 30000, Y = Y + 40000) %>%
  st_as_sf(coords = c("X","Y"), crs = 5070)
label_pt$label <- "EF5"

mf_init(joplin_bbox, expandBB = c(0.3, 0.3, 0.3, 0.3))

map_title <- glue("The Joplin EF5 Tornado carved {round(joplin_length_mi, 1)} miles through\nsouthwest Missouri on May 22, 2011")

mf_map(states, type = "base", col = "#EBEBDF", border = "#BBBBAA", 
       lwd = 0.5, add = TRUE)

mf_map(counties, type = "base", col = NA, border = "#D5D5C8", 
       lwd = 0.2, add = TRUE)

# rivers as blue lines
mf_map(rivers, col = "#A8C8E0", lwd = 0.8, add = TRUE)

mf_map(joplin_buff2, type = "base", 
       col = "#f6e27a",
       border = NA,
       alpha = 0.4,
       add = TRUE
       )


mf_map(joplin_buff1, type = "base", 
       col = "#e08c3a",
       border = NA,
       alpha = 0.4,
       add = TRUE
       )



mf_map(paths_filt, type = "base", 
       col = "#c0392b",
       lwd = 8,  # Joplin EF5 should be thick
       add = TRUE)


# city dots
mf_map(cities, type = "base", pch = 20, cex = 0.8, 
       col = "#444444", add = TRUE)

mf_label(cities, var = "name", 
         halo = TRUE, r = 0.15,
         cex = 0.5, col = "#222222", overlap = FALSE)

# state abbreviation labels
mf_label(states, var = "STUSPS",
         col = "#999988", cex = 0.5, overlap = FALSE)

mf_layout(
  title = map_title,
  credits = "Data: NOAA Storm Prediction Center, 1950–2024 Tornado Paths",
  scale = TRUE,
  frame = FALSE
)
  # EF5 annotation
mf_label(label_pt, var = "label",
         col = "#7b0000",
         cex = 0.9,
         font = 3,  # italic
         halo = TRUE, r = 0.1)

mf_label(paths_filt, var = "label",
         pch  = 20,
       cex  = 1.2,
       add  = TRUE,
       leg_title    = "",           # removes legend title
       leg_pos      = "bottomright", # moves to bottom right
       leg_box_col  = "#FFFFFF",    # white bounding box fill
       leg_box_border = "#AAAAAA")   # grey border around box))

## Warning in text.default(x = cc[, 1], y = cc[, 2], labels = words, cex = cex, :
## "add" is not a graphical parameter

## Warning in text.default(x = cc[, 1], y = cc[, 2], labels = words, cex = cex, :
## "leg_title" is not a graphical parameter

## Warning in text.default(x = cc[, 1], y = cc[, 2], labels = words, cex = cex, :
## "leg_pos" is not a graphical parameter

## Warning in text.default(x = cc[, 1], y = cc[, 2], labels = words, cex = cex, :
## "leg_box_col" is not a graphical parameter

## Warning in text.default(x = cc[, 1], y = cc[, 2], labels = words, cex = cex, :
## "leg_box_border" is not a graphical parameter

1.4 Reflection

I researched the correct project for a Southern U.S., and tested between ESPG:5070 and ESPG 26915. I ended up selecting the former. I had issues with considering how many states to include. Since I eventually settled on setting a bounding box for the Joplin surrounding region anyway, I probably could re-add information from Tennessee. The only data I downloaded was line path data from NOAA. I also have to keep in mind that they are using MapBox and I am using map_sf for my annotations. So other than that, I did a great job referencing different geographical components.Could have incorporated more estimation information. My goal was to document the path of the Joplin tornado, and maybe in a project I would have a section describing furhter statistics about the sheer power of this storm.

2 Map 2: Projection as a Design Choice

My practicum project revolves around South African population projection, so I am using some preliminary data to showcase distortion on South Africa. The census is based on 2011 ward-level (akin to say, a county).

2.1 Maps

# Load data
wards <- read_csv("data/Untitled_2026-04-04-0953.csv") %>%
  st_as_sf(coords = c("LONGITUDE", "LATITUDE"), crs = 4326)

# aggregate est ward pop to province
province_pop <- wards %>%
  st_drop_geometry() %>%
  group_by(PROVINCE) %>%
  summarise(total_pop = sum(EST_WARD_POP, na.rm = TRUE))

sa_provinces <- ne_download(
  scale = 50, 
  type = "states", 
  category = "cultural", 
  returnclass = "sf"
) %>%
  filter(admin == "South Africa")

neighbors_sa <- ne_countries(continent = "africa", returnclass = "sf") %>%
  filter(admin %in% c("Botswana", "Namibia", "Lesotho", "Eswatini", "Mozambique"))

province_pop <- wards %>%
  st_drop_geometry() %>%
  group_by(PROVINCE) %>%
  summarise(total_pop = sum(EST_WARD_POP, na.rm = TRUE))

sa_provinces <- sa_provinces %>%
  left_join(province_pop, by = c("name" = "PROVINCE"))

table <- sa_provinces %>% st_drop_geometry() %>% select(name, total_pop)
table(table)

province_centroids <- sa_provinces %>%
  st_centroid() %>%
  st_transform("ESRI:102022")

# you need these objects ready
sa_provinces_albers <- sa_provinces %>% 
  st_transform("ESRI:102022") %>%
  st_make_valid() %>%
  mutate(
    area_km2   = as.numeric(st_area(.)) / 1e6,
    pop_density = total_pop / area_km2
  )
sa_provinces_wgs <- sa_provinces %>%
  st_make_valid() %>%
  mutate(
    area_km2    = as.numeric(st_area(sa_provinces_albers)) / 1e6,
    pop_density = total_pop / area_km2
  )

neighbors_sa_albers <- neighbors_sa %>% st_transform("ESRI:102022")
neighbors_sa_wgs <- neighbors_sa

# water — Indian + Atlantic Ocean context
ocean <- ne_download(scale = 50, type = "ocean", 
                     category = "physical", returnclass = "sf")

# neighbor centroids for labels
neighbor_labels_albers <- neighbors_sa %>%
  st_centroid() %>%
  st_transform("ESRI:102022")

sa_bbox_albers <- st_bbox(sa_provinces_albers)

expanded_bbox <- st_bbox(c(
  xmin = as.numeric(sa_bbox_albers["xmin"]) - 200000,  # less west padding
  xmax = as.numeric(sa_bbox_albers["xmax"]) + 200000,  # much less east — drops Madagascar
  ymin = as.numeric(sa_bbox_albers["ymin"]) - 8000,  # less south
  ymax = as.numeric(sa_bbox_albers["ymax"]) + 300000   # enough north for Botswana label
), crs = st_crs("ESRI:102022"))


map_theme <- theme_void() +
  theme(
    plot.title    = element_text(face = "bold", size = 11, color = "#222222"),
    plot.subtitle = element_text(size = 8, color = "#666655"),
    legend.position = "right"
  )

fill_scale <- scale_fill_viridis_c(
  name = "People per km²",
  labels = scales::comma,
  option = "viridis"
)
# WGS84 version
# WGS84
map_wgs84 <- ggplot() +
  geom_sf(data = neighbors_sa, fill = "#E8E8E0", color = "#BBBBAA", linewidth = 0.3) +
  geom_sf(data = sa_provinces_wgs, aes(fill = pop_density), color = "#BBBBAA", linewidth = 0.3) +
geom_sf_text(data = sa_provinces_wgs %>% st_centroid(),
             aes(label = name),
             size = 1.8, color = "white", fontface = "italic") +
  geom_sf_text(data = neighbors_sa %>% 
                 st_transform("EPSG:4326") %>% st_centroid(),
               aes(label = admin),
               size = 1.8, color = "#666655") +
  fill_scale +
labs(title = "WGS84 Projection",
     subtitle = "Distorts province shapes")+
  theme_void()

## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries

map_wgs84 + map_theme

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

fill_scale <- scale_fill_viridis_c(
  name = "People per km²",
  labels = scales::comma,
  option = "viridis"
)

map_albers <- ggplot() +
  geom_sf(data = neighbors_sa %>% st_transform("ESRI:102022"),
          fill = "#E8E8E0", color = "#BBBBAA", linewidth = 0.3) +
  geom_sf(data = sa_provinces_albers, 
          aes(fill = pop_density), color = "#BBBBAA", linewidth = 0.3) +
  geom_sf_text(data = sa_provinces_albers %>% st_centroid(),
               aes(label = name),
               size = 1.8, color = "white", fontface = "italic") +
  geom_sf_text(data = neighbors_sa %>% 
                 st_transform("ESRI:102022") %>% st_centroid(),
               aes(label = admin),
               size = 1.8, color = "#666655") +
  fill_scale +
  labs(title = "Africa Albers Equal Area") +
  theme_void()

## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries

map_albers + map_theme

2.2 Reflection

I selected WGS84 and Africa Albers Equal Area to look at population data for South Africa. WGS84 elongates South Africa and its neighbors, but Africa Albers Eqaul Area lets the geographies sit a bit more comfortably on the map. I feel like WGS84 projection could make African countries seem smaller than they are, creating sort of a sense of unimportant: something that many Westerners subconsciously feel about South Africa.

3 Map 3: Redesign a Past Map

3.1 Context

This map was made to show the # of licenses for a study area (SA) in East Passyunk neighborhood in Philadelphia. In fact, the dotted red outline is the distance in which my team and I canvassed the neighborhood. Since the class was related to nighttime culture, policy, and amenities, we wanted to show both a) active business licenses in our East Passyunnk SA, and b) # of active Pennsylvania Liquor Control Board (PLCB) licenses. Restauarnts, bars, and pubs must be in possession of a PLCB license to serve and distribute alcohol.

Let’s dive in to what we are working with.

3.2 Before

Original map submitted for [CPLN 6890: The City at Night]

3.3 After

# 1. Load data from old RMD

bl <- st_read("data/business_licenses.geojson", quiet = TRUE) %>%  #Business Licenses
  st_transform(2272)

ep_boundary <- st_read("data/study_areas_2025_1.shp", quiet = TRUE) %>% # Custom study area boundary
  st_transform(2272)

plcb_raw <- readr::read_csv("data/PHL_PLCB_geocoded.csv")

plcb_sf <- plcb_raw %>%
  mutate(lon = as.numeric(lon), lat = as.numeric(lat)) %>%
  filter(!is.na(lon), !is.na(lat)) %>%
  st_as_sf(coords = c("lon","lat"), crs = 4326) %>%
  st_transform(2272)

# 2. Clip to study area
bl_ep    <- st_intersection(bl, ep_boundary)
plcb_ep  <- st_intersection(plcb_sf, ep_boundary)

# 3. Aggregate data more cleanly than last semester
bl_ep_clean <- bl_ep %>%
  filter(licensestatus == "Active") %>%
  mutate(category = case_when(
    licensetype %in% c(
      "Food Preparing and Serving",
      "Food Preparing and Serving (30+ SEATS)",
      "Sidewalk Cafe",
      "Food Establishment, Retail Permanent Location",
      "Food Caterer",
      "Food Estab, Retail Non-Permanent Location (Annual)"
    ) ~ "Food Service",
    licensetype %in% c(
      "Amusement",
      "Special Assembly Occupancy"
    ) ~ "Entertainment",
    TRUE ~ NA_character_
  )) %>%
  filter(!is.na(category)) %>%
  arrange(addressobjectid, category) %>%
  group_by(addressobjectid) %>%
  slice(1) %>%
  ungroup() %>%
  st_zm(drop = TRUE)

# verify immediately after
nrow(bl_ep_clean)  # should be 48
table(bl_ep_clean$category)

plcb_ep_clean <- plcb_ep %>%
  filter(Status == "Active") %>%
  mutate(category = "Liquor License") %>%
  st_zm(drop = TRUE)

# 4. Combine East Passyunk Business Licenses and PLCB East Passyunk licenses together
licenses_clean <- bind_rows(
  bl_ep_clean %>%
    transmute(category = category, name = business_name, geometry = geometry),
  plcb_ep_clean %>%
    transmute(category = category, name = Licensee, geometry = geometry)
)
# ggplot draws factors in order, so Restaurant renders on top

# 5. Grab roads via tigris, clip to SA boundary
options(tigris_use_cache = TRUE)
roads <- roads(state = "PA", county = "Philadelphia") %>%
  st_transform(2272)

ep_bbox <- st_bbox(ep_boundary)

Things to note:

I aggregated the active business licenses types a bit differently than I should have for class, so I made a few changes here.
Anything related to Food Preparing and Serving to Food Establishment and Sidewalk Cafe will fall under “Food Service”. The Amusement and Special Assembly Occupancy is thrust under “Entertainment”. The other licenses will be labeled as Other and then dropped entirely from mapping.

# philly boundary for inset
philly <- counties(state = "PA", cb = TRUE) %>%
  st_transform(2272) %>%
  filter(NAME == "Philadelphia")

# Load neighborhoods
neighborhoods_url <- "https://raw.githubusercontent.com/opendataphilly/open-geo-data/master/philadelphia-neighborhoods/philadelphia-neighborhoods.geojson"
neighborhoods <- st_read(neighborhoods_url) %>%
  st_transform(2272)

# get the neighborhood polygons you want first
ep_neighborhoods <- neighborhoods %>%
  st_transform(2272) %>%
  filter(NAME %in% c("PASSYUNK_SQUARE", "EAST_PASSYUNK")) %>%
  st_union()

ep_highlight <- neighborhoods %>%
  filter(NAME == "EAST_PASSYUNK")

# Roads for context
roads_context <- st_intersection(
  roads %>% filter(RTTYP %in% c("M","S")),
  philly  
) 


roads_aoi <- st_intersection(
  roads %>% filter(RTTYP %in% c("M", "S")),
  ep_neighborhoods
)

# Rivers
philly_rivers <- ne_download(scale = 10, type = "rivers_lake_centerlines",
                             category = "physical", returnclass = "sf") %>%
  st_transform(2272) %>%
  st_crop(st_bbox(philly))

# Neighboring Counties
neighbors <- counties(state = c("PA", "NJ"), cb = TRUE) %>%
  st_transform(2272) %>%
  filter(NAME %in% c("Delaware", "Montgomery", "Bucks", "Camden", 
                     "Gloucester", "Burlington"))

# Title Count
n_liquor     <- sum(licenses_clean$category == "Liquor License")
n_restaurant <- sum(licenses_clean$category == "Food Service")
n_entertain  <- sum(licenses_clean$category == "Entertainment")

map_title <- glue(
  "{n_liquor} Liquor Licenses and {n_restaurant} Food Establishments\nanchor East Passyunk's Corridor Identity"
)

# Palette almost replicates the East Passyunk Business Improvement District
category_palette <- c(
  "Food Service"    = "#B000EC",
  "Liquor License"  = "#E87301",
  "Entertainment"   = "#7a6652"
)

# ── main map ──────────────────────────────────────────────────────────────────
main_map <- ggplot() +
  geom_sf(data = roads_aoi, color = "#BBBBAA", linewidth = 0.25) +
  geom_sf(data = ep_boundary, fill = NA, color = "#060606",
          linewidth = 0.5, linetype = "longdash") +
  geom_sf(data = licenses_clean, aes(color = category),
          size = 2, alpha = 0.85) +
  scale_color_manual(values = category_palette) +
  coord_sf(
    xlim = c(ep_bbox["xmin"] - 50, ep_bbox["xmax"] + 50),
    ylim = c(ep_bbox["ymin"] - 50, ep_bbox["ymax"] + 50)
  ) +
  labs(color = NULL, caption = "Data: OpenDataPhilly; PLCB via M. Fichman; Tigris") +
  theme_void() +
  theme(
    plot.background  = element_rect(fill = "#EBEBDF", color = NA),
    panel.background = element_rect(fill = "#EBEBDF", color = NA),
    plot.caption  = element_text(size = 7, color = "#999988", margin = margin(t=8)),
    plot.margin   = margin(10, 10, 10, 10),
    legend.position  = "bottom",
    legend.text      = element_text(size = 9, color = "#333333")
  ) +
  annotation_north_arrow(location = "tr",
    style = north_arrow_minimal(),
    height = unit(0.6,"cm"), width = unit(0.6,"cm")) +
  annotation_scale(location = "br", text_col = "#666655", line_col = "#666655")

inset_map <- ggplot() +
  geom_sf(data = philly, fill = "#DEDED4", color = "#BBBBAA", linewidth = 0.4) +
  geom_sf(data = neighborhoods %>% st_transform(2272),
          fill = NA, color = "#CCCCBB", linewidth = 0.1) +
  geom_sf(data = philly_rivers, color = "#A8C8E0", linewidth = 0.6) +
  geom_sf(data = ep_highlight, fill = "#B000EC", color = "#B000EC",
          linewidth = 1, alpha = 0.8) +
 annotate("text", 
           x = st_coordinates(st_centroid(ep_highlight))[1] + 8000,
           y = st_coordinates(st_centroid(ep_highlight))[2] - 4000,
           label = "East Passyunk", 
           size = 2.2, color = "#EF1D1D", fontface = "italic",
           hjust = 0) +
  theme_void() +
  theme(
    plot.background  = element_rect(fill = "#EBEBDF",
                                    color = "#111111",  # bold black border
                                    linewidth = 1.2),
    panel.background = element_rect(fill = "#EBEBDF", color = NA)
  )

## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries

combined_map <- ggdraw() +
  draw_plot(main_map) +
  draw_plot(inset_map,
            x = 0.68, y = 0.67,   # pushed right, slightly lower
            width = 0.30, height = 0.30) +
  draw_label(map_title,
             x = 0.05, y = 0.99,
             hjust = 0, vjust = 1,
             size = 13, fontface = "bold", color = "#222222",
             lineheight = 1.2) +
  draw_label("Active business and PLCB liquor licenses, East Passyunk AOI",
             x = 0.05, y = 0.88,   # more gap from title
             hjust = 0, vjust = 1,
             size = 9, color = "#666655")
# save combined


#ggsave2("ep_combined_map.png", plot = combined_map,
        #width = 8, height = 9, dpi = 300)

combined_map

3.4 Reflection

One of the stark problems with my original map is the usage of the dark/black background. I wanted to try and emulate the course theme, but now I am aware it is not the greatest styling technique. My takeaway title was very lack-luster and does not give the reader much to go off of. Furthermore, the original map contained too many obsolete licenses that the OG map gives you a heart attack trying to interpret big cirlce from small circle. The legend did not offer much help in this case.

Here are some of the changes I chose to make:

White + tan background sort of emulates a CARTO/ArcGIS grey background, but I personally like the almost greenish reflection.
Dropped meaningless categories when filtering license types. All food service types go in one bin, active liquor types go in another. If people are really curious about the various types of liquor licenses, take a look at the code.
Developed an inset map for spatial context. If someone looked at this map w/o the title present, they would have no clue where we are. If I had more time, I would annotate where the Acme and names of various establishments on the primary map for further world building.

Big Idea Statement:

East Passyunk’s identity as a dining and nightlife destination is supported by its license landscape; 21 liquor licenses and 46 food establishments are composed in about a 0.7 mile stretch.

3.4.1 AI Declaration

I used Claude for formatting recommendations and with questions regarding some calculations. I also asked Cortex, Snowflake’s native AI for South African data visualization tips, since this data is for a different project.

MUSA 6310 Map Upgrades

Alex Stauffer

2026-04-04