day-25-hexagons.R

# Day 25: Hexagons - US Wind Energy Infrastructure
# Multi-resolution H3 hexagonal aggregation of 76,000+ wind turbines
# Aggregating total installed capacity (MW) per hexagon
# Points fade in at high zoom for individual turbine exploration

library(mapgl)
library(h3o)
library(sf)
library(dplyr)
library(readr)

# Load USGS Wind Turbine Database (76,000+ turbines with precise GPS coordinates)
# Source: https://energy.usgs.gov/uswtdb/data/
turbines_raw <- read_csv("uswtdb_V8_1_20250522.csv", show_col_types = FALSE)

# Filter to continental US with valid coordinates and capacity
turbines <- turbines_raw |>
  filter(
    !is.na(xlong),
    !is.na(ylat),
    !is.na(t_cap)
  ) |>
  mutate(
    capacity_mw = t_cap / 1000 # Convert kW to MW
  ) |>
  select(
    lon = xlong,
    lat = ylat,
    capacity_mw,
    manufacturer = t_manu,
    model = t_model,
    project = p_name,
    year = p_year,
    state = t_state,
    hub_height = t_hh,
    rotor_dia = t_rd
  )

cat("Turbines loaded:", nrow(turbines), "\n")
cat("Total capacity:", round(sum(turbines$capacity_mw), 1), "MW\n")

# Convert to sf
turbines_sf <- st_as_sf(turbines, coords = c("lon", "lat"), crs = 4326)

# Inspect
maplibre_view(turbines_sf)

# Generate H3 indices at multiple resolutions
turbines_h3 <- turbines_sf |>
  mutate(
    h3_res3 = h3_from_points(geometry, 3),
    h3_res4 = h3_from_points(geometry, 4),
    h3_res5 = h3_from_points(geometry, 5),
    h3_res6 = h3_from_points(geometry, 6)
  )

# Aggregation function - sum capacity (MW) per hexagon
aggregate_h3 <- function(data, h3_col) {
  agg <- data |>
    st_drop_geometry() |>
    mutate(h3_str = as.character(.data[[h3_col]])) |>
    filter(!is.na(h3_str)) |>
    group_by(h3_str) |>
    summarise(
      total_mw = round(sum(capacity_mw, na.rm = TRUE), 1),
      turbine_count = n(),
      .groups = "drop"
    )

  agg |>
    mutate(
      h3 = h3_from_strings(h3_str),
      geometry = st_as_sfc(h3)
    ) |>
    select(total_mw, turbine_count, geometry) |>
    st_as_sf(crs = 4326)
}

hex_res3 <- aggregate_h3(turbines_h3, "h3_res3")
hex_res4 <- aggregate_h3(turbines_h3, "h3_res4")
hex_res5 <- aggregate_h3(turbines_h3, "h3_res5")
hex_res6 <- aggregate_h3(turbines_h3, "h3_res6")

cat("\nHexagon counts by resolution:\n")
cat(
  "Res 3:",
  nrow(hex_res3),
  "hexagons, MW range:",
  round(range(hex_res3$total_mw)),
  "\n"
)
cat(
  "Res 4:",
  nrow(hex_res4),
  "hexagons, MW range:",
  round(range(hex_res4$total_mw)),
  "\n"
)
cat(
  "Res 5:",
  nrow(hex_res5),
  "hexagons, MW range:",
  round(range(hex_res5$total_mw)),
  "\n"
)
cat(
  "Res 6:",
  nrow(hex_res6),
  "hexagons, MW range:",
  round(range(hex_res6$total_mw)),
  "\n"
)

# Wind energy color palette - clean, saturated for light basemap
wind_colors <- c(
  "#e0f3db",
  "#a8ddb5",
  "#4eb3d3",
  "#2b8cbe",
  "#0868ac",
  "#084081"
)

# Create color expressions based on MW capacity
make_color_expr <- function(data) {
  max_val <- max(data$total_mw, na.rm = TRUE)
  interpolate(
    column = "total_mw",
    values = c(
      0,
      max_val * 0.05,
      max_val * 0.15,
      max_val * 0.35,
      max_val * 0.6,
      max_val
    ),
    stops = wind_colors
  )
}

color_expr3 <- make_color_expr(hex_res3)
color_expr4 <- make_color_expr(hex_res4)
color_expr5 <- make_color_expr(hex_res5)
color_expr6 <- make_color_expr(hex_res6)

# Build the multi-resolution hexagon map
maplibre(
  style = maptiler_style("dataviz", variant = "light"),
  center = c(-98.5, 39.5),
  zoom = 4,
  pitch = 40,
  bearing = 0
) |>
  # Resolution 3: Continental overview (zoom 0-5)
  add_fill_layer(
    id = "hex-res3",
    source = hex_res3,
    fill_color = color_expr3,
    fill_opacity = interpolate(
      property = "zoom",
      type = list("exponential", 1.5),
      values = c(0, 4.5, 5.5),
      stops = c(0.9, 0.9, 0)
    ),
    fill_outline_color = "#08408180",
    max_zoom = 6,
    tooltip = concat(
      number_format("turbine_count", use_grouping = TRUE),
      " turbines | ",
      number_format(
        "total_mw",
        maximum_fraction_digits = 1,
        use_grouping = TRUE
      ),
      " MW"
    )
  ) |>
  # Resolution 4: Regional patterns (zoom 4-7)
  add_fill_layer(
    id = "hex-res4",
    source = hex_res4,
    fill_color = color_expr4,
    fill_opacity = interpolate(
      property = "zoom",
      type = list("exponential", 1.5),
      values = c(4.5, 5.5, 6.5, 7.5),
      stops = c(0, 0.9, 0.9, 0)
    ),
    fill_outline_color = "#08408160",
    min_zoom = 4.5,
    max_zoom = 8,
    tooltip = concat(
      number_format("turbine_count", use_grouping = TRUE),
      " turbines | ",
      number_format(
        "total_mw",
        maximum_fraction_digits = 1,
        use_grouping = TRUE
      ),
      " MW"
    )
  ) |>
  # Resolution 5: State/regional detail (zoom 6-9)
  add_fill_layer(
    id = "hex-res5",
    source = hex_res5,
    fill_color = color_expr5,
    fill_opacity = interpolate(
      property = "zoom",
      type = list("exponential", 1.5),
      values = c(6.5, 7.5, 8.5, 9.5),
      stops = c(0, 0.9, 0.9, 0)
    ),
    fill_outline_color = "#08408150",
    min_zoom = 6.5,
    max_zoom = 10,
    tooltip = concat(
      number_format("turbine_count", use_grouping = TRUE),
      " turbines | ",
      number_format(
        "total_mw",
        maximum_fraction_digits = 1,
        use_grouping = TRUE
      ),
      " MW"
    )
  ) |>
  # Resolution 6: Wind farm detail (zoom 8-11)
  add_fill_layer(
    id = "hex-res6",
    source = hex_res6,
    fill_color = color_expr6,
    fill_opacity = interpolate(
      property = "zoom",
      type = list("exponential", 1.5),
      values = c(8.5, 9.5, 10.5, 11.5),
      stops = c(0, 0.9, 0.9, 0)
    ),
    fill_outline_color = "#08408140",
    min_zoom = 8.5,
    max_zoom = 12,
    tooltip = concat(
      number_format("turbine_count", use_grouping = TRUE),
      " turbines | ",
      number_format(
        "total_mw",
        maximum_fraction_digits = 1,
        use_grouping = TRUE
      ),
      " MW"
    )
  ) |>
  # Individual turbine points - fade in past zoom 9
  add_circle_layer(
    id = "turbine-glow",
    source = turbines_sf,
    circle_radius = interpolate(
      property = "zoom",
      values = c(9, 12, 16),
      stops = c(6, 12, 20)
    ),
    circle_color = "#0868ac",
    circle_blur = 1,
    circle_opacity = interpolate(
      property = "zoom",
      values = c(9, 10.5, 12),
      stops = c(0, 0.15, 0.25)
    ),
    min_zoom = 9
  ) |>
  add_circle_layer(
    id = "turbine-points",
    source = turbines_sf,
    circle_radius = interpolate(
      property = "zoom",
      values = c(9, 12, 16),
      stops = c(3, 6, 12)
    ),
    circle_color = interpolate(
      column = "capacity_mw",
      values = c(0, 2, 4, 8, 13),
      stops = c("#a8ddb5", "#4eb3d3", "#2b8cbe", "#0868ac", "#084081")
    ),
    circle_opacity = interpolate(
      property = "zoom",
      values = c(9, 10.5, 12),
      stops = c(0, 0.85, 1)
    ),
    circle_stroke_color = "#ffffff",
    circle_stroke_width = interpolate(
      property = "zoom",
      values = c(9, 14),
      stops = c(1, 2)
    ),
    circle_stroke_opacity = interpolate(
      property = "zoom",
      values = c(9, 11, 13),
      stops = c(0, 0.7, 1)
    ),
    min_zoom = 9,
    popup = concat(
      "<div style='font-family: -apple-system, BlinkMacSystemFont, \"Segoe UI\", Roboto, sans-serif; ",
      "padding: 16px; background: #fff; border-radius: 8px; color: #1a1a1a; ",
      "min-width: 260px; box-shadow: 0 4px 16px rgba(0,0,0,0.12);'>",
      "<div style='color: #0868ac; font-size: 11px; font-weight: 600; ",
      "margin-bottom: 4px; text-transform: uppercase; letter-spacing: 1.5px;'>",
      "Wind Turbine</div>",
      "<div style='font-size: 18px; font-weight: 700; color: #084081; margin-bottom: 12px; ",
      "padding-bottom: 12px; border-bottom: 2px solid #e0f3db;'>",
      get_column("project"),
      "</div>",
      "<div style='display: grid; grid-template-columns: 1fr 1fr; gap: 12px;'>",
      "<div><span style='color: #666; font-size: 10px; text-transform: uppercase; letter-spacing: 0.5px;'>Capacity</span><br>",
      "<span style='font-size: 16px; font-weight: 700; color: #0868ac;'>",
      get_column("capacity_mw"),
      " MW</span></div>",
      "<div><span style='color: #666; font-size: 10px; text-transform: uppercase; letter-spacing: 0.5px;'>Year</span><br>",
      "<span style='font-size: 16px; font-weight: 600; color: #333;'>",
      get_column("year"),
      "</span></div>",
      "<div><span style='color: #666; font-size: 10px; text-transform: uppercase; letter-spacing: 0.5px;'>Manufacturer</span><br>",
      "<span style='font-size: 14px; color: #333;'>",
      get_column("manufacturer"),
      "</span></div>",
      "<div><span style='color: #666; font-size: 10px; text-transform: uppercase; letter-spacing: 0.5px;'>Model</span><br>",
      "<span style='font-size: 14px; color: #333;'>",
      get_column("model"),
      "</span></div>",
      "<div><span style='color: #666; font-size: 10px; text-transform: uppercase; letter-spacing: 0.5px;'>Hub Height</span><br>",
      "<span style='font-size: 14px; color: #333;'>",
      get_column("hub_height"),
      " m</span></div>",
      "<div><span style='color: #666; font-size: 10px; text-transform: uppercase; letter-spacing: 0.5px;'>Rotor Diameter</span><br>",
      "<span style='font-size: 14px; color: #333;'>",
      get_column("rotor_dia"),
      " m</span></div>",
      "</div></div>"
    )
  ) |>

  add_control(
    html = paste0(
      "<div style='font-family: -apple-system, BlinkMacSystemFont, \"Segoe UI\", Roboto, sans-serif; ",
      "background: rgba(255,255,255,0.95); padding: 16px 20px; border-radius: 8px; ",
      "box-shadow: 0 2px 12px rgba(0,0,0,0.1); max-width: 320px;'>",
      "<div style='display: flex; align-items: center; gap: 12px; margin-bottom: 8px;'>",
      "<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"32\" height=\"32\" viewBox=\"0 0 24 24\" ",
      "fill=\"none\" stroke=\"#0868ac\" stroke-width=\"1.5\" stroke-linecap=\"round\" stroke-linejoin=\"round\">",
      "<path d=\"M12 11m-2 0a2 2 0 1 0 4 0a2 2 0 1 0 -4 0\" />",
      "<path d=\"M10 11v-2.573c0 -.18 .013 -.358 .04 -.536l.716 -4.828c.064 -.597 .597 -1.063 1.244 -1.063s1.18 .466 1.244 1.063l.716 4.828c.027 .178 .04 .357 .04 .536v2.573\" />",
      "<path d=\"M13.01 9.28l2.235 1.276c.156 .09 .305 .19 .446 .3l3.836 2.911c.487 .352 .624 1.04 .3 1.596c-.325 .556 -1 .782 -1.548 .541l-4.555 -1.68a3.624 3.624 0 0 1 -.486 -.231l-2.235 -1.277\" />",
      "<path d=\"M13 12.716l-2.236 1.277a3.624 3.624 0 0 1 -.485 .23l-4.555 1.681c-.551 .241 -1.223 .015 -1.548 -.54c-.324 -.557 -.187 -1.245 .3 -1.597l3.836 -2.91a3.41 3.41 0 0 1 .446 -.3l2.235 -1.277\" />",
      "<path d=\"M7 21h10\" /><path d=\"M10 21l1 -7\" /><path d=\"M13 14l1 7\" />",
      "</svg>",
      "<div>",
      "<div style='font-size: 18px; font-weight: 700; color: #084081; line-height: 1.2;'>",
      "U.S. Wind Energy Infrastructure</div>",
      "<div style='font-size: 12px; color: #666; font-weight: 500;'>",
      "76,000+ turbines | 154 GW installed capacity</div>",
      "</div></div>",
      "<div style='font-size: 11px; color: #888; line-height: 1.4; border-top: 1px solid #e0e0e0; ",
      "padding-top: 8px; margin-top: 4px;'>",
      "Hexagons show total MW capacity per area. Zoom in to explore individual turbines.",
      "</div>",
      "<div style='font-size: 10px; color: #aaa; margin-top: 8px;'>",
      "Source: <a href='https://energy.usgs.gov/uswtdb/' target='_blank' ",
      "style='color: #0868ac; text-decoration: none;'>USGS Wind Turbine Database</a> (2025)</div>",
      "</div>"
    ),
    position = "top-left",
    id = "title-control"
  ) |>
  add_legend(
    legend_title = "Installed Capacity (MW)",
    values = c("Low", "", "", "", "", "High"),
    colors = wind_colors,
    type = "continuous",
    position = "bottom-left",
    width = "180px"
  ) |>
  add_fullscreen_control() |>
  add_navigation_control(position = "top-right")