# Create a data directory if it doesn't exist
if (!dir.exists("data")) {
dir.create("data")
}
# Download Yellowknife municipal boundary
boundary_url <- "https://opendata.yellowknife.ca/home/ServeFile/8812810a-8f07-4e83-b086-9f43e2f93aa4?FileType=7"
# Perform GET request with httr
res <- httr::GET(
url = boundary_url,
httr::write_disk("data/yellowknife_boundary.zip", overwrite = TRUE),
httr::verbose()
)Geospatial Analysis of Canadian Wildfire Impacts: A Multi-sensor Approach
Overview
In this lesson, we will analyze the 2023 Yellowknife wildfire evacuation using multiple satellite-derived datasets. The analysis integrates active fire detections, smoke plumes, wind patterns, and infrastructure data to understand the complex conditions that led to the evacuation of approximately 20,000 residents. You will learn how to access and process data from NASA’s FIRMS (Fire Information for Resource Management System), NOAA’s GOES-18 satellite, and VIIRS (Visible Infrared Imaging Radiometer Suite) burned area products. The lesson demonstrates how to combine these datasets to assess evacuation route vulnerability and regional smoke impacts.
Through hands-on analysis, you will work with both near real-time and post-processed satellite products. The integration of these datasets shows how emergency managers can use publicly available data to support evacuation decisions and understand hazard conditions. You will also learn key techniques for processing and visualizing spatial data, including clustering algorithms for fire perimeter detection, smoke plume analysis, and the creation of integrated hazard visualizations.
This lesson emphasizes practical applications of remote sensing data for emergency management, while also exploring the strengths and limitations of different satellite products. By analyzing the Yellowknife evacuation, you will gain experience with workflows that can be applied to other wildfire events and natural hazards.
Preface: A Note on Data Interpretation and Emergency Management
This lesson uses real data from an evacuation that impacted thousands of residents. While we focus on technical analysis, it’s important to remember that wildfire evacuation decisions involve many factors beyond what can be captured in satellite data. The workflows demonstrated here are meant to complement, not replace, the expertise of emergency managers, fire behavior specialists, and local authorities.
The integration of multiple data sources can provide valuable insights, but satellite data has limitations: - Cloud cover can obstruct observations - Temporal and spatial resolution may miss important changes - Data latency may affect real-time decision making - Local conditions may differ from satellite-based measurements
Always consider these limitations when using remote sensing data to support emergency management decisions.
Learning Objectives
After completing this lesson, students will be able to:
Technical Competencies
- Configure and authenticate API access to multiple NASA Earth data products
- Process and harmonize multi-source remote sensing data using R’s spatial packages (sf, terra, stars)
- Implement spatial analysis workflows combining vector and raster data for wildfire assessment
- Create publication-quality static visualizations and basic interactive maps of wildfire impacts
Analytical Skills
- Evaluate the progression and impact of wildfires using multiple remote sensing products
- Assess infrastructure vulnerability and damage using spatial analysis techniques
- Integrate multiple spatial datasets to conduct a comprehensive wildfire impact assessment
- Critical analyze the strengths and limitations of different remote sensing products for fire analysis
Practical Applications
- Apply real-world workflows for accessing and processing wildfire-related remote sensing data
- Demonstrate best practices for reproducible spatial analysis in R
- Implement methods for combining satellite-derived fire products with infrastructure data
- Document and communicate spatial analysis results effectively
Advanced Understanding (Optional Extensions)
- Develop automated workflows for multi-sensor data integration
- Apply uncertainty assessment in wildfire impact analysis
- Create advanced visualizations for temporal pattern analysis
- Implement validation approaches using multiple data sources
Introduction
Wildfires and the 2023 Season
Wildfires are a natural and essential component of North American forest ecosystems, particularly in the boreal forests that stretch across Canada. These fires play a crucial role in maintaining forest health, recycling nutrients, and promoting biodiversity. However, climate change has begun to alter historical fire regimes, leading to more frequent, intense, and extensive wildfires that pose increasing challenges to both ecosystems and human communities.
In Canada’s boreal forests, fire has historically occurred in natural cycles ranging from 50 to 200 years, helping to maintain a mosaic of forest ages and compositions. Indigenous peoples have long understood and managed these fire regimes, using controlled burns as a land management tool. Today, changing climate conditions are disrupting these traditional patterns. Longer fire seasons, increased lightning activity, and more severe drought conditions are creating unprecedented fire behavior and risk scenarios.
The human impacts of these changing fire regimes are particularly acute in Canada’s northern communities, where limited road networks often leave few evacuation options. Additionally, smoke from these fires can affect air quality across vast regions, creating public health challenges even for communities far from the actual flames. The economic implications are also significant, affecting industries from tourism to mining, and straining emergency response resources.
The 2023 wildfire season emerged as a stark illustration of these evolving challenges. The Northwest Territories experienced one of its most devastating fire seasons on record, characterized by record-breaking temperatures and severe drought conditions. The season saw over 236 fires that burned approximately 4 million hectares of land, forcing the first complete evacuation of Yellowknife in its history. The scale of displacement was unprecedented, with an estimated 70% of NWT residents forced to leave their homes. Critical infrastructure, including communities, mines, and transportation networks, faced significant threats, highlighting the vulnerability of northern communities to extreme fire events.
These events showcased both the immediate and long-term challenges of managing wildfires in a changing climate, while also demonstrating the critical importance of robust monitoring and early warning systems. The season served as a powerful example of why enhanced understanding of fire behavior and improved prediction capabilities are essential for protecting communities and ecosystems in an increasingly fire-prone world.
Yellowknife
The Northwest Territories (NWT) encompasses a vast expanse of northern Canada, covering approximately 1.14 million square kilometers of diverse terrain dominated by boreal forest ecosystems. This region represents one of Canada’s most fire-prone areas, where the interplay of climate, vegetation, and geography creates conditions conducive to regular fire activity. The landscape is fundamentally shaped by the Canadian Shield, an ancient bedrock formation characterized by exposed Precambrian rock, countless lakes, and wetlands that create a complex mosaic of fire-susceptible and fire-resistant terrain.
At the heart of this region lies Yellowknife, the territorial capital, situated on the northern shore of Great Slave Lake. With approximately 20,000 residents, it serves as the primary urban center and economic hub of the NWT. The city’s boreal forest setting is dominated by black spruce and jack pine, species that have evolved with fire and even depend on it for regeneration. Black spruce, in particular, has adapted to release its seeds from sealed cones when exposed to the intense heat of forest fires, while jack pine stands often require fire to maintain their ecological dominance.
The region experiences a subarctic climate, characterized by long, cold winters where temperatures regularly drop below -30°C, and short, warm summers with extended daylight hours. This climate regime creates a distinct fire season that, historically, has been concentrated in the brief summer months. However, climate change is increasingly extending this fire season at both ends, creating new challenges for fire management and community safety.
Infrastructure in the NWT is notably limited, with Yellowknife connected to southern Canada by a single highway (Highway 3) and air travel serving as a crucial transportation link. This isolation is even more pronounced for the numerous remote communities scattered throughout the territory, many of which are accessible only by air or seasonal ice roads. This limited transportation network creates particular challenges for emergency management and evacuation scenarios during severe fire events, as demonstrated during recent fire seasons.
The 2023 Yellowknife Wildfire Crisis
The summer of 2023 brought unprecedented wildfire activity to the Northwest Territories, culminating in the evacuation of Yellowknife in mid-August. The immediate threat came from the Behchokǫ̀ wildfire (designated ZF015), which began southwest of Yellowknife and rapidly expanded under extreme fire weather conditions. By August 14th, the fire had grown to over 163,000 hectares and was advancing approximately 1.2 kilometers per hour toward the territorial capital.
On August 16th, territorial officials issued a mandatory evacuation order for all 20,000 Yellowknife residents, marking the largest evacuation in NWT history. The evacuation presented significant logistical challenges due to limited egress routes. Highway 3, the only road connecting Yellowknife to southern Canada, became a critical lifeline with evacuees forming long convoys south toward Alberta. Meanwhile, military and civilian aircraft conducted one of the largest air evacuations in Canadian history, prioritizing elderly residents, those with medical needs, and individuals without personal vehicles.
The fire’s behavior was particularly concerning due to the extreme conditions: temperatures reaching 30°C, relative humidity below 30%, and strong, erratic winds gusting over 40 km/h. These conditions, combined with the region’s unusually dry summer, created what fire behavior specialists termed “crown fire conditions” - where fires spread rapidly through the upper canopy of the forest. By August 17th, the fire had reached within 15 kilometers of Yellowknife’s city limits.
Emergency response efforts focused on establishing protective measures around the city, including fire breaks and sprinkler systems. Fire crews, supported by the Canadian military, worked to protect critical infrastructure including the airport, power plants, and communications facilities. The nearby communities of Dettah and Ndilǫ also faced evacuation orders, highlighting the fire’s broad regional impact.
The evacuation order remained in effect for 21 days, with residents unable to return until September 6th. This extended displacement period had significant social and economic impacts on the community, while also demonstrating the challenges of protecting isolated northern communities from wildfire threats. The event served as a stark reminder of the increasing vulnerability of northern communities to extreme fire events in a changing climate.
Wildfire and Climate Change
The unprecedented scale and intensity of the 2023 NWT fire season exemplifies the broader impacts of climate change on northern fire regimes. Canada’s northern regions are warming at approximately three times the global average rate, fundamentally altering the conditions that have historically governed wildfire behavior. This accelerated warming is creating longer fire seasons, with snow melting earlier in spring and winter arriving later in fall, extending the period during which fires can ignite and spread. The warming trend is particularly pronounced in the Northwest Territories, where average temperatures have increased by 2.3°C since 1948, compared to a 1.1°C increase globally.
These changing conditions are challenging traditional approaches to fire management and community protection. Historical fire behavior models, based on decades of past observations, are becoming less reliable as new climate patterns emerge. Fire intensity is increasing due to drier fuels and more frequent extreme weather events, while fire spread patterns are becoming more erratic and less predictable. For northern communities like Yellowknife, these changes represent an existential challenge: traditional fire breaks may prove insufficient, evacuation routes may become more frequently threatened, and the resources required for fire suppression may exceed historical requirements. The 2023 evacuation of Yellowknife serves as a harbinger of the complex challenges that northern communities face in adapting to this new reality of fire risk in a warming climate.
Geospatial Analysis in Emergency Response and Resource Management
Modern wildfire management increasingly relies on geospatial analysis and remote sensing technologies to support critical decision-making. Natural resource managers and emergency response teams leverage multiple satellite platforms, weather data, and spatial analysis techniques to monitor fire progression, assess risks, and coordinate response efforts in near-real-time. During the Yellowknife evacuation, these tools proved instrumental in tracking the Behchokǫ̀ wildfire’s progression, predicting its behavior, and making informed evacuation decisions.
The integration of multiple data sources - from active fire detections to smoke forecasts - allows emergency managers to develop a comprehensive understanding of evolving fire situations. NASA’s FIRMS system provides near-real-time fire locations, while GOES satellites track smoke plume development and movement. When combined with infrastructure data and analyzed through modern spatial techniques, these datasets enable rapid assessment of threats to communities, critical infrastructure, and evacuation routes. Post-fire, these same tools support damage assessment and recovery planning through the analysis of burn severity and infrastructure impacts.
The Yellowknife wildfire crisis of 2023 serves as an ideal case study for demonstrating these analytical approaches. In the following analysis, we will explore how different remote sensing products can be integrated to monitor and assess wildfire impacts, particularly in northern communities where traditional ground-based monitoring may be limited. We will demonstrate techniques for processing and analyzing multiple satellite data streams, assessing infrastructure vulnerability, and quantifying fire impacts - the same types of analyses that support real-world emergency response and resource management decisions. This practical application of geospatial analysis techniques provides valuable insights into both the technical methods and their critical role in modern fire management.
Analysis
Data Science Review
Before beginning, please note that this lesson uses the R programming language and the following R packages:
sf: Essential for handling vector spatial data and manipulation of geographic objects.stars: Extends spatial functionality to handle raster data and spatiotemporal arrays, particularly useful for satellite imagery.osmdata: Provides interface for OpenStreetMap data, enabling programmatic access to road networks and infrastructure data.aws.s3: Enables direct access to AWS S3 storage, essential for retrieving GOES satellite data.httr2: Modern HTTP client for R, used for making API requests and handling responses.ncdf4: Interface to Unidata netCDF format, crucial for working with satellite data files.concaveman: Implements concave hull algorithms for creating boundaries around point sets.dbscan: Provides density-based clustering algorithms for spatial point pattern analysis.dplyr: Part of tidyverse, essential for data manipulation and transformation.ggplot2: Creates elegant data visualizations and is the primary plotting tool in this lesson.tidyterra: Bridges terra spatial objects with tidyverse tools, particularly for visualization.terra: Advanced package for spatial data handling, particularly raster data processing.knitr: Generates dynamic reports and handles table formatting.kableExtra: Enhances table formatting capabilities with additional styling options.units: Handles unit conversions and dimensional analysis in spatial calculations.
If you’d like to learn more about the functions used in this lesson, you can refer to the documentation on their respective websites.
The sf package is fundamental for spatial data operations, while stars handles raster and spatiotemporal data. osmdata enables access to OpenStreetMap data for infrastructure analysis. aws.s3 and httr2 are crucial for data access via cloud storage and APIs respectively.
The ncdf4 package is essential for working with satellite data files, while concaveman and dbscan provide specialized spatial analysis capabilities. dplyr and ggplot2 form the backbone of data manipulation and visualization, enhanced by tidyterra for spatial data integration.
terra provides advanced spatial data processing capabilities, particularly for raster analysis. knitr and kableExtra handle report generation and table formatting, while units ensures proper handling of physical quantities and conversions.
Make sure these packages are installed before you begin working with the code in this document.
Initial Data Processing & Study Area Definition
To begin our analysis, we’ll need to define our study area around Yellowknife. The City of Yellowknife provides municipal boundary data through their open data portal. We’ll use this boundary as our primary study area and create buffered zones for different aspects of our analysis.
First, let’s set up our data directory and download the boundary file:
Next, we’ll unzip the downloaded file and read the municipal boundary layer. The data comes in a geodatabase format and is already in UTM Zone 11N (EPSG:26911):
# Unzip the file
utils::unzip("data/yellowknife_boundary.zip",
exdir = "data/yellowknife_boundaries")
# Read the municipal boundary layer
yk_boundary <- sf::st_read(
"data/yellowknife_boundaries/yk_Municipal_Boundary_gdb.gdb",
layer = "yk_municipal_boundary"
) |>
identity()Reading layer `yk_municipal_boundary' from data source
`C:\Users\jbrin\Documents\packages\TOPSTSCHOOL-disasters\data\yellowknife_boundaries\yk_Municipal_Boundary_gdb.gdb'
using driver `OpenFileGDB'
Simple feature collection with 1 feature and 1 field
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 627843.7 ymin: 6922154 xmax: 639196 ymax: 6937449
Projected CRS: NAD83 / UTM zone 11NLet’s examine the structure of our boundary object to ensure we have the data we need:
# Examine the structure of the boundary object
utils::str(yk_boundary)Classes 'sf' and 'data.frame': 1 obs. of 2 variables:
$ OBJECTID: int 2
$ SHAPE :sfc_MULTIPOLYGON of length 1; first list element: List of 1
..$ :List of 1
.. ..$ : num [1:9, 1:2] 627844 632135 631950 638805 639196 ...
..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
- attr(*, "sf_column")= chr "SHAPE"
- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
..- attr(*, "names")= chr "OBJECTID"Finally, we’ll create a basic map of our study area using ggplot2. This will help us verify that we’ve properly loaded the boundary data:
# Create a basic map of our study area
ggplot2::ggplot() +
ggplot2::geom_sf(data = yk_boundary,
fill = "lightblue",
color = "black") +
ggplot2::theme_minimal() +
ggplot2::labs(
title = "Yellowknife Municipal Boundary",
subtitle = paste0("Municipal Area: ",
round(sf::st_area(yk_boundary)/1000000, 2),
" sq km"),
caption = "Data source: City of Yellowknife Open Data Portal"
) +
ggplot2::coord_sf()
This boundary will serve as our primary reference for analyzing the wildfire impacts around Yellowknife. Throughout our analysis, we’ll create different buffer zones around this boundary to capture various aspects of the fire event, such as: - A 50km buffer for analyzing local fire detections - A 500km buffer for analyzing smoke plume dispersion - Custom buffers for specific analyses of evacuation routes and infrastructure
Data Science Review
When working with spatial data it’s important to understand a few key concepts:
Coordinate Reference Systems (CRS): Different spatial datasets might use different coordinate systems. The Yellowknife boundary data comes in UTM Zone 11N (EPSG:26911), which is appropriate for this region of Canada.
Spatial Objects in R: We’re using the
sf(Simple Features) package, which is the modern standard for handling vector spatial data in R. Thesfpackage integrates well with the tidyverse and provides efficient tools for spatial operations.Data Formats: Spatial data comes in many formats. Here we’re working with a geodatabase (.gdb), but you might also encounter shapefiles (.shp), GeoJSON, and other formats. The
sfpackage can handle most common spatial data formats.
Fire Detection
Satellite-based fire detection has transformed our ability to monitor and respond to wildfires, particularly in remote regions like Canada’s Northwest Territories. These systems use thermal anomaly detection algorithms that identify areas significantly warmer than their surroundings, enabling near real-time monitoring of fire activity across vast landscapes. When combined with ground observations and other data sources, satellite fire products provide crucial information for fire managers, emergency responders, and researchers studying fire behavior and impacts.
Several major satellite systems provide operational fire detection products. NASA’s MODIS (Moderate Resolution Imaging Spectroradiometer) instruments aboard the Terra and Aqua satellites have been providing global fire detections at 1km resolution since 2000, making them valuable for long-term studies. The VIIRS (Visible Infrared Imaging Radiometer Suite) sensors on the Suomi-NPP and NOAA-20 satellites offer improved detection capabilities with 375m resolution, able to detect smaller fires and provide more precise fire perimeter mapping. Geostationary satellites like GOES-16 and GOES-17 provide rapid fire detection updates every 5-15 minutes, though at coarser spatial resolution. The European Space Agency’s Sentinel-3 satellites also provide global fire detection through their Sea and Land Surface Temperature Radiometer (SLSTR) instrument.
In this analysis, we’ll focus on two complementary fire products that helped track the Yellowknife fires. The VIIRS 375m active fire product from NOAA-20, accessed through NASA’s FIRMS (Fire Information for Resource Management System), provides detailed fire detection points with associated fire radiative power measurements. We’ll combine this with the VIIRS Burned Area product (VNP64A1), which uses changes in surface reflectance to map the final extent of burned areas at 500m resolution. Together, these products allow us to analyze both the progression of active burning and the ultimate footprint of the fires.
Fire Information for Resource Management System (FIRMS)
NASA’s Fire Information for Resource Management System (FIRMS) distributes near real-time active fire data within 3 hours of satellite observation. The system synthesizes fire detections from both MODIS and VIIRS instruments, though we’ll focus on the VIIRS 375m data from NOAA-20 for its enhanced spatial resolution. Each FIRMS record provides not just the location of detected fires, but also their intensity measured as Fire Radiative Power (FRP) in megawatts. FRP serves as a proxy for fire intensity, with higher values generally indicating more intense burning and greater fuel consumption. This combination of location, timing, and intensity data makes FIRMS an invaluable tool for tracking fire progression and identifying areas of most intense fire activity. During the Yellowknife evacuation, FIRMS data provided critical information about fire location and intensity to emergency managers and the public.
First, we need to create our search area and acquire relevant basemap data. We’ll create a 50km buffer around Yellowknife to capture the local fire activity:
# Create search area and get bounding box
search_area <- sf::st_buffer(yk_boundary, 50000)
search_area_wgs84 <- sf::st_transform(search_area, 4326)
bbox <- sf::st_bbox(search_area_wgs84) |>
round(4)
# Download satellite imagery for our area
basemap_tiles <- maptiles::get_tiles(
search_area_wgs84,
provider = "Esri.WorldImagery",
zoom = 11 # Adjust zoom level as needed
)Next, we’ll query the FIRMS API to get fire detections for our area. The FIRMS API requires registration for an API key, which should be stored securely in your environment variables:
# Get FIRMS data
firms_url <- sprintf(
"https://firms.modaps.eosdis.nasa.gov/api/area/csv/%s/VIIRS_NOAA20_NRT/%s,%s,%s,%s/10/2023-08-10",
Sys.getenv("FIRMS_API_KEY"),
bbox["xmin"], bbox["ymin"],
bbox["xmax"], bbox["ymax"]
)
resp <- httr2::request(firms_url) |>
httr2::req_perform()We’ll process the FIRMS data into a spatial object and add proper datetime fields:
# Process FIRMS data
firms_points <- httr2::resp_body_string(resp) |>
textConnection() |>
utils::read.csv() |>
sf::st_as_sf(
coords = c("longitude", "latitude"),
crs = 4326
) |>
sf::st_transform(sf::st_crs(yk_boundary))
# Add proper datetime field and date
firms_points$datetime <- as.POSIXct(
paste(firms_points$acq_date,
sprintf("%04d", firms_points$acq_time)),
format = "%Y-%m-%d %H%M",
tz = "UTC"
)
firms_points$date <- as.Date(firms_points$acq_date)To understand the temporal progression of fire activity, we’ll calculate daily statistics:
# Calculate daily statistics
daily_stats <- firms_points |>
as.data.frame() |>
dplyr::group_by(date) |>
dplyr::summarise(
detections = dplyr::n(),
mean_frp = mean(frp),
max_frp = max(frp),
total_frp = sum(frp)
)
Understanding Fire Radiative Power (FRP)
Fire Radiative Power (FRP) is a measure of the rate of radiant heat output from a fire, measured in megawatts (MW). Higher FRP values generally indicate more intense fires. The VIIRS sensor can detect:
- Smaller/cooler fires: ~1-10 MW
- Moderate fires: ~10-100 MW
- Large/intense fires: >100 MW
FRP helps fire managers assess fire intensity and potential severity, though it’s important to note that factors like cloud cover can affect detection capabilities.
Now we can create a series of visualizations to understand the spatial and temporal patterns of fire activity. First, let’s create a map showing all fire detections colored by their intensity (FRP):
# Create main static visualization
p1 <- ggplot2::ggplot() +
tidyterra::geom_spatraster_rgb(data = basemap_tiles) +
ggplot2::geom_sf(data = search_area,
fill = NA,
color = "gray60",
linetype = "dashed") +
ggplot2::geom_sf(data = yk_boundary,
fill = NA,
color = "#00BFFF",
linewidth = 1.2) +
ggplot2::geom_sf(data = firms_points,
ggplot2::aes(color = frp),
alpha = 0.7,
size = 2) +
ggplot2::scale_color_gradientn(
colors = c("orange", "orange red", "red", "dark red", "purple"),
name = "Fire Radiative\nPower (MW)"
) +
ggplot2::theme_minimal() +
ggplot2::labs(
title = "Yellowknife Fire Detections",
subtitle = paste("August 10-23, 2023\nTotal detections:", nrow(firms_points)),
caption = "Data: VIIRS 375m from NOAA-20 satellite"
)
print(p1)
To understand how fire activity changed over time, we’ll create two temporal plots - one showing the number of daily detections and another showing fire intensity:
# Define evacuation date
evac_date <- as.Date("2023-08-16")
# Create temporal progression plot
p2 <- ggplot2::ggplot(daily_stats, ggplot2::aes(x = date)) +
ggplot2::geom_vline(xintercept = evac_date,
color = "purple",
linetype = "dashed",
alpha = 0.85,
linewidth = 1.5) +
ggplot2::geom_line(ggplot2::aes(y = detections), color = "#00BFFF", linewidth = 1) +
ggplot2::geom_point(ggplot2::aes(y = detections), color = "#00BFFF", size = 3) +
ggplot2::theme_minimal() +
ggplot2::labs(
title = "Daily Fire Detections",
subtitle = "Purple line indicates evacuation date (Aug 16)",
y = "Number of detections",
x = "Date"
)
print(p2)
Finally, let’s examine how fire intensity changed throughout the period:
# Create FRP intensity plot with proper legend
p3 <- ggplot2::ggplot(daily_stats, ggplot2::aes(x = date)) +
ggplot2::geom_vline(xintercept = evac_date,
color = "purple",
linetype = "dashed",
alpha = 0.85,
linewidth = 1.5) +
ggplot2::geom_line(ggplot2::aes(y = mean_frp, color = "Mean FRP"), linewidth = 1) +
ggplot2::geom_point(ggplot2::aes(y = mean_frp, color = "Mean FRP"), size = 3) +
ggplot2::geom_line(ggplot2::aes(y = max_frp, color = "Max FRP"), linewidth = 1) +
ggplot2::geom_point(ggplot2::aes(y = max_frp, color = "Max FRP"), size = 3) +
ggplot2::theme_minimal() +
ggplot2::labs(
title = "Fire Radiative Power Over Time",
subtitle = "Purple line indicates evacuation date (Aug 16)",
y = "Fire Radiative Power (MW)",
x = "Date"
) +
ggplot2::scale_color_manual(
values = c("Mean FRP" = "#00BFFF", "Max FRP" = "#0040ff"),
name = "Metric"
)
print(p3)
Let’s summarize our findings in a clear table:
# Create summary statistics table using kableExtra
stats_df <- data.frame(
Metric = c(
"Date Range",
"Total Detections",
"Mean FRP",
"Max FRP",
"Peak Detection Date",
"Peak Daily Detections"
),
Value = c(
paste(min(firms_points$date), "to", max(firms_points$date)),
as.character(nrow(firms_points)),
paste(round(mean(firms_points$frp), 2), "MW"),
paste(round(max(firms_points$frp), 2), "MW"),
as.character(daily_stats$date[which.max(daily_stats$detections)]),
as.character(daily_stats$detections[which.max(daily_stats$detections)])
)
)
# Create formatted HTML table
knitr::kable(stats_df,
format = "html",
caption = "Yellowknife Fire Analysis Summary") |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE,
position = "left"
) |>
kableExtra::row_spec(0, bold = TRUE)| Metric | Value |
|---|---|
| Date Range | 2023-08-10 to 2023-08-19 |
| Total Detections | 8896 |
| Mean FRP | 10.26 MW |
| Max FRP | 1643.92 MW |
| Peak Detection Date | 2023-08-15 |
| Peak Daily Detections | 2347 |
Data Science Review
Several key data science concepts are demonstrated in this analysis:
API Interaction: We use the FIRMS API to programmatically access fire detection data, demonstrating how to construct and execute API requests.
Spatial Data Processing: The analysis shows how to:
- Transform between coordinate systems
- Create spatial buffers
- Join spatial and temporal data
Data Visualization Best Practices:
- Using appropriate color scales for fire intensity
- Including reference lines for key events (evacuation date)
- Providing clear labels and legends
- Creating multi-panel visualizations for different aspects of the data
Fire Complex Perimeters
While individual fire detections give us point-based information about fire activity, understanding the overall fire extent and progression requires analyzing these points as continuous complexes. We’ll use spatial clustering and concave hull techniques to create meaningful fire perimeters from our point detections.
First, let’s create a function that will generate fire perimeters for a given date using concave hull geometry:
# Create perimeters with adjusted parameters
# The concavity parameter (default = 1) controls how tightly the hull wraps around points:
# - Lower values (e.g., 0.5) create more detailed, tighter boundaries
# - Higher values (e.g., 2) create smoother, more generalized boundaries
# The threshold parameter sets the minimum segment length (in coordinate units)
# to prevent very small segments in the hull
create_fire_perimeters <- function(points,
date,
# Controls how tightly hull fits points
concavity = 1,
# Minimum segment length
threshold = 1) {
daily_points <- points |>
dplyr::filter(date == !!date)
if (nrow(daily_points) < 3) {
return(NULL)
}
coords <- sf::st_coordinates(daily_points)
hull <- concaveman::concaveman(coords, concavity = concavity, length_threshold = threshold)
polygon_sf <- sf::st_polygon(list(hull)) |>
sf::st_sfc(crs = sf::st_crs(points)) |>
sf::st_sf()
perimeter <- polygon_sf |>
dplyr::mutate(
date = date,
n_points = nrow(daily_points),
mean_frp = mean(daily_points$frp),
area_km2 = as.numeric(units::set_units(sf::st_area(polygon_sf), "km^2"))
)
return(perimeter)
}
Understanding Concave Hulls
A concave hull (also known as a concave closure or α-shape) provides a more realistic boundary around a set of points compared to a convex hull. While a convex hull creates the smallest convex polygon that contains all points, a concave hull can “wrap” more tightly around the points, capturing the true shape of the point distribution. This is particularly useful for fire perimeter mapping, where we want to represent the actual extent of fire activity rather than just the outer bounds.
Next, let’s create daily perimeters for our analysis period:
# Create daily perimeters
unique_dates <- unique(firms_points$date)
fire_perimeters <-
do.call(rbind,
lapply(unique_dates, function(d)
create_fire_perimeters(firms_points, d, 1, 0.01)))
# Find dates of maximum area and point counts
max_stats <- list(
area_date = fire_perimeters |>
dplyr::arrange(dplyr::desc(area_km2)) |>
dplyr::slice(1) |>
dplyr::pull(date),
points_date = fire_perimeters |>
dplyr::arrange(dplyr::desc(n_points)) |>
dplyr::slice(1) |>
dplyr::pull(date)
)However, single perimeters can sometimes oversimplify complex fire situations. To better represent distinct fire complexes, we’ll create a clustering function that can identify separate fire groups:
# eps (epsilon) is the maximum distance between two points for them to be considered neighbors
# eps = 5000 means points within 5km of each other are considered part of the same cluster
# This value was chosen based on:
# - Typical fire spread distances in boreal forests
# - Spatial resolution of VIIRS fire detections (375m)
# - Desire to identify distinct fire complexes while avoiding over-segmentation
create_clustered_perimeters <- function(points,
date,
eps = 5000) {
# Filter points for the given date
daily_points <- points |>
dplyr::filter(date == !!date)
if (nrow(daily_points) < 3) {
return(NULL)
}
# Get coordinates for clustering
coords <- sf::st_coordinates(daily_points)
# Perform DBSCAN clustering
clusters <- dbscan::dbscan(coords, eps = eps, minPts = 3)
# Add cluster information to points
daily_points$cluster <- clusters$cluster
# Create separate hull for each cluster
hulls <- daily_points |>
dplyr::filter(cluster > 0) |> # Remove noise points (cluster = 0)
dplyr::group_by(cluster) |>
dplyr::group_map(function(cluster_points, cluster_id) {
if (nrow(cluster_points) < 3) {
return(NULL)
}
coords <- sf::st_coordinates(cluster_points)
hull <- concaveman::concaveman(coords,
concavity = 1,
length_threshold = 1)
current_poly <- sf::st_polygon(list(hull)) |>
sf::st_sfc(crs = sf::st_crs(points)) |>
sf::st_sf()
current_area <- as.numeric(
units::set_units(sf::st_area(current_poly), "km^2")
)
current_poly |>
dplyr::mutate(
date = date,
cluster = cluster_id$cluster,
n_points = nrow(cluster_points),
mean_frp = mean(cluster_points$frp),
area_km2 = current_area
)
})
# Combine all hulls for this date
hulls_combined <- do.call(rbind, hulls)
return(hulls_combined)
}
Understanding DBSCAN Clustering
DBSCAN (Density-Based Spatial Clustering of Applications with Noise) is particularly well-suited for identifying fire complexes because it:
- Can find clusters of arbitrary shape
- Automatically determines the number of clusters
- Can identify outliers as noise
- Uses two intuitive parameters:
- eps: The maximum distance between two points to be considered neighbors
- minPts: The minimum number of points to form a cluster
Now we can create clustered perimeters for all dates and visualize the progression:
# Create clustered perimeters for all dates
all_clustered_perimeters <- do.call(rbind,
lapply(unique_dates, function(d) {
perims <- create_clustered_perimeters(firms_points, d)
if (!is.null(perims)) {
perims$date <- d
}
return(perims)
})
)
# Create progression plot with layered hulls
p4_progression <- ggplot2::ggplot() +
tidyterra::geom_spatraster_rgb(data = basemap_tiles) +
ggplot2::geom_sf(
data = search_area,
fill = NA,
color = "#00BFFF",
linetype = "dashed"
) +
ggplot2::geom_sf(
data = yk_boundary,
fill = NA,
color = "#00BFFF",
linewidth = 1.2
) +
# Add hulls in reverse chronological order
ggplot2::geom_sf(
data = all_clustered_perimeters |>
dplyr::arrange(dplyr::desc(date)),
ggplot2::aes(fill = date),
alpha = 0.4,
linewidth = 0.8
) +
ggplot2::scale_fill_viridis_c(
name = "Date",
option = "magma",
labels = function(x) format(as.Date(x), "%B %d"),
breaks = unique_dates[seq(1, length(unique_dates), by = 2)] # Show every other date
) +
ggplot2::theme_minimal() +
ggplot2::labs(
title = "Fire Complex Progression",
subtitle = paste(format(min(unique_dates), "%B %d"),
"to",
format(max(unique_dates), "%B %d")),
caption = "Note: Clustered using DBSCAN algorithm (eps = 3km)\nEarlier dates displayed on top"
)
print(p4_progression)
This is a bit messy, so let’s also examine the fire complexes on days of peak activity by comparing the maximum area date with the maximum detection count date:
# Create dataset for faceted plot
comparison_perimeters <- fire_perimeters |>
dplyr::filter(date %in% c(max_stats$area_date, max_stats$points_date)) |>
dplyr::mutate(
type = dplyr::case_when(
date == max_stats$area_date ~ "Maximum Hull Area",
date == max_stats$points_date ~ "Maximum Detection Count"
)
)
# Get corresponding VIIRS points for these dates
comparison_points <- firms_points |>
dplyr::filter(date %in% c(max_stats$area_date, max_stats$points_date)) |>
dplyr::mutate(
type = dplyr::case_when(
date == max_stats$area_date ~ "Maximum Hull Area",
date == max_stats$points_date ~ "Maximum Detection Count"
)
)
# Get clustered perimeters for comparison dates
comparison_clusters <- rbind(
create_clustered_perimeters(firms_points, max_stats$area_date) |>
dplyr::mutate(type = "Maximum Hull Area"),
create_clustered_perimeters(firms_points, max_stats$points_date) |>
dplyr::mutate(type = "Maximum Detection Count")
)Update the labels to include dates for better context:
# Modify the facet labels to include dates
comparison_clusters <- comparison_clusters |>
dplyr::mutate(
type = dplyr::case_when(
type == "Maximum Hull Area" ~
paste("Maximum Hull Area -", format(max_stats$area_date, "%B %d")),
type == "Maximum Detection Count" ~
paste(
"Maximum Detection Count -",
format(max_stats$points_date, "%B %d")
)
)
)
# Update points labels to match
comparison_points <- comparison_points |>
dplyr::mutate(
type = dplyr::case_when(
type == "Maximum Hull Area" ~
paste("Maximum Hull Area -", format(max_stats$area_date, "%B %d")),
type == "Maximum Detection Count" ~
paste(
"Maximum Detection Count -",
format(max_stats$points_date, "%B %d")
)
)
)Finally, create a comparative visualization:
# Create faceted comparison plot with clustered hulls
p4c <- ggplot2::ggplot() +
# Base layers
tidyterra::geom_spatraster_rgb(data = basemap_tiles) +
ggplot2::geom_sf(
data = search_area,
fill = NA,
color = "#00BFFF",
linetype = "dashed"
) +
ggplot2::geom_sf(
data = yk_boundary,
fill = NA,
color = "#00BFFF",
linewidth = 1.2
) +
# Add clustered hull polygons
ggplot2::geom_sf(
data = comparison_clusters,
ggplot2::aes(fill = factor(cluster)),
alpha = 0.60,
linewidth = 1.2,
show.legend = FALSE
) +
# Add VIIRS points
ggplot2::geom_sf(
data = comparison_points,
color = "darkred",
size = 1,
alpha = 0.25
) +
# Facet by type
ggplot2::facet_wrap( ~ type) +
# Style
ggplot2::scale_color_viridis_c(name = "Fire Radiative\nPower (MW)",
option = "inferno") +
ggplot2::theme_minimal() +
ggplot2::labs(
title = "Comparison of Fire Complexes",
caption = "Note: Points clustered using DBSCAN algorithm (eps = 5km)"
)
print(p4c)
# Print statistics about clusters
cat("\nCluster Statistics:\n")
Cluster Statistics:for (date_type in unique(comparison_clusters$type)) {
clusters_date <- comparison_clusters |>
dplyr::filter(type == date_type)
cat("\n", date_type, ":\n")
cat("Number of distinct fire complexes:",
length(unique(clusters_date$cluster)),
"\n")
cat("Total area across all complexes:",
round(sum(clusters_date$area_km2), 2),
"km²\n")
cat("Largest single complex:",
round(max(clusters_date$area_km2), 2),
"km²\n")
}
Maximum Hull Area - August 10 :
Number of distinct fire complexes: 8
Total area across all complexes: 363.22 km²
Largest single complex: 227.14 km²
Maximum Detection Count - August 15 :
Number of distinct fire complexes: 9
Total area across all complexes: 355.42 km²
Largest single complex: 215.12 km²Interpreting Fire Complex Analysis
When analyzing fire complexes, several key metrics help us understand the fire situation:
Number of Distinct Complexes: Multiple fire complexes can indicate either natural fire spread patterns or multiple ignition sources.
Complex Size: The size of individual complexes helps fire managers assess resource needs and potential threats.
Complex Distribution: The spatial arrangement of complexes relative to infrastructure and communities informs evacuation planning.
Temporal Progression: Understanding how complexes grow and merge over time helps predict future fire behavior.
Knowledge Check: FIRMS
What is the primary advantage of FIRMS data for wildfire monitoring?
- It provides weather forecasts for fire-prone areas
- It delivers near real-time fire detections within 3 hours of satellite observation
- It measures ground-level temperature
- It predicts future fire spread
VIIRS Burned Area
The VIIRS Burned Area product (VNP64A1) provides a comprehensive assessment of burned areas by detecting changes in surface reflectance characteristics that occur after fire. Unlike FIRMS data which captures active burning, VNP64A1 identifies the full extent of burned areas at 500m resolution, including areas that may have burned without being detected as active fires due to cloud cover or timing of satellite overpasses. The product assigns each pixel a “burn date” indicating when the burn was first detected, allowing for temporal analysis of fire progression. This approach is particularly valuable for post-fire assessment and validation of active fire detections, though it has limitations in near real-time monitoring due to its processing latency. For the Yellowknife fires, VNP64A1 data helps us understand the total area impacted and verify the progression patterns suggested by FIRMS detections.
The VNP64A1 data used in this analysis was acquired through NASA’s Earthdata Cloud Search Portal (earthdata.nasa.gov/search), which allows users to search, preview, and download specific granules based on location and date range. After identifying granules covering the Yellowknife region during August 2023, the data was downloaded in HDF5 format for processing into a folder called “data”.
First, let’s locate and read our VIIRS burned area files. We’ll focus on files from August 1 (day 213) to September 1 (day 244):
# Get a list of VIIRS Burned Area files
viirs_files <- list.files(
"data/VNP64A1.002/VNP64A1_002-20241031_164531/",
pattern = "A2023(213|244).*\\.hdf$",
full.names = TRUE
)
# Print the files we'll be working with
print("VIIRS Burned Area files to process:")[1] "VIIRS Burned Area files to process:"print(viirs_files) [1] "data/VNP64A1.002/VNP64A1_002-20241031_164531/VNP64A1.A2023213.h11v02.002.2023276101150.hdf"
[2] "data/VNP64A1.002/VNP64A1_002-20241031_164531/VNP64A1.A2023213.h11v03.002.2023276102309.hdf"
[3] "data/VNP64A1.002/VNP64A1_002-20241031_164531/VNP64A1.A2023213.h12v02.002.2023276103300.hdf"
[4] "data/VNP64A1.002/VNP64A1_002-20241031_164531/VNP64A1.A2023213.h12v03.002.2023276101306.hdf"
[5] "data/VNP64A1.002/VNP64A1_002-20241031_164531/VNP64A1.A2023213.h13v02.002.2023276103154.hdf"
[6] "data/VNP64A1.002/VNP64A1_002-20241031_164531/VNP64A1.A2023244.h11v02.002.2023312150112.hdf"
[7] "data/VNP64A1.002/VNP64A1_002-20241031_164531/VNP64A1.A2023244.h11v03.002.2023312150411.hdf"
[8] "data/VNP64A1.002/VNP64A1_002-20241031_164531/VNP64A1.A2023244.h12v02.002.2023312150240.hdf"
[9] "data/VNP64A1.002/VNP64A1_002-20241031_164531/VNP64A1.A2023244.h12v03.002.2023312150107.hdf"
[10] "data/VNP64A1.002/VNP64A1_002-20241031_164531/VNP64A1.A2023244.h13v02.002.2023312145922.hdf"
Understanding VIIRS Burned Area Data
The VNP64A1 product contains several data layers, but the primary “Burn Date” layer records the day of year when a burn was detected. Special values include: - 0: Unburned - -1: Unmapped - -2: Water
Now we’ll process each HDF file, extracting the burn date information:
# Create empty list to store processed burn date layers
burn_dates_list <- list()
# Process each HDF file
for (i in seq_along(viirs_files)) {
# Read the HDF file using terra
burn_rast <- terra::rast(viirs_files[i])
# Print layer names for first file to verify structure
if (i == 1) {
print("\nLayer names in the HDF file:")
print(names(burn_rast))
}
# Extract burn date layer (first layer)
burn_date <- burn_rast[[1]] # Select "Burn Date" layer
# Process burn dates
# Convert to more meaningful values and handle special codes
values <- terra::values(burn_date)
values[values <= 0] <- NA
terra::values(burn_date) <- values
# Store processed layer
burn_dates_list[[i]] <- burn_date
print(paste("Processed file", i, "of", length(viirs_files)))
}[1] "\nLayer names in the HDF file:"
[1] "\"Burn Date\"" "\"Burn Date Uncertainty\""
[3] "QA" "\"First Day\""
[5] "\"Last Day\""
[1] "Processed file 1 of 10"
[1] "Processed file 2 of 10"
[1] "Processed file 3 of 10"
[1] "Processed file 4 of 10"
[1] "Processed file 5 of 10"
[1] "Processed file 6 of 10"
[1] "Processed file 7 of 10"
[1] "Processed file 8 of 10"
[1] "Processed file 9 of 10"
[1] "Processed file 10 of 10"Next, we’ll combine our processed tiles into a single mosaic and align it with our study area:
# Mosaic all tiles together in one step
burn_mosaic <- terra::mosaic(burn_dates_list[[1]],
burn_dates_list[[2]],
burn_dates_list[[3]],
burn_dates_list[[4]],
fun = "last") # Use most recent valid observation
# Reproject to match our other data (EPSG:4326)
burn_mosaic_4326 <- terra::project(burn_mosaic, "EPSG:4326")
# Add a buffer to capture the full extent of the fire
study_area_buffer <- sf::st_buffer(yk_boundary, 50000)
# Convert study area buffer to SpatVector for terra operations
study_area_vect <- terra::vect(study_area_buffer)
# Project study area to match burn mosaic's original CRS (sinusoidal)
study_area_proj <- terra::project(study_area_vect, terra::crs(burn_mosaic))
# Crop and mask the burn data
burn_clipped <- terra::crop(burn_mosaic, study_area_proj)
burn_clipped <- terra::mask(burn_clipped, study_area_proj)
# Now project everything to EPSG:4326 for visualization
burn_clipped_4326 <- terra::project(burn_clipped, "EPSG:4326")
study_area_4326 <- terra::project(study_area_vect, "EPSG:4326")
yk_boundary_4326 <- terra::project(terra::vect(yk_boundary), "EPSG:4326")To prepare for visualization, we’ll convert our data to a format suitable for ggplot and create meaningful time periods:
# Convert raster to data frame for plotting
burn_df <- terra::as.data.frame(burn_clipped_4326, xy = TRUE)
names(burn_df)[3] <- "burn_date"
# Convert burn dates to actual dates
burn_df$date <- as.Date("2023-01-01") + burn_df$burn_date - 1
# Create categories focusing on August
burn_df$period <- cut(burn_df$date,
breaks = as.Date(c("2023-08-01", "2023-08-08",
"2023-08-15", "2023-08-22",
"2023-09-01")),
labels = c("Aug 1-7", "Aug 8-14",
"Aug 15-21", "Aug 22-31"),
include.lowest = TRUE
)Now we can create our first visualization showing the progression of burned areas:
# Create the visualization with discrete categories
p1 <- ggplot2::ggplot() +
ggplot2::geom_tile(data = burn_df,
ggplot2::aes(x = x, y = y, fill = period)) +
ggplot2::geom_sf(data = study_area_4326,
fill = NA,
color = "gray60",
linetype = "dashed") +
ggplot2::geom_sf(data = yk_boundary_4326,
fill = NA,
color = "#00BFFF",
linewidth = 1) +
# Using a sequential red color scheme to show progression
ggplot2::scale_fill_manual(
values = c("Aug 1-7" = "#FEE5D9",
"Aug 8-14" = "#FCAE91",
"Aug 15-21" = "#FB6A4A", # Evacuation period
"Aug 22-31" = "#CB181D"),
na.value = "transparent",
name = "Burn Period",
drop = FALSE
) +
ggplot2::theme_minimal() +
ggplot2::labs(
title = "Burn Progression Around Yellowknife",
subtitle = paste("August 2023\n",
"(Evacuation ordered August 16th)"),
caption = "Data: VIIRS VNP64A1 Burned Area Product"
)
print(p1)
# Print summary of area burned in each period
print("\nArea burned by period (km²):")[1] "\nArea burned by period (km²):"burn_summary <- table(burn_df$period) * 0.25 # Convert pixel count to km²
print(burn_summary)
Aug 1-7 Aug 8-14 Aug 15-21 Aug 22-31
21.25 119.75 30.25 3.50 For validation and comparison, let’s create a visualization that overlays our FIRMS active fire detections on the burn areas:
# Create comparison with FIRMS data
firms_points_4326 <- terra::project(terra::vect(firms_points), "EPSG:4326")
p2 <- ggplot2::ggplot() +
ggplot2::geom_tile(data = burn_df,
ggplot2::aes(x = x, y = y, fill = period)) +
ggplot2::geom_sf(data = study_area_4326,
fill = NA,
color = "gray60",
linetype = "dashed") +
ggplot2::geom_sf(data = yk_boundary_4326,
fill = NA,
color = "#00BFFF",
linewidth = 1) +
ggplot2::geom_sf(data = firms_points,
ggplot2::aes(color = frp),
alpha = 0.6,
size = 0.8) +
ggplot2::scale_fill_manual(
values = c("Aug 1-7" = "#FEE5D9",
"Aug 8-14" = "#FCAE91",
"Aug 15-21" = "#FB6A4A",
"Aug 22-31" = "#CB181D"),
na.value = "transparent",
name = "Burn Period",
drop = FALSE
) +
ggplot2::scale_color_viridis_c(
name = "Fire Radiative\nPower (MW)"
) +
ggplot2::theme_minimal() +
ggplot2::labs(
title = "Burn Scars and Active Fire Detections",
subtitle = "VIIRS Burned Area with FIRMS Hotspots",
caption = "Data: VNP64A1 and FIRMS VIIRS"
)
print(p2)
Interpreting Burned Area Analysis
When analyzing burned areas, several key considerations help validate and interpret the results:
Temporal Progression: The timing of detected burns should generally align with active fire detections
Spatial Agreement: Burn scars should correspond with areas where we detected active fires
Resolution Effects: The 500m resolution of the burned area product means small fires might be missed
Edge Effects: Areas near the edge of satellite swaths or scene boundaries may show artificial patterns
The combination of active fire detections (FIRMS) and burned area products (VNP64A1) provides a more complete picture of fire impacts than either product alone.
Knowledge Check: VIIRS Burned Area
Why is the VIIRS Burned Area product (VNP64A1) complementary to FIRMS data?
- It provides higher temporal resolution
- It measures ground temperature
- It detects changes in surface reflectance to identify burned areas, even those that may have been missed by active fire detection
- It predicts where fires will occur
GOES-18 Smoke Plume Analysis
Tracking smoke plume extent and movement is crucial for public health and emergency management during wildfires. While several satellite products can detect smoke, aerosol optical depth (AOD) measurements from geostationary satellites provide particularly valuable data for monitoring smoke plumes due to their high temporal frequency. AOD quantifies the degree to which aerosols (including smoke particles) prevent the transmission of light through the atmosphere, serving as a proxy for smoke density and air quality impacts.
The GOES-18 ABI (Advanced Baseline Imager) Level 2 Aerosol Optical Depth (AOD) product provides near real-time monitoring of aerosol loading across North America at 15-minute intervals. This high temporal resolution allows for detailed tracking of smoke plume development and transport, though at a coarser spatial resolution than some polar-orbiting satellites. The product provides AOD values ranging from 0 (clear sky) to 5 (very thick aerosol), with values above 1 typically indicating significant smoke impacts. During the Yellowknife evacuation, this data helped emergency managers assess potential impacts on communities and transportation routes as smoke complicated both ground and air evacuation efforts.
I’ll help break down the “Visualizing a Snapshot” code section into meaningful chunks with clear narrative:
Visualizing a Snapshot
Let’s begin by selecting a timestamp during peak fire activity to examine smoke conditions. We’ll focus on August 15th, the day before the evacuation order:
# Pick a timestamp during peak fire activity
test_date <- as.Date("2023-08-15") # Day before evacuation
test_hour <- 18 # Mid-afternoon local timeThe GOES-18 data is stored in Amazon Web Services (AWS) S3 buckets, organized by product, year, day of year, and hour. We’ll construct the appropriate path and retrieve the file list:
# Get GOES file list
goes_url <- sprintf(
"%s/%s/%s/%02d",
"ABI-L2-AODF",
format(test_date, "%Y"),
format(test_date, "%j"),
test_hour
)
files_df <- aws.s3::get_bucket_df(
bucket = "noaa-goes18",
prefix = goes_url,
region = "us-east-1",
max = 20
)
Understanding GOES Data Organization
GOES satellite data follows a consistent naming convention: - ABI-L2-AODF: Product name (Aerosol Optical Depth Full Disk) - Year/DayOfYear/Hour: Temporal organization - Multiple files per hour represent different scan times
Now we’ll download and read the most recent file for our selected timestamp:
# Get most recent file
latest_file <- files_df |>
dplyr::arrange(dplyr::desc(LastModified)) |>
dplyr::slice(1)
# Download and read file
url <- sprintf("https://noaa-goes18.s3.amazonaws.com/%s", latest_file$Key)
temp_file <- tempfile(fileext = ".nc")
download.file(url, temp_file, mode = "wb", quiet = TRUE)
aod_data <- stars::read_stars(temp_file, sub = "AOD", quiet = TRUE)
unlink(temp_file)To focus our analysis on the Yellowknife area, we’ll create two buffer zones around the city:
# Create visualization extent (100km buffer for context)
yk_buffer_big <- sf::st_buffer(yk_boundary, 100000)
yk_buffer_small <- sf::st_buffer(yk_boundary, 50000) # 50km for analysis
# Transform to GOES projection
yk_buffer_goes <- sf::st_transform(yk_buffer_big, sf::st_crs(aod_data))
yk_small_goes <- sf::st_transform(yk_buffer_small, sf::st_crs(aod_data))
# Crop to our area of interest
aod_local <- sf::st_crop(aod_data, yk_buffer_goes)Finally, we can create a visualization of the smoke conditions:
# Create visualization
ggplot2::ggplot() +
stars::geom_stars(data = aod_local) +
ggplot2::geom_sf(data = sf::st_transform(yk_boundary, sf::st_crs(aod_data)),
fill = NA, color = "#00BFFF", linewidth = 1) +
ggplot2::geom_sf(data = yk_small_goes,
fill = NA, color = "blue", linetype = "dashed") +
ggplot2::scale_fill_viridis_c(
name = "AOD",
na.value = NA,
option = "inferno",
limits = c(0, 3)
) +
ggplot2::theme_minimal() +
ggplot2::labs(
title = "GOES-18 AOD Data Around Yellowknife",
subtitle = format(test_date, "%B %d, %Y %H:00 UTC"),
caption = "Red: City boundary\nBlue dashed: 50km analysis buffer"
)
Interpreting AOD Values
Aerosol Optical Depth values can be roughly interpreted as: - 0.0-0.1: Clear sky - 0.1-0.3: Light haze - 0.3-0.5: Moderate haze - 0.5-1.0: Heavy haze - >1.0: Very heavy smoke/haze
Higher values indicate more particles in the atmospheric column, suggesting thicker smoke plumes and potentially worse air quality.
Temporal Analysis of AOD Values
To understand how smoke conditions evolved during the evacuation period, we’ll analyze AOD values over several days. First, let’s create some helper functions to streamline our data processing:
# Functions for GOES data handling
goes_date_format <- function(date) {
list(
year = format(date, "%Y"),
doy = format(date, "%j")
)
}
# Function to construct https URL from Key
construct_url <- function(key, bucket = "noaa-goes18") {
sprintf("https://%s.s3.amazonaws.com/%s", bucket, key)
}
# Function to get GOES files
get_goes_files <- function(date, hour, product = "ABI-L2-AODF", satellite = "goes18") {
# Get date components
date_parts <- goes_date_format(date)
# Construct prefix
prefix <- sprintf(
"%s/%s/%s/%02d/",
product,
date_parts$year,
date_parts$doy,
hour
)
message("Searching prefix: ", prefix, " in ", satellite)
# Get file list from bucket
files_df <- aws.s3::get_bucket_df(
bucket = paste0("noaa-", satellite),
prefix = prefix,
region = "us-east-1",
max = 20
)
if (nrow(files_df) > 0) {
files_df$LastModified <- as.POSIXct(
files_df$LastModified,
format = "%Y-%m-%dT%H:%M:%S.000Z",
tz = "UTC"
)
files_df$satellite <- satellite
}
return(files_df)
}Now we’ll create a function to process AOD data for specific timestamps:
# Function to get and process AOD data for a specific timestamp
get_aod_data <- function(date, hour) {
# Get file list
files_df <- get_goes_files(
date = date,
hour = hour,
product = "ABI-L2-AODF",
satellite = "goes18"
)
if (nrow(files_df) == 0) {
message("No files found for: ", date, " ", hour, ":00 UTC")
return(NULL)
}
# Get most recent file
latest_file <- files_df |>
dplyr::arrange(dplyr::desc(LastModified)) |>
dplyr::slice(1)
# Download and read file
url <- construct_url(latest_file$Key, bucket = "noaa-goes18")
temp_file <- tempfile(fileext = ".nc")
utils::download.file(
url = url,
destfile = temp_file,
mode = "wb",
quiet = TRUE
)
# Read AOD data
aod_data <- stars::read_stars(
temp_file,
sub = "AOD",
quiet = TRUE
)
# Clean up
unlink(temp_file)
return(aod_data)
}
GOES Data Processing Tips
The GOES AWS bucket provides the most recent 7 days of data at no cost. When downloading multiple files: - Use appropriate time intervals to avoid redundant data - Consider file size limitations and download times - Remember that UTC timestamps are used for data organization
We’ll also need a function to calculate statistics for our area of interest:
# Modify calculate_local_stats function to strip units
calculate_local_stats <- function(aod_data, local_area) {
# Transform and crop to local area
tryCatch({
local_goes <- sf::st_transform(local_area, sf::st_crs(aod_data))
aod_local <- sf::st_crop(aod_data, local_goes)
# Check if we have any valid data
if (sum(!is.na(aod_local[[1]])) > 0) {
# Convert to numeric to strip units
# GOES AOD values come with units attached (dimensionless optical depth units)
# We strip these for easier calculation and compatibility with other functions
# Raw AOD values typically range from 0-5, where:
# 0-0.1: Clear sky
# 0.1-0.3: Light haze
# >1.0: Heavy smoke
values <- as.numeric(aod_local[[1]])
# Convert to numeric to strip units
values <- as.numeric(aod_local[[1]])
stats <- list(
mean_aod = mean(values, na.rm = TRUE),
max_aod = max(values, na.rm = TRUE),
valid_pixels = sum(!is.na(values)),
total_pixels = length(values)
)
} else {
stats <- list(
mean_aod = NA_real_,
max_aod = NA_real_,
valid_pixels = 0L,
total_pixels = length(aod_local[[1]])
)
}
return(stats)
}, error = function(e) {
message("Error in processing: ", e$message)
return(list(
mean_aod = NA_real_,
max_aod = NA_real_,
valid_pixels = 0L,
total_pixels = 0L
))
})
}Now we can define our analysis period and process multiple timestamps:
# Define analysis period (5 days before to 3 days after evacuation)
analysis_dates <- seq(
as.Date("2023-08-11"), # 5 days before evacuation
as.Date("2023-08-19"), # 3 days after evacuation
by = "day"
)
# Define daily hours to analyze (daylight hours in UTC)
analysis_hours <- c(16, 18, 20, 22) # Approximately 10am-4pm local time
# Create analysis timestamp combinations
analysis_times <- expand.grid(
date = analysis_dates,
hour = analysis_hours
)
# Process data for all analysis times
aod_results <- lapply(1:nrow(analysis_times), function(i) {
date <- analysis_times$date[i]
hour <- analysis_times$hour[i]
message("\nProcessing: ", date, " ", hour, ":00 UTC")
# Get AOD data
aod_data <- get_aod_data(date, hour)
if (is.null(aod_data)) {
return(list(
datetime = as.POSIXct(
paste(date, sprintf("%02d:00:00", hour)),
tz = "UTC"
),
mean_aod = NA_real_,
max_aod = NA_real_,
valid_pixels = 0L,
total_pixels = 0L
))
}
# Calculate statistics
stats <- calculate_local_stats(aod_data, yk_buffer_small)
# Add timestamp
stats$datetime <- as.POSIXct(
paste(date, sprintf("%02d:00:00", hour)),
tz = "UTC"
)
return(stats)
})
# Convert to data frame and handle NaN/Inf values
aod_df <- dplyr::bind_rows(aod_results) |>
dplyr::mutate(
mean_aod = replace(mean_aod, !is.finite(mean_aod), NA),
max_aod = replace(max_aod, !is.finite(max_aod), NA)
)For context, let’s provide a reference table for interpreting AOD values in terms of air quality provided by the World Health Organization:
# Create WHO threshold reference table
who_thresholds <- data.frame(
Category = c("Good", "Moderate", "Unhealthy for Sensitive Groups",
"Unhealthy", "Very Unhealthy", "Hazardous"),
`PM2.5 (µg/m³)` = c("0-12", "12.1-35.4", "35.5-55.4",
"55.5-150.4", "150.5-250.4", ">250.5"),
`Approximate AOD` = c("<0.1", "0.1-0.3", "0.3-0.5",
"0.5-1.0", "1.0-2.0", ">2.0"),
check.names = FALSE
)Now we can create visualizations to show how smoke conditions evolved. First, a time series of AOD values:
# Create time series visualization with WHO thresholds
p1 <- ggplot2::ggplot() +
# Add WHO threshold reference lines
ggplot2::geom_hline(yintercept = 0.5, linetype = "dashed", color = "yellow", alpha = 0.5) +
ggplot2::geom_hline(yintercept = 1.0, linetype = "dashed", color = "orange", alpha = 0.5) +
ggplot2::geom_hline(yintercept = 2.0, linetype = "dashed", color = "red", alpha = 0.5) +
# Add evacuation reference line
ggplot2::geom_vline(
xintercept = as.POSIXct("2023-08-16 00:00:00", tz = "UTC"),
color = "red",
linetype = "dashed",
alpha = 0.5
) +
# Add AOD lines
ggplot2::geom_line(data = aod_df,
ggplot2::aes(x = datetime, y = mean_aod, color = "Mean AOD")) +
ggplot2::geom_point(data = aod_df,
ggplot2::aes(x = datetime, y = mean_aod, color = "Mean AOD")) +
ggplot2::geom_line(data = aod_df,
ggplot2::aes(x = datetime, y = max_aod, color = "Max AOD")) +
ggplot2::geom_point(data = aod_df,
ggplot2::aes(x = datetime, y = max_aod, color = "Max AOD")) +
# Customize appearance
ggplot2::scale_color_manual(
name = "Metric",
values = c("Mean AOD" = "blue", "Max AOD" = "red")
) +
ggplot2::theme_minimal() +
ggplot2::labs(
title = "Aerosol Optical Depth Around Yellowknife",
subtitle = "50km radius, GOES-18 observations",
x = "Date/Time (UTC)",
y = "Aerosol Optical Depth",
caption = "Red dashed vertical line: Evacuation date\nHorizontal lines: WHO thresholds (yellow=unhealthy, orange=very unhealthy, red=hazardous)"
)
print(p1)
Finally, let’s summarize our findings:
# Create summary statistics table
valid_coverage <- aod_df |>
dplyr::filter(total_pixels > 0) |>
dplyr::summarise(
coverage = mean(valid_pixels/total_pixels * 100, na.rm = TRUE)
) |>
dplyr::pull(coverage)
stats_table <- data.frame(
Metric = c("Analysis Period", "Mean AOD", "Maximum AOD", "Average Valid Coverage"),
Value = c(
paste(format(min(aod_df$datetime), "%B %d %H:%M UTC"), "to",
format(max(aod_df$datetime), "%B %d %H:%M UTC")),
sprintf("%.3f", mean(aod_df$mean_aod, na.rm = TRUE)),
sprintf("%.3f", max(aod_df$max_aod, na.rm = TRUE)),
sprintf("%.1f%%", valid_coverage)
)
)
# Display tables
knitr::kable(stats_table,
caption = "AOD Analysis Summary Statistics")| Metric | Value |
|---|---|
| Analysis Period | August 11 16:00 UTC to August 19 22:00 UTC |
| Mean AOD | 0.751 |
| Maximum AOD | 4.997 |
| Average Valid Coverage | 17.0% |
knitr::kable(who_thresholds,
caption = "WHO Air Quality Categories and Approximate AOD Values")| Category | PM2.5 (µg/m³) | Approximate AOD |
|---|---|---|
| Good | 0-12 | <0.1 |
| Moderate | 12.1-35.4 | 0.1-0.3 |
| Unhealthy for Sensitive Groups | 35.5-55.4 | 0.3-0.5 |
| Unhealthy | 55.5-150.4 | 0.5-1.0 |
| Very Unhealthy | 150.5-250.4 | 1.0-2.0 |
| Hazardous | >250.5 | >2.0 |
Understanding Temporal AOD Patterns
When interpreting temporal AOD patterns, consider: 1. Diurnal cycles - smoke often settles overnight and lifts during the day 2. Wind patterns - can cause rapid changes in local smoke conditions 3. Cloud interference - can cause gaps in the data 4. Vertical distribution - AOD measures the entire column, not just surface conditions
Visualizing the Maximum Plume Spread
To understand the full spatial extent of smoke impacts, we’ll analyze a larger area around Yellowknife during peak smoke conditions. First, let’s create a function to process plume data for specific timestamps:
# Function to process plume data for a specific timestamp
analyze_plume <- function(date, hour, plume_buffer) {
# Get GOES file
files_df <- get_goes_files(
date = date,
hour = hour,
product = "ABI-L2-AODF",
satellite = "goes18"
)
if (nrow(files_df) == 0) {
message("No files found for: ", date, " ", hour, ":00 UTC")
return(NULL)
}
# Get most recent file for timestamp
latest_file <- files_df |>
dplyr::arrange(desc(LastModified)) |>
dplyr::slice(1)
# Download and read data
url <- construct_url(latest_file$Key, bucket = "noaa-goes18")
temp_file <- tempfile(fileext = ".nc")
utils::download.file(url, temp_file, mode = "wb", quiet = TRUE)
aod_data <- stars::read_stars(temp_file, sub = "AOD", quiet = TRUE)
unlink(temp_file)
# Transform and crop to plume analysis extent
plume_goes <- sf::st_transform(plume_buffer, sf::st_crs(aod_data))
plume_aod <- sf::st_crop(aod_data, plume_goes)
# Calculate plume statistics
values <- as.numeric(plume_aod[[1]])
stats <- list(
datetime = as.POSIXct(
paste(date, sprintf("%02d:00:00", hour)),
tz = "UTC"
),
mean_aod = mean(values, na.rm = TRUE),
max_aod = max(values, na.rm = TRUE),
heavy_smoke_area = sum(values > 2, na.rm = TRUE) * 0.0081, # Approx km² per pixel
moderate_smoke_area = sum(values > 1, na.rm = TRUE) * 0.0081,
data = plume_aod
)
return(stats)
}
Quantifying Smoke Extent
We classify smoke coverage into two categories: - Moderate smoke: AOD > 1 (significant visibility reduction) - Heavy smoke: AOD > 2 (severe air quality impacts)
The area calculation (0.0081 km² per pixel) is based on the GOES-R ABI resolution at nadir.
Now let’s analyze a focused period around the evacuation:
# Define analysis period (3 days centered on evacuation)
analysis_dates <- seq(
as.Date("2023-08-15"), # Day before evacuation
as.Date("2023-08-17"), # Day after evacuation
by = "day"
)
# Define daily hours to analyze (focusing on peak hours)
analysis_hours <- c(16, 18, 20, 22) # ~10am-4pm local time
# Create timestamp combinations
analysis_times <- expand.grid(
date = analysis_dates,
hour = analysis_hours
)
# Define larger analysis extent (500km radius for plume tracking)
plume_buffer <- sf::st_buffer(yk_boundary, 500000) # 500km
# Process all timestamps
plume_results <- lapply(1:nrow(analysis_times), function(i) {
date <- analysis_times$date[i]
hour <- analysis_times$hour[i]
message("\nProcessing: ", date, " ", hour, ":00 UTC")
analyze_plume(date, hour, plume_buffer)
})Let’s extract statistics and create a time series of plume extent:
# Extract statistics into data frame
plume_stats <- dplyr::bind_rows(lapply(plume_results, function(x) {
if (!is.null(x)) {
data.frame(
datetime = x$datetime,
mean_aod = x$mean_aod,
max_aod = x$max_aod,
heavy_smoke_area = x$heavy_smoke_area,
moderate_smoke_area = x$moderate_smoke_area
)
}
}))
# Create time series visualization of plume metrics
p1_plume <- ggplot2::ggplot(plume_stats) +
ggplot2::geom_vline(
xintercept = as.POSIXct("2023-08-16 00:00:00", tz = "UTC"),
color = "red",
linetype = "dashed",
alpha = 0.5
) +
ggplot2::geom_line(ggplot2::aes(x = datetime, y = moderate_smoke_area,
color = "Moderate Smoke (AOD > 1)")) +
ggplot2::geom_point(ggplot2::aes(x = datetime, y = moderate_smoke_area,
color = "Moderate Smoke (AOD > 1)")) +
ggplot2::geom_line(ggplot2::aes(x = datetime, y = heavy_smoke_area,
color = "Heavy Smoke (AOD > 2)")) +
ggplot2::geom_point(ggplot2::aes(x = datetime, y = heavy_smoke_area,
color = "Heavy Smoke (AOD > 2)")) +
ggplot2::scale_color_manual(
name = "Smoke Category",
values = c("Moderate Smoke (AOD > 1)" = "orange",
"Heavy Smoke (AOD > 2)" = "red")
) +
ggplot2::theme_minimal() +
ggplot2::labs(
title = "Smoke Plume Coverage Around Yellowknife",
subtitle = "500km radius analysis extent",
x = "Date/Time (UTC)",
y = "Affected Area (km²)",
caption = "Red dashed line indicates evacuation date"
)
print(p1_plume)
Now let’s examine conditions during the peak smoke period:
# Identify peak smoke period
peak_time <- plume_stats |>
dplyr::arrange(desc(moderate_smoke_area)) |>
dplyr::slice(1)
peak_data <- plume_results[[which(
sapply(plume_results, function(x) {
!is.null(x) && x$datetime == peak_time$datetime
})
)]]
# Create peak plume visualization
p2_plume <- ggplot2::ggplot() +
stars::geom_stars(data = peak_data$data) +
ggplot2::geom_sf(data = sf::st_transform(yk_boundary,
sf::st_crs(peak_data$data)),
fill = NA, color = "#00BFFF", linewidth = 1) +
ggplot2::geom_sf(data = sf::st_transform(plume_buffer,
sf::st_crs(peak_data$data)),
fill = NA, color = "blue", linetype = "dashed") +
ggplot2::scale_fill_viridis_c(
name = "AOD",
na.value = NA,
option = "inferno",
limits = c(0, 3)
) +
ggplot2::theme_minimal() +
ggplot2::labs(
title = "Peak Smoke Plume Extent",
subtitle = format(peak_data$datetime, "%B %d, %Y %H:%M UTC"),
caption = paste(
"Red: City boundary\nBlue dashed: 500km analysis extent\n",
"Affected area (AOD > 1):",
round(peak_time$moderate_smoke_area), "km²"
)
)
print(p2_plume)
# Create summary statistics table
stats_summary <- data.frame(
Metric = c(
"Analysis Period",
"Peak AOD Value",
"Maximum Affected Area (AOD > 1)",
"Maximum Heavy Smoke Area (AOD > 2)",
"Average Affected Area (AOD > 1)",
"Peak Observation Time"
),
Value = c(
paste(format(min(plume_stats$datetime), "%B %d %H:%M UTC"), "to",
format(max(plume_stats$datetime), "%B %d %H:%M UTC")),
sprintf("%.2f", max(plume_stats$max_aod)),
sprintf("%.0f km²", max(plume_stats$moderate_smoke_area)),
sprintf("%.0f km²", max(plume_stats$heavy_smoke_area)),
sprintf("%.0f km²", mean(plume_stats$moderate_smoke_area)),
format(peak_time$datetime, "%B %d %H:%M UTC")
)
)
# Print summary table
knitr::kable(stats_summary,
caption = "Smoke Plume Analysis Summary",
format = "pipe")| Metric | Value |
|---|---|
| Analysis Period | August 15 16:00 UTC to August 17 22:00 UTC |
| Peak AOD Value | 5.00 |
| Maximum Affected Area (AOD > 1) | 28 km² |
| Maximum Heavy Smoke Area (AOD > 2) | 5 km² |
| Average Affected Area (AOD > 1) | 11 km² |
| Peak Observation Time | August 16 16:00 UTC |
Interpreting Plume Extent
When analyzing smoke plume extent: 1. Consider both intensity (AOD values) and spatial coverage 2. Remember that AOD measures total column aerosol - smoke may be at various heights 3. Cloud cover can create gaps in the data 4. The 500km analysis radius helps capture regional-scale impacts
GOES-18 Wind Patterns
Understanding atmospheric transport is crucial for predicting smoke impacts during wildfire events. Wind patterns not only influence fire behavior but also determine how smoke will affect surrounding communities. During the Yellowknife evacuation, wind conditions played a dual role: affecting both fire progression and creating hazardous air quality conditions that complicated evacuation efforts.
Traditionally, smoke transport has been monitored through a combination of ground-based weather stations and atmospheric models. However, these methods can have significant limitations in remote regions like the Northwest Territories, where weather stations are sparse and model validation data is limited. Satellite-derived wind products help fill these gaps by providing broad spatial coverage and frequent updates. The GOES-18 Derived Motion Winds (DMW) product is particularly valuable as it provides wind measurements at multiple atmospheric levels by tracking the movement of features (such as clouds and water vapor) in successive satellite images.
For our analysis, we’ll focus on winds during the period of peak smoke intensity, combining GOES-18 DMW data with our aerosol optical depth measurements. This integration helps us understand: 1. The atmospheric conditions that influenced smoke transport 2. Potential evacuation routes that may have been impacted 3. Areas where convergent winds may have concentrated smoke 4. The broader regional transport patterns that affected air quality across the Northwest Territories
First, let’s create a function to retrieve and process wind data.
# Function to get wind data for a specific timestamp
get_wind_data <- function(date, hour) {
# Get DMW (Derived Motion Winds) file
files_df <- get_goes_files(
date = date,
hour = hour,
product = "ABI-L2-DMWVF", # Full disk product
satellite = "goes18"
)
if (nrow(files_df) == 0) {
message("No wind files found for: ", date, " ", hour, ":00 UTC")
return(NULL)
}
# Get most recent file
latest_file <- files_df |>
dplyr::arrange(desc(LastModified)) |>
dplyr::slice(1)
# Download and read data
url <- construct_url(latest_file$Key, bucket = "noaa-goes18")
temp_file <- tempfile(fileext = ".nc")
utils::download.file(url, temp_file, mode = "wb", quiet = TRUE)
# Read with ncdf4
nc <- ncdf4::nc_open(temp_file)
# Read wind data
tryCatch({
wind_data <- list(
speed = ncdf4::ncvar_get(nc, "wind_speed"),
direction = ncdf4::ncvar_get(nc, "wind_direction"),
lat = ncdf4::ncvar_get(nc, "lat"),
lon = ncdf4::ncvar_get(nc, "lon"),
quality = ncdf4::ncvar_get(nc, "DQF") # Quality flag
)
# Create data frame of valid winds (where quality == 0)
valid_idx <- which(wind_data$quality == 0)
wind_df <- data.frame(
longitude = wind_data$lon[valid_idx],
latitude = wind_data$lat[valid_idx],
speed = wind_data$speed[valid_idx],
direction = wind_data$direction[valid_idx]
)
# Remove any NA or invalid values
wind_df <- wind_df[complete.cases(wind_df) &
wind_df$speed >= 0 &
wind_df$speed < 150 & # reasonable wind speed limits
wind_df$direction >= 0 &
wind_df$direction <= 360, ]
# Create spatial object
wind_sf <- sf::st_as_sf(wind_df,
coords = c("longitude", "latitude"),
crs = 4326)
# Convert to stars
wind_stars <- stars::st_as_stars(wind_sf)
ncdf4::nc_close(nc)
unlink(temp_file)
return(wind_stars)
}, error = function(e) {
message("Error reading wind data: ", e$message)
ncdf4::nc_close(nc)
unlink(temp_file)
return(NULL)
})
}
Wind Analysis
The GOES-18 Derived Motion Winds (DMW) product: - Tracks features in successive satellite images to derive wind vectors - Provides wind speed and direction at multiple atmospheric levels - Includes quality flags to identify reliable measurements - Helps understand atmospheric transport patterns
To visualize wind vectors, we need to convert speed and direction into vector components:
# Function to convert wind speed and direction to u,v components
calculate_wind_components <- function(speed, direction) {
# Convert direction from meteorological to mathematical angle
math_angle <- (270 - direction) * pi / 180
# Calculate u and v components
u <- speed * cos(math_angle)
v <- speed * sin(math_angle)
return(list(u = u, v = v))
}
# Get wind data for peak plume time
peak_winds <- get_wind_data(
as.Date(peak_time$datetime),
lubridate::hour(peak_time$datetime)
)Now we’ll process the wind data to match our analysis area:
# Transform wind data to match AOD projection
wind_transformed <- sf::st_transform(peak_winds, sf::st_crs(peak_data$data))
# Extract wind components
wind_speed <- as.numeric(wind_transformed[[1]])
wind_direction <- as.numeric(wind_transformed[[2]])
wind_components <- calculate_wind_components(wind_speed, wind_direction)
# Convert to spatial features and clip to study area
wind_sf <- sf::st_as_sf(wind_transformed)
# Properly handle intersections using st_intersects
plume_goes <- sf::st_transform(plume_buffer, sf::st_crs(peak_data$data))
intersects <- sf::st_intersects(wind_sf, plume_goes)
intersecting_indices <- which(lengths(intersects) > 0)
wind_sf_cropped <- wind_sf[intersecting_indices, ]To create a clear visualization, we’ll subsample the wind vectors:
# Get coordinates of cropped winds
wind_coords <- sf::st_coordinates(wind_sf_cropped)
# Create regular subsample indices (every nth point)
n_points <- nrow(wind_coords)
subsample_interval <- max(1, floor(n_points / 300)) # Aim for ~300 vectors
subsample <- seq(1, n_points, by = subsample_interval)
# Create wind dataframe with subsampled points
wind_df <- data.frame(
x = wind_coords[subsample, 1],
y = wind_coords[subsample, 2],
u = wind_components$u[intersecting_indices][subsample],
v = wind_components$v[intersecting_indices][subsample],
speed = wind_speed[intersecting_indices][subsample]
)Finally, let’s create our integrated visualization of smoke conditions and wind patterns:
# Create the plot with subsampled wind vectors
p3_plume <- ggplot2::ggplot() +
ggplot2::geom_sf(data = sf::st_transform(yk_boundary,
sf::st_crs(peak_data$data)),
fill = NA, color = "#00BFFF", linewidth = 0.75) +
ggplot2::geom_sf(data = sf::st_transform(plume_buffer,
sf::st_crs(peak_data$data)),
fill = NA, color = "blue", linetype = "dashed", linewidth = 2) +
# Add wind vectors
# The scaling factor of 2500 was chosen to create visible vectors at this map scale:
# - Raw u,v components are in m/s
# - Map coordinates are in meters (EPSG:3857)
# - 2500 means a 1 m/s wind creates a 2.5km vector on the map
# - This scaling provides good visibility while avoiding overcrowding
# - Adjust this value based on your map extent and wind speeds
ggplot2::geom_segment(data = wind_df,
ggplot2::aes(x = x, y = y,
xend = x + u*2500, # Scale factor for vector length
yend = y + v*2500,
color = speed),
arrow = grid::arrow(length = ggplot2::unit(0.2, "cm")),
alpha = 1,
linewidth = 1.25) +
ggplot2::scale_color_viridis_c(
name = "Wind Speed\n(m/s)",
option = "turbo"
) +
ggplot2::theme_minimal() +
ggplot2::labs(
title = "Yellowknife GOES Wind Vectors",
subtitle = format(peak_data$datetime, "%B %d, %Y %H:%M UTC"),
caption = paste(
"Grey: City boundary\nBlue dashed: 500km analysis extent\n",
"Affected area (AOD > 1):",
round(peak_time$moderate_smoke_area), "km²"
)
) +
# Add coord_sf to ensure proper extent
ggplot2::coord_sf(xlim = sf::st_bbox(plume_goes)[c("xmin", "xmax")],
ylim = sf::st_bbox(plume_goes)[c("ymin", "ymax")])
print(p3_plume)
# Create wind statistics summary
wind_stats <- data.frame(
Metric = c(
"Mean Wind Speed",
"Max Wind Speed",
"Predominant Wind Direction",
"Wind Speed at Yellowknife"
),
Value = c(
sprintf("%.1f m/s", mean(wind_speed, na.rm = TRUE)),
sprintf("%.1f m/s", max(wind_speed, na.rm = TRUE)),
sprintf("%.0f°", median(wind_direction, na.rm = TRUE)),
sprintf("%.1f m/s", wind_speed[which.min(sf::st_distance(
sf::st_as_sf(wind_transformed, as_points = TRUE),
sf::st_transform(yk_boundary, sf::st_crs(wind_transformed))
))])
)
)
# Print statistics
knitr::kable(wind_stats,
caption = "Wind Conditions During Peak Plume Extent",
format = "pipe")| Metric | Value |
|---|---|
| Mean Wind Speed | 13.7 m/s |
| Max Wind Speed | 84.2 m/s |
| Predominant Wind Direction | 226° |
| Wind Speed at Yellowknife | 23.0 m/s |
Interpreting Wind Patterns
When analyzing wind vectors: 1. Vector length indicates wind speed 2. Vector direction shows where the wind is blowing towards 3. Color scaling highlights areas of stronger winds 4. Pattern variations can indicate: - Convergence zones where smoke may accumulate - Transport corridors for smoke movement - Areas of atmospheric mixing
Knowledge Check: GOES Wind Analysis
How does the GOES-18 Derived Motion Winds (DMW) product determine wind patterns?
- Using ground-based weather stations
- Through computer modeling
- By tracking the movement of features like clouds in successive satellite images
- By measuring atmospheric pressure
Evacuation Route Analysis
During the Yellowknife evacuation, the majority of residents traveled south along Highway 3 (Yellowknife Highway) to Highway 1 (Mackenzie Highway). These routes were critically important as they represented the only road-based evacuation options for most residents. Analyzing these evacuation corridors in relation to fire and smoke conditions helps us understand the challenges faced during the evacuation.
OpenStreetMap (OSM) provides a comprehensive source of road network data that is particularly valuable for remote regions like the Northwest Territories. As a collaborative mapping project, OSM relies on contributions from local knowledge, government data imports, and satellite imagery digitization. This crowdsourced approach often results in more up-to-date and detailed road networks than traditional commercial maps, especially in areas where official mapping resources may be limited. The OSM database includes not just road geometries, but also attributes like road classification, names, and reference numbers that help identify major evacuation routes.
For northern communities like Yellowknife, OSM’s road network data is crucial for emergency planning because it captures the limited transportation options available. The database accurately reflects the sparse road network characteristic of northern Canada, where communities are often connected by single highways that serve as critical lifelines during emergencies. Through the osmdata R package, we can programmatically access this road network data and focus our analysis on the key evacuation routes.
First, let’s obtain the road network data using OpenStreetMap (OSM) and identify our key evacuation routes:
# Transform plume buffer to WGS84 and get bbox for osm
plume_buffer_wgs84 <- sf::st_transform(plume_buffer, 4326)
bbox_wgs84 <- as.vector(sf::st_bbox(plume_buffer_wgs84))
# Get main evacuation route (Highway 3)
evac_route <- osmdata::opq(bbox = bbox_wgs84) |>
osmdata::add_osm_feature(
key = "highway",
value = c("primary", "trunk", "motorway")
) |>
osmdata::osmdata_sf()
OpenStreetMap Road Classifications
- motorway: highest grade highways
- trunk: major arterial routes
- primary: major regional routes In the Northwest Territories, most intercity routes are classified as trunk or primary roads.
Now we’ll filter for our specific evacuation routes:
# Filter for both Highway 3 and Highway 1 segments
evacuation_route <- evac_route$osm_lines |>
sf::st_transform(sf::st_crs(yk_boundary)) |>
dplyr::filter(
# Highway 3 segments
grepl("3", ref, fixed = TRUE) |
grepl("Highway 3", name, fixed = TRUE) |
grepl("Yellowknife Highway", name, fixed = TRUE) |
# Highway 1 segments
grepl("1", ref, fixed = TRUE) |
grepl("Highway 1", name, fixed = TRUE) |
(name == "Mackenzie Highway" & highway == "trunk")
) |>
sf::st_transform(sf::st_crs(peak_data$data))
# Transform evacuation routes if needed
evacuation_route <- evacuation_route |>
sf::st_transform(sf::st_crs(peak_data$data))
# Intersect routes to get those in the plume buffer
local_routes <- sf::st_intersection(evacuation_route, plume_goes)Let’s create a visualization of the evacuation routes in relation to our study area:
# Create local routes visualization
p4_local <- ggplot2::ggplot() +
# Base map with routes
ggplot2::geom_sf(data = plume_buffer,
fill = NA, color = "gray70", linetype = "dashed") +
ggplot2::geom_sf(data = local_routes,
color = "red",
linewidth = 1.5) +
ggplot2::geom_sf(data = sf::st_transform(yk_boundary,
sf::st_crs(peak_data$data)),
fill = NA, color = "#00BFFF", linewidth = 0.75) +
ggplot2::theme_minimal() +
ggplot2::labs(
title = "Evacuation Routes from Yellowknife",
subtitle = format(peak_data$datetime, "%B %d, %Y %H:%M UTC"),
caption = "Red: Primary evacuation routes\nLight blue: City boundary\nGray dashed: 50km buffer"
)
print(p4_local)
Evacuation Route Considerations
When analyzing evacuation routes, several factors are important: 1. Limited Options: In northern communities, there are often few alternate routes 2. Distance: Evacuees may need to travel long distances to reach safety 3. Vulnerability: Routes can be impacted by both fire and smoke 4. Capacity: Limited routes may face congestion during mass evacuation
Knowledge Check: Evacuation Routes
Why is OpenStreetMap particularly valuable for analyzing evacuation routes in remote regions like Northwest Territories?
- It provides real-time traffic updates
- It includes detailed elevation data
- It offers more up-to-date and detailed road networks through crowdsourced local knowledge
- It shows historical road conditions
Integrated Evacuation Conditions
By combining our analyses of smoke plumes, wind patterns, and evacuation routes, we can better understand the challenging conditions faced by evacuees and emergency managers. During wildfire evacuations, the interaction between smoke transport and evacuation routes is particularly critical - heavy smoke can reduce visibility, impact air quality, and potentially force evacuees to seek alternate routes.
Let’s start by getting the most relevant FIRMS data for our peak smoke period:
# Get FIRMS data for peak time with larger extent
bbox <- sf::st_bbox(sf::st_transform(plume_buffer, 4326)) |>
round(4)
firms_url <- sprintf(
"https://firms.modaps.eosdis.nasa.gov/api/area/csv/%s/VIIRS_NOAA20_NRT/%s,%s,%s,%s/10/%s",
Sys.getenv("FIRMS_API_KEY"),
bbox["xmin"], bbox["ymin"],
bbox["xmax"], bbox["ymax"],
format(peak_time$datetime, "%Y-%m-%d")
)
# Get the data
resp <- httr2::request(firms_url) |>
httr2::req_perform()
# Process into spatial data
firms_points_peak <- httr2::resp_body_string(resp) |>
textConnection() |>
utils::read.csv() |>
sf::st_as_sf(
coords = c("longitude", "latitude"),
crs = 4326
) |>
sf::st_transform(sf::st_crs(peak_data$data))Now we’ll create fire perimeters to show the extent of active burning:
# Create fire perimeters using our earlier clustering function
fire_perims <- create_clustered_perimeters(firms_points_peak,
date = as.Date(peak_time$datetime))
Understanding the Integrated Visualization
The following visualization combines several key elements: - Smoke plume extent (AOD values) - Active fire perimeters - Wind vectors showing transport patterns - Primary evacuation routes
This integration helps visualize the multiple hazards evacuees faced.
Now we can create our comprehensive visualization showing the interaction of all these factors:
# Create integrated visualization with larger extent
p4_integrated <- ggplot2::ggplot() +
# Smoke plume base
stars::geom_stars(data = peak_data$data, alpha = 0.50) +
# Evacuation routes
ggplot2::geom_sf(data = evacuation_route,
color = "black",
linewidth = 1.1) +
# Boundaries
ggplot2::geom_sf(data = sf::st_transform(yk_boundary,
sf::st_crs(peak_data$data)),
fill = NA, color = "#00BFFF", linewidth = 2) +
ggplot2::geom_sf(data = sf::st_transform(plume_buffer,
sf::st_crs(peak_data$data)),
fill = NA, color = "grey50", linetype = "dashed", linewidth = 1.5) +
# Fire perimeters
ggplot2::geom_sf(data = fire_perims,
fill = "orange",
alpha = 0.75,
color = "red",
linewidth = 0.5) +
# Wind vectors
ggplot2::geom_segment(data = wind_df,
ggplot2::aes(x = x, y = y,
xend = x + u*2500,
yend = y + v*2500,
color = speed),
arrow = grid::arrow(length = ggplot2::unit(0.2, "cm")),
alpha = 0.65,
linewidth = 1.1) +
# Scales
ggplot2::scale_fill_viridis_c(
name = "AOD (Plume Hazard)",
na.value = NA,
option = "inferno",
limits = c(0, 3)
) +
ggplot2::scale_color_gradientn(
name = "Wind Speed (m/s)",
colors = c("blue", "purple")
)+
ggplot2::theme_minimal() +
ggplot2::labs(
title = "Yellowknife, CA Wildfire Crisis: Evacuation Conditions",
subtitle = format(peak_data$datetime, "%B %d, %Y %H:%M UTC"),
caption = paste(
"Black: Primary evacuation routes (OpenStreetMap) | Light blue: City boundary\n",
"Orange areas: Active fire complexes (NASA FIRMS VIIRS 375m)\n",
"Wind vectors: GOES-18 Derived Motion Winds | Background: GOES-18 Aerosol Optical Depth\n",
"Gray dashed: 500km analysis extent"
)
)+
ggplot2::coord_sf(xlim = sf::st_bbox(plume_goes)[c("xmin", "xmax")],
ylim = sf::st_bbox(plume_goes)[c("ymin", "ymax")])
print(p4_integrated)
Analyzing Evacuation Conditions
When interpreting this integrated view, consider:
- Route Vulnerability:
- Proximity of fires to evacuation routes
- Smoke intensity along evacuation corridors
- Wind patterns that could drive fire or smoke toward routes
- Temporal Factors:
- This represents conditions at peak smoke intensity
- Conditions likely varied throughout the evacuation period
- Wind patterns influenced both fire spread and smoke transport
- Spatial Relationships:
- Limited evacuation options (single main route)
- Regional extent of smoke impacts
- Potential for route closures due to fire proximity
- Evacuation Management Implications:
- Need for real-time monitoring of multiple hazards
- Importance of alternative evacuation plans
- Value of integrated spatial analysis for decision support
This comprehensive visualization demonstrates the complex challenges faced during the Yellowknife evacuation. The combination of active fires, extensive smoke plumes, and limited evacuation routes created a particularly challenging scenario for both emergency managers and evacuees. The wind patterns during this period were especially critical, as they influenced both the movement of smoke and the potential for fire spread toward evacuation routes.
Knowledge Check: Comprehensive Analysis
Which combination of data sources provides the most complete understanding of evacuation conditions during a wildfire event?
- FIRMS data and weather forecasts only
- Wind patterns and smoke plumes only
- Integration of fire detections, smoke plumes, wind patterns, and evacuation routes
- Burned area analysis only
Conclusion
This lesson demonstrated how multiple remote sensing products can be integrated to analyze wildfire impacts, taking the 2023 Yellowknife evacuation as a case study. Through hands-on analysis of NASA’s FIRMS active fire detections, GOES-18 aerosol optical depth and wind data, and VIIRS burned area products, we gained insights into both the technical aspects of geospatial analysis and the real-world applications for emergency management.
This case study illustrates the vital role of Earth observation data in:
- Supporting emergency management decisions
- Documenting and analyzing extreme events
- Understanding the progression of natural disasters
- Validating ground-based observations with satellite data
The methods demonstrated here can be applied to other wildfire events or adapted for different types of natural hazards where multiple spatial datasets need to be integrated for comprehensive analysis.
Future Directions
While this analysis focused on basic integration of available datasets, future work could:
- Include additional data sources such as weather conditions and fuel moisture
- Develop automated processing chains for near-real-time analysis
- Incorporate ground validation data
- Analyze longer time series to understand historical context
- Integrate socioeconomic data to assess community impacts