A national forest’s fire history¶
Forest Service data from the Enterprise Data Warehouse, with georest.esri.edw
The Forest Service’s Enterprise Data Warehouse (EDW) publishes about 200 ArcGIS services: forest boundaries, roads, trails, timber harvests, and the national record of large fires from Monitoring Trends in Burn Severity (MTBS). This notebook puts two of them together to answer a simple question about the Arapaho and Roosevelt National Forests in Colorado, home of the Cameron Peak fire: how much has burned, and when?
You will:
search EDW for the services you need,
pull a forest boundary and the fires around it,
map them, and
chart the acres burned by year.
pip install georest plotly folium
import folium
import pandas as pd
import plotly.express as px
import plotly.io as pio
from georest.esri import edw
pd.set_option("display.max_colwidth", None) # show text in full, never cut with "..."
pio.renderers.default = "notebook_connected" # interactive here and on georest.io
1. Find the services¶
search_edw_services matches a keyword against service names and descriptions, so you do not
need to know EDW’s service names in advance.
pd.DataFrame(edw.search_edw_services("fire"))[["name", "theme", "description"]]
| name | theme | description | |
|---|---|---|---|
| 0 | EDW_AerialFireRetardantAvoidanceAreas_Aquatic_01 | inland_waters | This data depicts aquatic aerial fire retardant avoidance areas delivered as part of the 2011 Nationwide Aerial Application of Fire Retardant on… |
| 1 | EDW_AerialFireRetardantAvoidanceAreas_Terrestrial_01 | environment | This data depicts terrestrial aerial fire retardant avoidance areas delivered as part of the 2011 Nationwide Aerial Application of Fire Retardant on… |
| 2 | EDW_BurnedAreaEmergencyResponse_01 | environment | An initial Burned Area Reflectance Classification (BARC) dataset is created by analyzing a differenced Normalized Burn Ratio (dNBR) image, which was… |
| 3 | EDW_CommWildfireDefenseGrant_01 | environment | Community Wildfire Defense Grants (CWDG) is a program administered by the Forest Service to help at-risk communities and Tribes plan for and reduce… |
| 4 | EDW_FireOccurrenceAndPerimeter_01 | environment | The FirePerimeterFinal polygon layer represents final mapped wildland fire perimeters. |
| 5 | EDW_FireOccurrenceCurrentEdition_01 | environment | EDW_FireOccurrenceCurrentEdition_01 provides the current edition of Karen C. Short’s Spatial wildfire occurrence data for the United States (FPA FOD).… |
| 6 | EDW_FireOccurrenceFIRESTAT_YRLY_01 | environment | The FIRESTAT (Fire Statistics System) Fire Occurrence point layer represents ignition points, or points of origin, from which individual wildland… |
| 7 | EDW_FireshedRegistry_01 | environment | The Fireshed Registry is a geospatial dashboard and decision tool built to organize information about wildfire transmission to communities and… |
| 8 | EDW_ForeststoFaucets_02 | inland_waters | Forests to Faucets 2.0 builds upon the national Forests to Faucets(2011) by updating base data and adding new threats including wildfire, invasive… |
| 9 | EDW_MTBS_01 | environment | The Monitoring Trends in Burn Severity (MTBS) Program assesses the frequency, extent, and magnitude (size and severity) of all large wildland fires… |
| 10 | EDW_RAVG_v2_01 | environment | The USDA Forest Service Rapid Assessment of Vegetation Condition after Wildfire (RAVG) program produces geospatial and related data representing… |
| 11 | RAVG_DataExtract_01 | uncategorized |
EDW_MTBS_01 holds MTBS. Its layers come per year, plus two that span every year:
layers = edw.get_service_info("EDW_MTBS_01")["layers"]
[(layer["id"], layer["name"]) for layer in layers if "All Years" in layer["name"]]
[(62, 'Fire Occurrence Locations (All Years)'),
(63, 'Burned Area Boundaries (All Years)')]
2. The forest, and the fires around it¶
query_features returns GeoJSON. The forest boundary comes from EDW_ForestSystemBoundaries_01;
then every MTBS burned area inside the forest’s bounding box, which takes in the forest and the
country right around it.
forest = edw.query_features("EDW_ForestSystemBoundaries_01", 1,
where="forestname = 'Arapaho and Roosevelt National Forests'",
out_fields="forestname,gis_acres")
outline = forest["features"][0]
print(f"{outline['properties']['forestname']}: {outline['properties']['gis_acres']:,.0f} acres")
def bounds(geometry):
xs, ys = [], []
def walk(coords):
if isinstance(coords[0], (int, float)):
xs.append(coords[0]); ys.append(coords[1])
else:
for c in coords:
walk(c)
walk(geometry["coordinates"])
return min(xs), min(ys), max(xs), max(ys)
west, south, east, north = bounds(outline["geometry"])
fires = edw.query_features_with_pagination(
"EDW_MTBS_01", 63, geometry=f"{west},{south},{east},{north}",
out_fields="fire_name,year,acres,fire_type,ig_date")
print(f"{len(fires['features'])} MTBS fires in and around the forest, "
f"{min(f['properties']['year'] for f in fires['features'])} to "
f"{max(f['properties']['year'] for f in fires['features'])}")
Arapaho and Roosevelt National Forests: 2,491,783 acres
47 MTBS fires in and around the forest, 1988 to 2024
3. The map¶
Perimeters arrive at survey detail, several megabytes for this many fires, which is more than a
web map needs. simplify below drops points that change a shape by less than about 150 meters
(the Douglas-Peucker method), shrinking the map to a fraction of the size with no visible
difference at this zoom.
Click a fire for its name, year and size. Darker red is more recent.
def simplify(geometry, tolerance=0.0015):
"""Thin a Polygon or MultiPolygon's rings (Douglas-Peucker, in degrees)."""
def thin(points):
if len(points) < 5:
return points
keep = [False] * len(points)
keep[0] = keep[-1] = True
stack = [(0, len(points) - 1)]
while stack:
a, b = stack.pop()
(ax, ay), (bx, by) = points[a][:2], points[b][:2]
dx, dy = bx - ax, by - ay
norm = (dx * dx + dy * dy) ** 0.5 or 1e-12
far, far_i = 0.0, None
for i in range(a + 1, b):
px, py = points[i][:2]
d = abs(dy * px - dx * py + bx * ay - by * ax) / norm
if d > far:
far, far_i = d, i
if far_i is not None and far > tolerance:
keep[far_i] = True
stack += [(a, far_i), (far_i, b)]
return [p for p, k in zip(points, keep) if k]
polygons = geometry["coordinates"] if geometry["type"] == "MultiPolygon" else [geometry["coordinates"]]
thinned = [[thin(ring) for ring in poly] for poly in polygons]
return {"type": "MultiPolygon", "coordinates": thinned}
years = [f["properties"]["year"] for f in fires["features"]]
first, last = min(years), max(years)
m = folium.Map(location=[(south + north) / 2, (west + east) / 2], zoom_start=8, tiles="https://server.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/{z}/{y}/{x}",
attr="Tiles © Esri")
folium.GeoJson({"type": "Feature", "geometry": simplify(outline["geometry"]), "properties": {}},
name="National forest",
style_function=lambda f: {"color": "#2f6f4f", "weight": 2, "fillOpacity": 0.08}).add_to(m)
for f in sorted(fires["features"], key=lambda f: f["properties"]["year"]):
p = f["properties"]
shade = (p["year"] - first) / max(last - first, 1) # 0 oldest .. 1 newest
color = f"rgb({int(250 - 60 * shade)}, {int(170 - 150 * shade)}, {int(90 - 80 * shade)})"
folium.GeoJson(
{"type": "Feature", "geometry": simplify(f["geometry"]), "properties": p},
style_function=lambda _, c=color: {"color": c, "weight": 1, "fillColor": c, "fillOpacity": 0.6},
tooltip=f"{p['fire_name'].title()} ({p['year']})",
popup=folium.Popup(f"<b>{p['fire_name'].title()}</b><br>{p['year']}, "
f"{p['acres']:,.0f} acres<br>{p['fire_type']}", max_width=220),
).add_to(m)
m
4. Acres burned, year by year¶
Most years see little. A few see most of the burning: hover over a bar for the year’s fires.
table = pd.DataFrame([f["properties"] for f in fires["features"]])
table["fire"] = table["fire_name"].str.title()
by_year = (table.groupby("year")
.agg(acres=("acres", "sum"), fires=("fire", lambda s: ", ".join(sorted(s))))
.reindex(range(first, last + 1)) # every year, including quiet ones
.reset_index())
by_year["acres"] = by_year["acres"].fillna(0)
by_year["fires"] = by_year["fires"].fillna("none")
fig = px.bar(by_year, x="year", y="acres", custom_data=["fires"],
labels={"acres": "Acres burned", "year": ""},
title="MTBS acres burned in and around the Arapaho and Roosevelt National Forests")
fig.update_traces(marker_color="#c2410c",
hovertemplate="<b>%{x}</b>: %{y:,.0f} acres<br>%{customdata[0]}<extra></extra>")
fig.show()
share = table.nlargest(1, "acres")["acres"].iloc[0] / table["acres"].sum()
print(f"The single largest fire accounts for {share:.0%} of all the acres burned here since {first}.")
The single largest fire accounts for 25% of all the acres burned here since 1988.
The ten largest:
(table.nlargest(10, "acres")[["fire", "year", "acres", "fire_type"]]
.style.format({"acres": "{:,.0f}"}).hide(axis="index"))
| fire | year | acres | fire_type |
|---|---|---|---|
| Cameron Peak | 2020 | 208,743 | Wildfire |
| East Troublesome | 2020 | 193,717 | Wildfire |
| Mullen | 2020 | 179,795 | Wildfire |
| High Park | 2012 | 90,698 | Wildfire |
| Last Chance | 2012 | 17,961 | Wildfire |
| Williams Fork | 2020 | 14,677 | Wildfire |
| Bobcat | 2000 | 10,437 | Wildfire |
| Calwood | 2020 | 9,832 | Wildfire |
| Alexander Mountain | 2024 | 9,691 | Wildfire |
| Picnic Rock | 2004 | 9,007 | Wildfire |
Where to go next¶
Another forest. Change
forestname;edw.query_features("EDW_ForestSystemBoundaries_01", 1, out_fields="forestname")lists them all.Burn severity. MTBS also maps how severely each fire burned. The imagery guide reads land-cover change rasters the same way.
Other EDW data.
edw.search_edw_services(theme="recreation"),"roads","timber"; andedw.query_features_analyticfor ranked summaries computed on the server.
Full reference: georest.esri.edw.