In this blog post I’ll show how you can use R to map and analyze spatial data on wildfires in the US and in Colorado, using packages like {sf} and {leaflet}. As I write this, there are several major wildfires burning nearby along the Colorado Front Range (luckily they all seem to be under control now).
I downloaded the shapefiles for the Best available perimeters for all known wildland fires in the United States from the HIFLD Open website on August 6, 2024. HIFLD stands for Homeland Infrastructure Foundation-Level Data, and the HIFLD Open website contains a lot of interesting spatial data sets for the US.
This post is about how to work with spatial data in R, and not an up-to-date or official source of verified information about the fires. See also the following caveats from the data set:
After downloading the data and unzipping the file, the data can be read into R with the st_read() function from the {sf} package (Pebesma and Bivand 2023). The output gives some information about the spatial data including the coordinate reference system and bounding box.
Reading layer `Perimeters' from data source
`/Users/andy/Projects/Andy_Pickering_Portfolio/posts/wildfires_mapping/data/WFIGS_Interagency_Perimeters_-8845918407708086874/Perimeters.shp'
using driver `ESRI Shapefile'
Simple feature collection with 25141 features and 111 fields (with 9 geometries empty)
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -176.6504 ymin: 0.0003362746 xmax: 144.8534 ymax: 69.32212
Geodetic CRS: WGS 84The data contains a lot (111) of fields in addition to the geometry:
names(fires) [1] "poly_Sourc" "poly_Incid" "poly_Featu" "poly_MapMe" "poly_GISAc"
[6] "poly_Creat" "poly_DateC" "poly_Polyg" "poly_IRWIN" "poly_FORID"
[11] "poly_Acres" "poly_Sou_1" "poly_Sou_2" "attr_Sourc" "attr_ABCDM"
[16] "attr_ADSPe" "attr_Calcu" "attr_Conta" "attr_Contr" "attr_Creat"
[21] "attr_Incid" "attr_Disco" "attr_Dispa" "attr_Estim" "attr_Final"
[26] "attr_FFRep" "attr_FFR_1" "attr_FFR_2" "attr_FireB" "attr_Fir_1"
[31] "attr_Fir_2" "attr_Fir_3" "attr_FireC" "attr_Fir_4" "attr_Fir_5"
[36] "attr_Fir_6" "attr_FireD" "attr_Fir_7" "attr_FireM" "attr_FireO"
[41] "attr_FireS" "attr_Fir_8" "attr_Fir_9" "attr_Fi_10" "attr_FSJob"
[46] "attr_FSOve" "attr_GACC" "attr_ICS20" "attr_ICS_1" "attr_ICS_2"
[51] "attr_ICS_3" "attr_Inc_1" "attr_Inc_2" "attr_Inc_3" "attr_Inc_4"
[56] "attr_Inc_5" "attr_Initi" "attr_Ini_1" "attr_Ini_2" "attr_Ini_3"
[61] "attr_Irwin" "attr_IsFir" "attr_IsF_1" "attr_IsFSA" "attr_IsMul"
[66] "attr_IsRei" "attr_IsTre" "attr_IsUni" "attr_Local" "attr_Modif"
[71] "attr_Perce" "attr_Per_1" "attr_POOCi" "attr_POOCo" "attr_POODi"
[76] "attr_POOFi" "attr_POOJu" "attr_POO_1" "attr_POO_2" "attr_POOLa"
[81] "attr_POO_3" "attr_POOLe" "attr_POO_4" "attr_POO_5" "attr_POO_6"
[86] "attr_POO_7" "attr_POO_8" "attr_POOPr" "attr_POO_9" "attr_PO_10"
[91] "attr_POOSt" "attr_Predo" "attr_Pre_1" "attr_Prima" "attr_Secon"
[96] "attr_Total" "attr_Uniqu" "attr_FORID" "attr_WFDSS" "attr_Est_1"
[101] "attr_Organ" "attr_Strat" "attr_Cre_1" "attr_Mod_1" "attr_Sou_1"
[106] "attr_IsCpx" "attr_CpxNa" "attr_CpxID" "attr_Sou_2" "GlobalID"
[111] "attr_Inc_6" "geometry" The data we read with st_read() is actually a sf and data.frame object, which we can see if we look at its class. sf stands for Simple Features, a standard for storing spatial features like points, lines, and polygons. The spatial features are stored in a column of the data frame named geometry.
class(fires)[1] "sf" "data.frame"Since it is a data frame we can filter the data using standard methods like {dplyr}. I filtered the data to fires in Colorado using the attr_POOOSt column (I assume that the POOSt stands for Point of Origin State). I also limited the data to fires greater than 200 acres in size.
If we want to see the top 20 fires by size (acres), we can filter and arrange the data frame using {dplyr}. The 9th fire on the list is the Alexander Mountain Fire, one of the several fires currently burning.
The st_drop_geometry() function from {sf} allows you to drop the geometry column from a sf/data.frame object so you can work with a simple data frame that is smaller in size.
co_fires |>
sf::st_drop_geometry() |>
slice_max(order_by = poly_Acres, n = 20) |>
ggplot(aes(poly_Acres, y = reorder(poly_Incid, poly_Acres))) +
geom_col(fill = "red", color = "red", width = 0.8) +
labs(x = "Size [Acres]",
y = NULL,
title = "20 Largest Wildfires in Colorado",
caption = "Data source: WFIGS Interagency Fire Perimeters") +
scale_x_continuous(labels = scales::comma) + # add commas to axes labels
theme_minimal()
The data can easily be visualized on an interactive map using the {leaflet} package (Cheng et al. 2024). You can zoom the map in/out and drag it around, and hovering over a polygon displays the name of the wildfire.
leaflet() |>
addTiles() |>
addPolygons(data = co_fires, label = ~poly_Incid, color = "red", weight = 2) Interactive map of wildfires in Colorado. Source - HIFLD Open data
Finally I’ll isolate the Quarry fire burning southwest of Denver and show how you can easily jazz up the map with some nice features from the {leaflet} and {leaflet.extras} packages.
Some nice features you can add include:
# isolate the quarry fire
quarry <- co_fires |> filter(poly_Incid == "Quarry")
leaflet() |>
addProviderTiles(providers$Esri.WorldImagery, group = "Esri.WorldImagery") |>
addProviderTiles(providers$OpenTopoMap, group = "OpenTopoMap") |>
addPolygons(data = quarry, color = "red", weight = 3, label = "Fire Perimeter") |>
addMarkers(lat = quarry$attr_Initi, lng = quarry$attr_Ini_1, label = "Start Location") |>
leaflet::setView(lat = quarry$attr_Initi, lng = quarry$attr_Ini_1, zoom = 13) |> # set initial view
leaflet::addScaleBar() |>
leaflet.extras::addResetMapButton() |>
leaflet::addMiniMap(
position = "topright",
tiles = providers$Esri.WorldStreetMap,
toggleDisplay = TRUE,
minimized = FALSE,
zoomLevelFixed = 8
) |>
addLayersControl(
baseGroups = c("Esri.WorldImagery", "OpenTopoMap"),
# toggle for layers on the topleft
position = "topleft")Interactive map of Quarry fire in Jefferson County, CO. Inset shows location of the main map relative to Denver.
I hope this post was interesting and gives you an idea of some ways you can get started working with and visualizing spatial data in R.
R version 4.4.1 (2024-06-14)
Platform: x86_64-apple-darwin20
Running under: macOS Sonoma 14.6.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/Denver
tzcode source: internal
attached base packages:
[1] stats graphics grDevices datasets utils methods base
other attached packages:
[1] leaflet.extras_2.0.0 leaflet_2.2.2 lubridate_1.9.3
[4] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[7] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1
[10] tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0
[13] sf_1.0-17
loaded via a namespace (and not attached):
[1] utf8_1.2.4 generics_0.1.3 renv_1.0.4
[4] class_7.3-22 KernSmooth_2.23-24 stringi_1.8.4
[7] hms_1.1.3 digest_0.6.36 magrittr_2.0.3
[10] timechange_0.3.0 evaluate_0.24.0 grid_4.4.1
[13] fastmap_1.2.0 jsonlite_1.8.8 e1071_1.7-14
[16] DBI_1.2.2 fansi_1.0.6 crosstalk_1.2.1
[19] scales_1.3.0 jquerylib_0.1.4 cli_3.6.3
[22] rlang_1.1.4 units_0.8-5 munsell_0.5.1
[25] withr_3.0.1 yaml_2.3.10 tools_4.4.1
[28] tzdb_0.4.0 colorspace_2.1-0 vctrs_0.6.5
[31] R6_2.5.1 proxy_0.4-27 lifecycle_1.0.4
[34] classInt_0.4-10 leaflet.providers_2.0.0 htmlwidgets_1.6.4
[37] pkgconfig_2.0.3 pillar_1.9.0 gtable_0.3.5
[40] glue_1.7.0 Rcpp_1.0.13 xfun_0.46
[43] tidyselect_1.2.1 rstudioapi_0.16.0 knitr_1.48
[46] farver_2.1.1 htmltools_0.5.8.1 labeling_0.4.3
[49] rmarkdown_2.27 compiler_4.4.1