Miniproject 5 — Japan at Night (VIIRS 2025) Julia Fritsch DST 234

Code of Ethics (required)

In order to complete this assignment, I reviewed previous in-class examples, materials, and homework for reference. In addition, I used GitHub Copilot to help me think/navigate through how to approach the project based on my idea. I first applied the coding concepts I had already learned in class/HW, then used Copilot for suggestions on revisions and troubleshooting errors. I also looked up some examples of other leaflets to get inspriration. I wanted to expand in a different way from my previous project and did a lot of outside research into different leaflet code/formats. I originally was going to use Shiny(more interactive) for this project, but had a lot of compiling errors that were time costly.

To go beyond the basic requirements, I found a the Annual VNL V2 which was perfect for my idea it is a new consistently processed time series of annual global VIIRS nighttime lights that have been produced from monthly cloud-free average radiance grids spanning 2012* to 2020. The new methodology is a modification of the original method based on nightly data (Annual VNL V1). https://eogdata.mines.edu/products/vnl/ This was interesting to me and I had wanted to spent more time expanding on the features of each leaflet; but still think I did a lot with the time and resources I had learning a new way to apply code.

---

Setup

library(tidyverse)
library(terra)
library(leaflet)
library(viridisLite)
library(ggplot2)
library(htmltools)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)

# Japan bounding box (lon/lat)
japan_bbox <- ext(122.5, 146.5, 24.0, 46.5)

# City markers for context
cities_jp <- tibble::tribble(
  ~city, ~lat, ~lng,
  "Sapporo", 43.0618, 141.3545,
  "Sendai", 38.2682, 140.8694,
  "Tokyo", 35.6762, 139.6503,
  "Nagoya", 35.1815, 136.9066,
  "Osaka", 34.6937, 135.5023,
  "Hiroshima", 34.3853, 132.4553,
  "Fukuoka", 33.5904, 130.4017,
  "Naha", 26.2124, 127.6809
)

cities_sf <- sf::st_as_sf(cities_jp, coords = c("lng", "lat"), crs = 4326)

---

This report uses small Japan-only rasters in data/ for speed purposes.

dir.create("data", showWarnings = FALSE)

in_brightness <- "data-raw/VNL_npp_2025_global_vcmslcfg_v2_c202604011200.median_masked.dat.tif"
in_litmask    <- "data-raw/VNL_npp_2025_global_vcmslcfg_v2_c202604011200.lit_mask.dat.tif"

out_brightness <- "data/viirs2025_japan_brightness_small.tif"
out_litmask    <- "data/viirs2025_japan_litmask_small.tif"

need_build <- !file.exists(out_brightness) || !file.exists(out_litmask)

if (need_build) {
  if (!file.exists(in_brightness) || !file.exists(in_litmask)) {
    stop(
      "Missing files.\n\n",
      "Provide either:\n",
      "A) raw inputs in data-raw/:\n",
      "   - ", in_brightness, "\n",
      "   - ", in_litmask, "\n\n",
      "OR\n",
      "B) processed outputs in data/:\n",
      "   - ", out_brightness, "\n",
      "   - ", out_litmask, "\n"
    )
  }

  message("Building Japan-only rasters (one-time). This may take a few minutes...")

  r_b <- rast(in_brightness) |> clamp(lower = 0, values = TRUE)
  r_m <- rast(in_litmask)

  r_b_jp <- crop(r_b, japan_bbox)
  r_m_jp <- crop(r_m, japan_bbox)

  # bigger = faster but less detail
  fact_speed <- 10
  r_b_small <- aggregate(r_b_jp, fact = fact_speed, fun = "mean", na.rm = TRUE)
  r_m_small <- aggregate(r_m_jp, fact = fact_speed, fun = "modal", na.rm = TRUE)

  # Trim outliers for nicer color scaling
  q <- quantile(values(r_b_small), probs = c(0.01, 0.995), na.rm = TRUE)
  r_b_small <- clamp(r_b_small, lower = q[1], upper = q[2], values = TRUE)

  writeRaster(r_b_small, out_brightness, overwrite = TRUE)
  writeRaster(r_m_small, out_litmask, overwrite = TRUE)

  message("Wrote:\n- ", out_brightness, "\n- ", out_litmask)
}

stopifnot(file.exists(out_brightness), file.exists(out_litmask))

r_b0 <- rast(out_brightness)
r_m0 <- rast(out_litmask)

---

Summary metrics

# Lit % (pixel-based)
m_vals <- values(r_m0)
m_vals <- m_vals[is.finite(m_vals)]
pct_lit <- mean(m_vals > 0, na.rm = TRUE) * 100

# Brightness summaries (raw + log)
b_raw <- values(r_b0)
b_raw <- b_raw[is.finite(b_raw)]
b_log <- log10(b_raw + 1)

summary_tbl <- tibble(
  metric = c("Lit % (mask > 0)",
             "Brightness median (raw)", "Brightness 90th (raw)", "Brightness 99th (raw)",
             "Brightness median (log10+1)", "Brightness 90th (log10+1)", "Brightness 99th (log10+1)"),
  value = c(
    sprintf("%.1f%%", pct_lit),
    sprintf("%.4f", quantile(b_raw, 0.50, na.rm = TRUE)),
    sprintf("%.4f", quantile(b_raw, 0.90, na.rm = TRUE)),
    sprintf("%.4f", quantile(b_raw, 0.99, na.rm = TRUE)),
    sprintf("%.4f", quantile(b_log, 0.50, na.rm = TRUE)),
    sprintf("%.4f", quantile(b_log, 0.90, na.rm = TRUE)),
    sprintf("%.4f", quantile(b_log, 0.99, na.rm = TRUE))
  )
)

summary_tbl

# A tibble: 7 × 2
  metric                      value 
  <chr>                       <chr> 
1 Lit % (mask > 0)            5.9%  
2 Brightness median (raw)     0.0000
3 Brightness 90th (raw)       0.0646
4 Brightness 99th (raw)       5.1349
5 Brightness median (log10+1) 0.0000
6 Brightness 90th (log10+1)   0.0272
7 Brightness 99th (log10+1)   0.7878

4) Geographic context (coastlines/borders)

world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")

bbox_sf <- sf::st_as_sfc(sf::st_bbox(
  c(xmin = 122.5, xmax = 146.5, ymin = 24.0, ymax = 46.5),
  crs = 4326
))

coast_context <- sf::st_intersection(world, bbox_sf)

# Simplify for faster leaflet rendering
coast_context_s <- sf::st_simplify(coast_context, dTolerance = 0.05)

---

1) Hotspot labels (top brightest points)

rb_hot <- log10(r_b0 + 1)

pts <- terra::as.points(rb_hot, values = TRUE, na.rm = TRUE)

# Convert to DF with coordinates + value
pts_df <- as.data.frame(pts, xy = TRUE)

# ---- robust coordinate columns ----
cn <- names(pts_df)
x_col <- cn[1]
y_col <- cn[2]

pts_df <- pts_df |>
  dplyr::rename(lng = !!x_col, lat = !!y_col)

# ---- robust value column ----
val_col <- names(pts_df)[ncol(pts_df)]
pts_df <- pts_df |>
  dplyr::mutate(b_log = as.numeric(.data[[val_col]])) |>
  dplyr::filter(is.finite(b_log))

top_n <- 20

if (nrow(pts_df) == 0) {
  hotspots_sf <- sf::st_as_sf(
    tibble(lng = numeric(), lat = numeric(), b_log = numeric(), rank = integer(), label = character()),
    coords = c("lng", "lat"),
    crs = 4326
  )
} else {
  o <- order(pts_df$b_log, decreasing = TRUE)
  o <- o[seq_len(min(top_n, length(o)))]
  hotspots <- pts_df[o, c("lng", "lat", "b_log"), drop = FALSE]

  hotspots$rank <- seq_len(nrow(hotspots))
  hotspots$label <- paste0(
    "Hotspot #", hotspots$rank, " (", round(hotspots$b_log, 2), ")"
  )

  hotspots_sf <- sf::st_as_sf(hotspots, coords = c("lng", "lat"), crs = 4326)
}

---

Helper: build a styled leaflet map (used by multiple sections)

badge <- function(txt, bg = "#111827") {
  tags$span(
    style = paste0(
      "display:inline-block; padding:6px 10px; margin-right:8px; margin-top:6px;",
      "border-radius:999px; font-size:12px; font-weight:700;",
      "background:", bg, "; color:white; border:1px solid rgba(255,255,255,0.15);"
    ),
    txt
  )
}

panel_html <- function(title, subtitle) {
  tags$div(
    style = paste0(
      "background: rgba(0,0,0,0.72); color: white; padding: 14px 14px; ",
      "border-radius: 14px; max-width: 560px; ",
      "box-shadow: 0 12px 30px rgba(0,0,0,0.35);"
    ),
    tags$div(style = "font-size:18px; font-weight:900; margin-bottom:6px;", title),
    tags$div(style = "font-size:13px; line-height:1.35; opacity:0.95;", subtitle),
    tags$div(
      badge(paste0("Lit area: ", sprintf("%.1f%%", pct_lit)), "#0b1220"),
      badge(paste0("Median (log): ", sprintf("%.3f", quantile(b_log, 0.50, na.rm = TRUE))), "#111827"),
      badge(paste0("99th% (log): ", sprintf("%.3f", quantile(b_log, 0.99, na.rm = TRUE))), "#1f2937")
    ),
    tags$div(
      style = "margin-top:8px; font-size:11px; opacity:0.85;",
      "Ethics note: radiance ≠ skyglow; do not infer causation from brightness patterns."
    )
  )
}

make_map <- function(rb, vals, pal_b, title, subtitle,
                     show_mask = TRUE, show_hotspots = TRUE, show_cities = TRUE) {

  # Lit mask palette
  m_vals2 <- values(r_m0)
  m_vals2 <- m_vals2[is.finite(m_vals2)]
  pal_m <- colorBin(
    palette = c("transparent", "#2ee9ff"),
    domain = m_vals2,
    bins = c(-Inf, 0, Inf),
    na.color = "transparent"
  )

  base_maps <- c(
    "Dark"      = "CartoDB.DarkMatter",
    "Light"     = "CartoDB.Positron",
    "Satellite" = "Esri.WorldImagery",
    "Topo"      = "Esri.WorldTopoMap"
  )

  m <- leaflet(options = leafletOptions(zoomControl = TRUE, preferCanvas = TRUE)) |>
    addProviderTiles(base_maps[["Dark"]], group = "Dark") |>
    addProviderTiles(base_maps[["Light"]], group = "Light") |>
    addProviderTiles(base_maps[["Satellite"]], group = "Satellite") |>
    addProviderTiles(base_maps[["Topo"]], group = "Topo") |>
    setView(lng = 138.2, lat = 36.2, zoom = 5) |>

    addPolygons(
      data = coast_context_s,
      color = "rgba(255,255,255,0.25)",
      weight = 1,
      fill = FALSE,
      group = "Coastlines"
    ) |>

    addRasterImage(
      rb,
      colors = pal_b,
      opacity = 0.88,
      project = TRUE,
      group = "Brightness"
    ) |>
    addLegend(
      pal = pal_b,
      values = vals,
      title = "Brightness",
      position = "bottomright",
      opacity = 1
    ) |>

    addScaleBar(position = "bottomleft") |>
    addControl(html = panel_html(title, subtitle), position = "topleft")

  if (show_mask) {
    m <- m |>
      addRasterImage(
        r_m0,
        colors = pal_m,
        opacity = 0.28,
        project = TRUE,
        group = "Lit mask"
      )
  }

  if (show_cities) {
    m <- m |>
      addCircleMarkers(
        data = cities_sf,
        radius = 6,
        stroke = TRUE, weight = 1,
        color = "rgba(255,255,255,0.95)",
        fillColor = "#ffd166",
        fillOpacity = 0.95,
        label = ~city,
        group = "Cities",
        clusterOptions = markerClusterOptions(
          showCoverageOnHover = FALSE,
          spiderfyOnMaxZoom = TRUE,
          zoomToBoundsOnClick = TRUE
        )
      )
  }

  if (show_hotspots && exists("hotspots_sf") && nrow(hotspots_sf) > 0) {
    m <- m |>
      addCircleMarkers(
        data = hotspots_sf,
        radius = 5,
        stroke = TRUE, weight = 1,
        color = "#ffffff",
        fillColor = "#ff4d6d",
        fillOpacity = 0.85,
        label = ~label,
        group = "Hotspots"
      )
  }

  m |>
    addLayersControl(
      baseGroups = names(base_maps),
      overlayGroups = c("Brightness", "Lit mask", "Cities", "Hotspots", "Coastlines"),
      options = layersControlOptions(collapsed = FALSE)
    )
}

---

2) Side-by-side maps: Raw vs Log (no sync dependency)

# Raw brightness map palette
rb_raw <- r_b0
vals_raw <- values(rb_raw)
vals_raw <- vals_raw[is.finite(vals_raw)]
nbins <- 16

breaks_raw <- unique(quantile(vals_raw, probs = seq(0, 1, length.out = nbins + 1), na.rm = TRUE))
pal_raw <- colorBin(
  palette = viridisLite::cividis(nbins),
  domain = vals_raw,
  bins = breaks_raw,
  na.color = "transparent"
)

# Log brightness map palette
rb_log <- log10(r_b0 + 1)
vals_log <- values(rb_log)
vals_log <- vals_log[is.finite(vals_log)]
breaks_log <- unique(quantile(vals_log, probs = seq(0, 1, length.out = nbins + 1), na.rm = TRUE))
pal_log <- colorBin(
  palette = viridisLite::magma(nbins),
  domain = vals_log,
  bins = breaks_log,
  na.color = "transparent"
)

m_raw <- make_map(
  rb = rb_raw, vals = vals_raw, pal_b = pal_raw,
  title = "Raw Radiance (Quantile-binned)",
  subtitle = "Higher values dominate; dim structure can be hard to see.",
  show_mask = TRUE, show_hotspots = TRUE, show_cities = TRUE
)

m_log <- make_map(
  rb = rb_log, vals = vals_log, pal_b = pal_log,
  title = "Log Brightness: log10(radiance + 1)",
  subtitle = "Reveals mid-level structure beyond major metros.",
  show_mask = TRUE, show_hotspots = TRUE, show_cities = TRUE
)

htmltools::browsable(
  htmltools::div(
    style = "display:flex; gap:12px; align-items:stretch; flex-wrap:wrap;",
    htmltools::div(
      style = "flex:1; min-width:380px; height:650px;",
      m_raw
    ),
    htmltools::div(
      style = "flex:1; min-width:380px; height:650px;",
      m_log
    )
  )
)

---

3) Animated story (chapter maps)

These are 3 “chapters” with different presets. Scroll through them like a story.

Chapter 1 — The big metros (high intensity)

make_map(
  rb = rb_log, vals = vals_log, pal_b = pal_log,
  title = "Chapter 1: The Big Metros",
  subtitle = "Tokyo–Nagoya–Osaka corridor dominates the radiance distribution.",
  show_mask = FALSE, show_hotspots = TRUE, show_cities = TRUE
)

Chapter 2 — The lit mask (where stable lighting exists)

make_map(
  rb = rb_log, vals = vals_log, pal_b = pal_log,
  title = "Chapter 2: Where Lights Exist",
  subtitle = "Lit mask overlay highlights stable lighting detection even when intensity is modest.",
  show_mask = TRUE, show_hotspots = FALSE, show_cities = FALSE
)

Chapter 3 — Context + infrastructure edges

make_map(
  rb = rb_log, vals = vals_log, pal_b = pal_log,
  title = "Chapter 3: Coastlines & Corridors",
  subtitle = "Coastal infrastructure and corridor patterns appear when intensity is log-scaled.",
  show_mask = TRUE, show_hotspots = TRUE, show_cities = TRUE
)

---

Original graphic (static; include for rubric)

ggplot(tibble(radiance = b_raw), aes(radiance)) +
  geom_histogram(bins = 70, fill = "gray70", color = "white") +
  scale_x_continuous(trans = "log10") +
  labs(
    title = "Original (static): Distribution of VIIRS radiance in Japan (2025)",
    x = "Radiance (log10 scale)",
    y = "Pixel count",
    caption = "Static distribution is informative but not spatial; revised maps add geography, context, and exploration."
  ) +
  theme_minimal(base_size = 13)

---

Bonus: Styled distribution (matches the map better :D)

ggplot(tibble(brightness_log = b_log), aes(brightness_log)) +
  geom_histogram(bins = 70, fill = viridisLite::magma(1), color = "white", alpha = 0.95) +
  labs(
    title = "Brightness distribution (Japan bbox)",
    subtitle = "Transformed as log10(radiance + 1)",
    x = "log10(radiance + 1)",
    y = "Pixel count"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

---

What I changed (explanation)

• Added hotspot labeling (top-brightness pixels) to guide exploration.
• Added side-by-side maps to compare raw vs log brightness (no extra packages needed).
• Added a scroll narrative (“chapters”) using multiple curated map views.
• Added geographic context (coastlines/borders) using Natural Earth.
• Kept ethical framing and clarified that radiance is not the same as skyglow.

---

Data source

Earth Observation Group (EOG), VIIRS Nighttime Lights (Annual VNL V2, 2025):
https://eogdata.mines.edu/products/vnl/