Quantifying the local impact of wildfire on demographic factors

Harrison DeFord

Introduction

Why this is important

  • Since the advent of aircraft, land management policy has been to put out wildfires immediately, even in the face of ecosystem needs
    • Protection of people, property, and prosperity (Westover 2023)
    • Species such as lodgepole pine (P. contorta) and wild lupine (L. perennis) require fire

Why this is important

  • Davies et al. (2018) found that majority Black, Hispanic, and Native American areas were disproportionately impacted by wildfire
    • About twice as vulnerable to wildfire impacts when compared to other census tracts

The Calf Canyon–Hermits Peak Complex Fire, as seen from Cimarron, N.M., June 15, 2022. The fire, the largest in state history, impacted the historically Latine communities of Las Vegas, Rociada, and Mora.

Initial goals

My initial goals for this project were to:

  • Identify sources of spatial data related to historic wildfire perimeters
  • Use ACS data to estimate impacts of disasters between Decennial Census years
  • Create functions generic enough to apply to any disaster (or other) boundary

Data sources

Wildfire perimeters

Historic wildfire GIS information is maintained by the National Interagency Fire Center (NIFC).

NIFC data are distributed through ArcGIS Online products:

  • Hosted feature layers
  • Esri Hub (“Open Data”) site

These sources can be queried using the {arcgislayers} package (Parry 2024).

Census information

The functions I’ve written for this project lean heavily on Kyle Walker’s {tigris} (2024) and {tidycensus} (2024) packages.

  • Increase performance by only downloading required counties and tracts
  • Provide an easy spatial filter for NIFC data, by using the filter_geom argument in arc_read()

Gathering and cleaning wildfire data

Define the New Mexico boundary

nm <- tigris::states() |>
  filter(STUSPS == "NM")

Get the wildfire perimeters

nm_fires <- arcgislayers::arc_read(r"(https://services3.arcgis.com/T4QMspbfLg3qTGWY/ArcGIS/rest/services/InterAgencyFirePerimeterHistory_All_Years_View/FeatureServer/0)",
  where = "FIRE_YEAR BETWEEN 2010 AND 2020",
  filter_geom = st_as_sfc(nm)
) |>
  st_make_valid(geos_method = "valid_linework") |>
  distinct(toupper(INCIDENT), .keep_all = T) |>
  slice_max(
    order_by = Shape__Area,
    n = 10
  ) |>
  select(-`toupper(INCIDENT)`) |>
  mutate(COMMENTS = str_trunc(COMMENTS, 35))

Note

The NIFC data is much less clean than the Census data; we have to validate geometry and de-duplicate incidents. In addition, many of the comment-writers had lots to say – we’ll truncate the COMMENTS field to 35 characters.

Put them together

The cumulative area of these fires is 3368.7 square miles. However, there are other impacts that are not apparent by only looking at the burned land area.

Relating demographic information to wildfire

Determining a temporal range

  • An analyst must define what is meant by “before” and “after” an event, particularly when using 5-year ACS estimates
    • Default setting is 2 years before, 3 years after
    • Can be changed using years_before and years_after parameters, but must still be >= 5 years

Finding intersecting counties and tracts

  • To calculate the difference in demographic characteristics before and after an event, we’ll use census tracts
    • Determine which tracts an event boundary intersects with before and after it occurs
    • get_counties_all_years() and get_tracts_all_years()
      • Use st_filter() to determine which tracts an event boundary intersects with

Relating tract data to events

  • Each event in the input sf object is spatially joined to the tracts it intersects with before and after
    • A set of unique state, county, and year combinations is generated to improve performance
    • map_data() is applied to iterate get_acs() across these combinations
      • Common variables as well as user-defined variables

Summarize ACS data and join to events

  • Most impactful events are large enough to affect more than one Census tract
    • Estimates are aggregated across event boundaries
      • Support for extensive and intensive variables is included

ACS Variables

Variable starts with… Function to be computed
mean mean()
median median()
All other strings sum()

Default variables include:

c(
      "pop_total" = "B03002_001", # total population
      "pop_nh_white" = "B03002_003", # non-Hispanic White
      "pop_nh_black" = "B03002_004", # non-Hispanic Black
      "pop_nh_aian" = "B03002_005", # non-Hispanic American Indian/Alaskan Native
      "pop_nh_asian" = "B03002_006", # non-Hispanic Asian
      "pop_nh_hipi" = "B03002_007", # non-Hispanic Native Hawaiian And Other Pacific Islander
      "pop_nh_other" = "B03002_008", # non-Hispanic Some Other Race Alone
      "pop_nh_two" = "B03002_009", # non-Hispanic Two or More Races
      "pop_hisp" = "B03002_012", # Hispanic or Latino
      "median_hhi" = "B19013_001", # median household income
      "mean_hh_size" = "B25010_001", # household size
      "sch_enrollment" = "B14001_002" # school enrollment
    )

Calculating difference

The final output object is an sf object, returned in either "tidy" (the default) or "wide" format.

output = "tidy" output = "wide"
Twice as long as nrow(data) * length(variables), with one field containing the point-in-time, one containing the corresponding ACS vintage year, one containing the variable, and one containing the estimate Same length as nrow(data), with a column returned for each variable requested and one denoting the years of the ACS vintages used
No differences are calculated; the user must group and perform calculations Absolute differences are calculated; this avoids divide by zero error
Can be used for side-by-side comparison maps (such as a slider map) comparing an event to itself Useful for single choropleth maps comparing events to one another

“tidy” output

nm_fires_impact <- get_incident_impacts(
  data = nm_fires,
  id_col = OBJECTID,
  year = FIRE_YEAR
)
OBJECTID MAP_METHOD DATE_CUR COMMENTS GEO_ID IRWINID UNQE_FIRE_ID FIRE_YEAR LOCAL_NUM INCIDENT UNIT_ID POO_RESP_I FEATURE_CA GIS_ACRES USER_NAME SOURCE AGENCY FIRE_YEAR_INT Shape__Area Shape__Length FORID OTHERID pt_in_time data_yr variable computed
24800 Digitized From Hardcopy Orthophoto quad 201202230000 NA {066DE785-7868-44DC-B1FF-AD0600D05F53} NA 2011-AZASF-000020 2011 2011-020 Wallow NA NA Wildfire Final Fire Perimeter 538051.6 NA USFS USFS 2011 3161860839 725941.9 NA NA before 2009 mean_hh_size 2.83
24800 Digitized From Hardcopy Orthophoto quad 201202230000 NA {066DE785-7868-44DC-B1FF-AD0600D05F53} NA 2011-AZASF-000020 2011 2011-020 Wallow NA NA Wildfire Final Fire Perimeter 538051.6 NA USFS USFS 2011 3161860839 725941.9 NA NA before 2009 median_hhi 33802.00
24800 Digitized From Hardcopy Orthophoto quad 201202230000 NA {066DE785-7868-44DC-B1FF-AD0600D05F53} NA 2011-AZASF-000020 2011 2011-020 Wallow NA NA Wildfire Final Fire Perimeter 538051.6 NA USFS USFS 2011 3161860839 725941.9 NA NA before 2009 pop_hisp 3446.00
24800 Digitized From Hardcopy Orthophoto quad 201202230000 NA {066DE785-7868-44DC-B1FF-AD0600D05F53} NA 2011-AZASF-000020 2011 2011-020 Wallow NA NA Wildfire Final Fire Perimeter 538051.6 NA USFS USFS 2011 3161860839 725941.9 NA NA before 2009 pop_nh_aian 4764.00
24800 Digitized From Hardcopy Orthophoto quad 201202230000 NA {066DE785-7868-44DC-B1FF-AD0600D05F53} NA 2011-AZASF-000020 2011 2011-020 Wallow NA NA Wildfire Final Fire Perimeter 538051.6 NA USFS USFS 2011 3161860839 725941.9 NA NA before 2009 pop_nh_asian 69.00

“wide” output

nm_fires_impact_wide <- get_incident_impacts(
  data = nm_fires,
  id_col = OBJECTID,
  year = FIRE_YEAR,
  output = "wide"
)
OBJECTID MAP_METHOD DATE_CUR COMMENTS GEO_ID IRWINID UNQE_FIRE_ID FIRE_YEAR LOCAL_NUM INCIDENT UNIT_ID POO_RESP_I FEATURE_CA GIS_ACRES USER_NAME SOURCE AGENCY FIRE_YEAR_INT Shape__Area Shape__Length FORID OTHERID acs_vintage mean_hh_size median_hhi pop_hisp pop_nh_aian pop_nh_asian pop_nh_black pop_nh_hipi pop_nh_other pop_nh_two pop_nh_white pop_total sch_enrollment
24800 Digitized From Hardcopy Orthophoto quad 201202230000 NA {066DE785-7868-44DC-B1FF-AD0600D05F53} NA 2011-AZASF-000020 2011 2011-020 Wallow NA NA Wildfire Final Fire Perimeter 538051.56 NA USFS USFS 2011 3161860839 725941.9 NA NA 2009 to 2014 0.1250000 7718 -434 392 -31 98 2 2 -101 -4266 -4338 -1479
29377 Infrared Interpretation 201212310000 NA {D542FB27-0045-4CF3-A49F-A5913F690B67} NA 2012-NMGNF-000200 2012 000200 Whitewater-Baldy Complex NMGNF NMGNF Wildfire Final Fire Perimeter 297801.12 NA USFS USFS 2012 1730419827 1116513.1 NA NA 2010 to 2015 0.1466667 -1097 -236 -101 12 20 0 0 50 -290 -545 215
13348 Other 201107210000 NA {7BF09AF5-9A22-4A7A-B64B-2C959495DAE4} NA 2011-NMSNF-000099 2011 000099 Las Conchas NMSNF NMSNF Wildfire Final Fire Perimeter 156593.62 NA USFS USFS 2011 967018464 440306.1 NA NA 2009 to 2014 -0.0175000 -4409 60 1068 -33 -3 19 108 -75 -27 1117 551
29401 Infrared Interpretation 201307112315 KCFAST Report Incorrectly Has Si... {DF4D9D65-A856-472F-A827-60C4E2268D61} NA 2013-NMGNF-000008 2013 000008 Silver NMGNF NMGNF Wildfire Final Fire Perimeter 138705.45 NA USFS USFS 2013 799345275 375737.0 NA NA 2011 to 2016 -0.1866667 -15806 52 -25 20 27 0 19 151 -960 -716 -266
46314 GPS - Uncorrected Data 201202010000 Data received from fire team member {77E85B4F-7FFB-44F5-B0DC-B95568BCDD16} NA 2011-NMLNF-000425 2011 000425 Donaldson NMN5S NMN5S Wildfire Final Fire Perimeter 96129.21 NA USFS USFS 2011 558682949 206076.9 NA NA 2009 to 2014 0.6400000 -167 737 12 -8 -26 -5 0 -114 1211 1807 670

And some maps!

Using output = "tidy" allows us to compare one variable easily over time.

And more maps!

Using output = "wide" allows us to compare differences in one variable easily across multiple events.

In Conclusion

Findings

The Wallow Fire (2011) had the most severe impact on the Hispanic and Latino populations based on this methodology.

  • Hispanic/Latino population decreased by 442 individuals between 2009 and 2014
  • Total population decreased by over 4,300 people
  • School enrollment decreased by 1,479 pupils.

Further applications

The way that I’ve constructed this project allows it to be easily extended to other topics, not just wildfire:

  • Other urban disasters, like flooding
  • Positive events (like a new transit station opening)
  • A novel polygon, to see change in a region over time

Difficulties and limitations

  • Data is returned by tigris with different column names depending on the vintage used
  • ACS estimates are only available from 2009, so impacts can only be calculated for events occurring in 2010 onwards
  • Calculated differences (when output = "wide") can only be absolute
    • If an initial state is 0, then a relative difference cannot be calculated

What the project does well

  • Written in such a way it can be converted to a package relatively easily
  • Provides a good counterpart to existing disaster planning tools
  • Does not calculate an impact index
    • What is a negative social characteristic?

One last example

The Washington Metropolitan Area Transit Authority (WMATA) opened a series of metro rail stations through Northern Virginia on July 26, 2014, in the Phase 1 extension of the Silver Line (“Dulles Corridor Metrorail Project” 2015). These stations include Greensboro, McLean, Spring Hill, Tysons, and Wiehle–Reston East.

Using the above methodology, we can investigate the question: what was the impact of opening these stations on the surrounding areas?

Using transit data

According to the U.S. Department of Transportation, most pedestrians are willing to walk a quarter to a half a mile to a transit stop (2008). We can compute a polygon for the upper bound of this estimate using the mapboxapi package.

wmata_stops <- st_read("data/wmata_silver_p1.gpkg", layer = "wmata_silver_p1")
hf_mi_stops <- mb_isochrone(wmata_stops, 
                            profile = "walking", 
                            distance = c(805), 
                            id_column = "stop_id") |>
  mutate(stop_yr = 2014) |>
  left_join(st_drop_geometry(select(wmata_stops, c(stop_id, stop_name))),
            by = join_by(id == stop_id))

Using transit data

Calculating differences over time

wmata_impact <- get_incident_impacts(hf_mi_stops,
                                     id,
                                     stop_yr,
                                     output = "tidy")

Calculating differences over time

wmata_impact_wide <- get_incident_impacts(hf_mi_stops,
                                     id,
                                     stop_yr,
                                     output = "wide")

References

Davies, Ian P., Ryan D. Haugo, James C. Robertson, and Phillip S. Levin. 2018. “The Unequal Vulnerability of Communities of Color to Wildfire.” Edited by Julia A. Jones. PLOS ONE 13 (11): e0205825. https://doi.org/10.1371/journal.pone.0205825.
“Dulles Corridor Metrorail Project.” 2015. https://www.transit.dot.gov/sites/fta.dot.gov/files/docs/2015-03_-_Dulles_Phase_1_-_Comprehensive_Monthly_Report.pdf.
Nabors, Dan, Robert Schneider, Dalia Leven, Kimberly Lieberman, and Colleen Mitchell. 2008. “Pedestrian Safety Guide for Transit Agencies.” safety.fhwa.dot.gov; U.S. Department of Transportation. https://safety.fhwa.dot.gov/ped_bike/ped_transit/ped_transguide/ch4.cfm.
Parry, Josiah. 2024. “Arcgislayers: An Interface to ArcGIS Data Services.” https://CRAN.R-project.org/package=arcgislayers.
Walker, Kyle. 2024. “Tigris: Load Census TIGER/Line Shapefiles.” https://CRAN.R-project.org/package=tigris.
Walker, Kyle, and Matt Herman. 2024. “Tidycensus: Load US Census Boundary and Attribute Data as ’Tidyverse’ and ’Sf’-Ready Data Frames.” https://CRAN.R-project.org/package=tidycensus.
Westover, Robert Hudson. 2023. “Only You | US Forest Service.” US Forest Service. https://www.fs.usda.gov/about-agency/features/only-you.