The world’s leading publication for data science, AI, and ML professionals.

🔥Wildfire Modeling in Yosemite National Park

This post will highlight the most important sections of our tutorial on modeling wildfire events.

Photo by Sheng Li on Unsplash
Photo by Sheng Li on Unsplash

This updated post (2024–06–02) will highlight the most important sections of our tutorial on modeling wildfire events using Google Earth Engine (GEE), geemap, GeoPandas, and GRASS GIS within a Jupyter notebook in the Google Colab platform. Through our tutorial, you can perform a fire simulation in Yosemite National Park California USA, during the summer season from 2020–06–20 through 2020–09–22. The tutorial has three main sections: Getting Data, Importing Data, and Wildfire Simulation.

This tutorial stems from the final work for __ the workshop on Geospatial Data processing and analysis using GRASS GIS software organized by Mario Gulich Institute for Advanced Space Studies (Argentina) and taught by Dr. Veronica Andreo.

To follow the tutorial, go to https://github.com/acoiman/wildfire_modeling, run the Jupyter Notebook in Google Colab by clicking the "Open in Colab" button.


1. Getting Data

Getting the data we need to feed our models could be daunting. Sometimes, data are not readily available because they are sparse or dispersed over various sources. In this section, we will show you an approach to gather data for wildfire modeling purposes readily. Table 1 shows the data used in the tutorial.

The tutorial is conceived in a fashion that it will not be required to install software or packages into your local machine. You only need to run it on the Google Colab platform.

1.1 Downloading Data from GEE

We will call the WDPA FeatureCollection (FC) and extract the Yosemite National Park polygon. We will save the extracted FC as a shapefile into our working directory. We also create an interactive map using using geemap package.

# get World Data Protected Areas (WDPA) FeatureCollection
pa = ee.FeatureCollection("WCMC/WDPA/current/polygons")

# Yosemite National Park polygon
filter = ee.Filter.inList('ORIG_NAME', ['Yosemite'])
yosemite = pa.filter(filter)

# define a function to filter polygons
def is_polygon(feature):
    # Return true if the feature's geometry type is 'Polygon' or 'MultiPolygon'
    geometry_type = feature.geometry().type()
    return ee.Algorithms.If(
        ee.Algorithms.IsEqual(geometry_type, 'Polygon'),
        True,
        ee.Algorithms.IsEqual(geometry_type, 'MultiPolygon')
    )

# apply the map method with a filter function to get only polygons
yosemite = yosemite.map(lambda feature: ee.Feature(feature).set('isPolygon', is_polygon(feature)))
                .filter(ee.Filter.eq('isPolygon', True))

# save Yosemite fc into shp
geemap.ee_to_shp(yosemite, filename="gisdata/shp/yosemite.shp")
Yosemite National Park with geemap package. Photo by the Author
Yosemite National Park with geemap package. Photo by the Author

To obtain 100-Hour Fuel Moisture, Wind Direction, and Wind Velocity data from GEE we need first to call the GRIDMET ImageCollection and filter by date. Second, we will create a rectangle that will be used to clip the data. Next, we will select the required data and calculate the mode. The resulting image will be clipped to the previously created polygon. After that, we will create 2000 random points that will be used to sample the Image using the mean. Finally, the resulting FC will be transformed into a GeoDataFrame and saved as a shapefile in our working directory.

# load GRIDMET ImageCollection
gridmet = ee.ImageCollection('IDAHO_EPSCOR/GRIDMET') 
    .filter(ee.Filter.date('2020-06-20', '2020-09-22'));

# create a rectangle
coords = [-120.10157709250639,37.36645735913506,-118.99122199244081,38.31280791348237]
bbox = ee.Geometry.Rectangle(coords)

# select 100-hour dead fuel moisture (fm100) and calculate the mode
fm100 = gridmet.select('fm100').mode()

# clip data to rectangle
fm100_clip =fm100.clip(bbox)

# create 2000 random points
fm100_points = ee.FeatureCollection.randomPoints(bbox, 2000);

# sample points using the mean
samplefm100 = fm100_clip.reduceRegions(**{
  'collection':fm100_points,
  'reducer':ee.Reducer.mean(),
  'scale': 4000,
  });

# download FeatureCollection as shp
out_dir = 'gisdata/shp/'
out_shp = os.path.join(out_dir, "samplefm100.shp")
geemap.ee_export_vector(samplefm100, out_shp, verbose=True)

We will download Enhanced vegetation index (EVI) data from GEE by calling the Landsat 8 Collection 1 Tier 1 32-Day EVI Composite, and filter by date and bounds using the Yosemite FC. Next, we will calculate a mean Image and clip it using the previously created polygon. Lastly, the clipped Image will be exported to our Google Drive storage, and from there we will move it to our working directory.

1.2. Downloading Fire Behavior Fuel Model

To run our fire simulation we need a fuel model defined by the USDA Forest Service. In this case, we will use the 13 Anderson Fire Behavior Fuel Model available at https://www.landfire.gov/fuel/fbfm13. This model has 13 classes of fuel which vary in load, distribution, and particle size (Petrasova, et.al., 2018).


2. Importing Data

Now, we will proceed to import the data we downloaded from GEE and other sources. This section involves setting up our GRASS GIS environment and importing the data itself.

2.1 Setting up GRASS GIS

When setting up GRASS in Jupyter Notebook we use Python to initialize GRASS GIS. First, we install GRASS GIS and then "ask" GRASS GIS where its Python packages are to be able to run it from the notebook. Next, we import GRASS GIS packages. Finally, we set up default font displays and the overwrite behavior so that we do not need to add the overwrite flag every time we execute a command.

# install GRASS GIS
%%bash
DEBIAN_FRONTEND=noninteractive
sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable
apt update
apt install grass subversion grass-dev
apt remove libproj22
# ask GRASS GIS where its Python packages are to be able to run it from the notebook
sys.path.append(
    subprocess.check_output(["grass", "--config", "python_path"], text=True).strip()
)
# import the GRASS GIS packages we need
import grass.script as gs
import grass.jupyter as gj

# default font displays
os.environ['GRASS_FONT'] = 'sans'
# overwrite existing maps
os.environ['GRASS_OVERWRITE'] = '1'
gs.set_raise_on_error(True)
gs.set_capture_stderr(True);

In GRASS GIS, projects are organized by areas stored in subdirectories called Locations. Each Location is defined by its coordinate system, projection, and geographical limits. To create a Location we will use the command grass with flags -e (exit after the creation of location or mapset), -c EPSG:3310 to create a new GRASS projected location in specified GISDBASE with given EPSG code (3310) (NAD83 / California Albers). We indicate to GRASS GIS that the initial database directory is grassdata and the Location name is wfs_3310, then we create a new Mapset called yosemite and check if we are working into this Mapset.

!grass -e -c EPSG:3310 grassdata/wfs_3310

# data directory
homedir = "/content/drive/MyDrive/WFS_updated"

# GRASS GIS database variables
grassdata = os.path.join(homedir, "grassdata")
location = "wfs_3310"
mapset = "PERMANENT" #"yosemite"

# Start the GRASS GIS Session
session = gj.init(grassdata, location, mapset)

# Create a new mapset
gs.run_command("g.mapset",
               mapset="yosemite",
               flags="c")
# show current GRASS GIS settings to check if we are in the created Location and Mapset
gs.gisenv()

2.2 Loading Data

In this subsection, we will upload all the datasets we need to run our fire simulation. First, we will import our raster containing the fuel model, then we will set the the geographic area in which GRASS should work (Region) to the fuel raster. Finally, we will proceed to import other required datasets.

# import fuel model raster
# your path could differ, so set you path accordingly
gs.run_command('r.import', input='gisdata/raster/LC20_F13_200.tif',output='fuel', 
    resolution='value', resolution_value=30.0);

# set region
g.region(raster='fuel')

# print region info
print(g.region(flags='p', stdout_=PIPE).outputs.stdout)
# import Yosemite National Park vector
gs.run_command('v.import', input='gisdata/shp/yosemite.shp',output='yosemite')

# import evi raster
gs.run_command('r.import', input='gisdata/raster/evi.tif', output='evi', 
    resolution='value', resolution_value=30.0)

# import fm 100h samples
gs.run_command('v.import', input='gisdata/shp/samplefm100.shp', 
                   output='samplefm100')

# import vs (wind velocity) samples
gs.run_command('v.import', input='gisdata/shp/samplevs.shp', 
                   output='samplevs')

# import th (wind direction) samples
gs.run_command('v.import', input='gisdata/shp/sampleth.shp', 
                   output='sampleth')

# import landsat image
gs.run_command('r.import', input='gisdata/raster/landsat.tif',output='landsat', 
                   resolution='value', resolution_value=30.0)

# import dem
gs.run_command('r.import', input='gisdata/raster/dem.tif',output='dem',
               resolution='value',resolution_value=30.0);

3. Wildfire Simulation

We will simulate wildfire spread in the Yosemite National Park. We will derive new data from the imported data, model fire events with different durations, and visualize wildfire propagation.

3.1 Creating a Map of the Study Area

To visualize the park polygon and the study area we will use the grass.jupyter package. Through the InteractiveMap() function we will create an interactive GRASS map. We first create an interactive map, set vector styles, and add vectors of the park and the study area.

# create interactive map
fig = gj.InteractiveMap(width=600, height=600, tiles='OpenStreetMap', use_region=True)

# set vector styles
def yose_style(vector):
  return {
      "fillOpacity": 0.0,
      "weight": 3,
      'color': "#000"
      }
def sa_style(vector):
  return {
      "fillOpacity": 0.0,
      "weight": 2,
      'color': "#FF0000"
      }

# add vectors and map control
fig.add_vector("yosemite", title='Yosemite Boundary', style_function=yose_style)
fig.add_vector("sa", title="Study Area",  style_function=sa_style)
fig.add_layer_control(position = "bottomright")

# add map title
display(HTML('<h1 style="font-size:150%";>
<strong>Study Area</strong></h1>'))

# show map
fig.show()
Study Area. Photo by the Author
Study Area. Photo by the Author

The study area is located in the center-east of the Yosemite National Park. The mean slope is 14.49 degrees, the mean wind speed is 467.37 feet/min, the mean dead fuel moisture 1 hour is 2.87%, and the area is 270.26 km². Predominant fuel classes are closed timber litter with a coverage of 29.70%, and timber (grass and understory) covering 21.67% of the study area. According to the NLCD: USGS National Land Cover Database, this area is mainly covered by evergreen forest (51%), shrub/scrub (34%), and grassland/herbaceous (5%).

Overview of the Study Area. Photo by https://tedmuller.us/Outdoor/Hiking/2015/150714-GaylorLakesGraniteLakesLoop.htm
Overview of the Study Area. Photo by https://tedmuller.us/Outdoor/Hiking/2015/150714-GaylorLakesGraniteLakesLoop.htm

The goal of the fire simulation in this Region is to determine how wildfire simulation behaves in a large area surrounded by mountains (valley) for a long time (32 hours).

3.2 Calculating Data

Modeling wildfire events is a two-step process. First, we execute the GRASS r.ros module that generates Rate of Spread (ROS) raster maps. Second, we run the GRASS r.spread module that simulates elliptically anisotropic spread.

As we assume windy and topographically sharp conditions, __ th_e r.ro_s module needs the following inputs:

  • moisture_1h: raster map containing the 1-hour fuel moisture (%).
  • moisture_live: raster map containing live fuel moisture (%).
  • velocity: raster map containing wind velocities (ft/min).
  • direction: raster map containing wind directions (degree).
  • slope: raster map containing slope (degree).
  • aspect: raster map containing aspect, counterclockwise from east (GRASS convention) in degrees.
  • elevation: raster map containing elevation (m).

To calculate the required inputs we will create and apply a Python function.

def caldata(regname, suffix):

    '''Function to calculate GRASS GIS raster for slope, aspect, fuel moisture 100h, fuel moisture 10h,
    fuel moisture 1h, wind velocity, wind direction, and live fuel moisture

    :param:str regname:region name
    :param:str suffix:GRASS raster filename suffix

    :return str: successfulness message
    '''
    try:
        # set region
        g.region(region=regname)

        # calculate slope and aspect
        gs.run_command('r.slope.aspect', elevation='dem', slope='slope'+suffix, aspect='aspect'+suffix)

        # calculate fuel moisture 100h raster from point data
        gs.run_command('v.surf.idw', input='samplefm100', output='moisture_100h'+suffix, 
                       column='mean')

        # calculate wind velocity raster from point data
        gs.run_command('v.surf.idw', input='samplevs', output='wind_speed'+suffix, 
                       column='vsfpm')

        # calculate wind direction raster from point data
        gs.run_command('v.surf.idw', input='sampleth', output='wind_dir'+suffix, 
                           column='mean')

        # define variable for string concatanation
        ss1 = 'moisture_1h'
        ss10 = 'moisture_10h'
        ss100 = 'moisture_100h'
        lfm='lfm'

        # moisture_100h = moisture_10h + 1 = moisture_1h + 2
        # moisture_10h = moisture_100h -1
        expm10 = ss10+suffix+'='+ss100+suffix+'-1'
        r.mapcalc(expression=expm10)

        # moisture_100h = moisture_10h + 1 = moisture_1h + 2
        # moisture_1h =  moisture_10h -2
        expm1 = ss1+suffix+'='+ss10+suffix+'-2'
        r.mapcalc(expression=expm1)

        # estimating live fuel moisture from evi
        explfm = lfm+suffix+'=(417.602 * evi) + 6.78061'
        r.mapcalc(expression=explfm)

        # rescale LFM to 0-100
        output = lfm+suffix+'_scaled'
        r.rescale(input='lfm'+suffix, output=output, to=(0,100));

        return "successfully calculated"

    except:
        print("Something went wrong")

3.3 Modeling

Now, we will run the r.ros module for our Study Area.

# generates rate of spread raster map
r.ros(model='fuel', moisture_1h='moisture_1h_sa', 
    moisture_live='lfm_sa_scaled', velocity='wind_speed_sa', 
    direction='wind_dir_sa', slope='slope_sa', aspect='aspect_sa', 
    elevation='dem', base_ros='out_base_ros', 
    max_ros='out_max_ros', direction_ros='out_dir_ros', 
    spotting_distance='out_spotting');

Next, we will create a point raster map that will be used as a source of our fire simulation.

# create a vector map from an ASCII points
gs.run_command('v.in.ascii',input='gisdata/txt/source.txt', 
                   output='source', separator='comma');

# rasterize vector map
gs.run_command('v.to.rast', input='source', output='source', type='point', use='cat');

Through the r.spread module, we will simulate a 32-hour wildfire event. Our first simulation period is 480 min (8 hours) because the initial time for the current simulation is 0 (default value) and our lag is 480 min. Note that we use the "s" flag because we are considering spotting distance.

# elliptically anisotropic spread simulation 8 hours
r.spread(flags="s", base_ros='out_dir_ros', max_ros='out_max_ros', 
    direction_ros='out_dir_ros', start='source', 
    spotting_distance='out_spotting', wind_speed='wind_speed_sa', 
    fuel_moisture='moisture_1h_sa', output='spread_8h', lag=480);

Having executed the r.ros and r.spread modules we can now visualize our results. In this case, we will use the grass.jupyter package to create four maps. Each map will contain the rate of spread raster for each time lag, a vector indicating the source of the simulation, a CIR (Color Infrared) composite of the Landsat image, and the corresponding title, and legend.

# create interactive map
fig = gj.InteractiveMap(width=600, height=350, tiles='OpenStreetMap', use_region=True)

# add fire source
fig.add_vector("source", title="Fire source")

# add CIR image
fig.add_raster("landsat_cir_sa", opacity=0.8, title="Landsat CIR")

# ad raster cumulative time of spread 8 hours hours
fig.add_raster("spread_8h", opacity=1.0, title="Fire spread 8 hours")

# add layer control
fig.add_layer_control(position = "bottomright")
Wildfire spread within 8 and 16 hours. Photo by the Author
Wildfire spread within 8 and 16 hours. Photo by the Author
Wildfire spread within 24 and 32 hours. Photo by the Author
Wildfire spread within 24 and 32 hours. Photo by the Author

4. Conclusions

This tutorial demonstrates we can integrate GRASS GIS, Google Earth Engine, and Jupyter Notebook in a unique development environment. However, more work should be done to pass GEE and Geopandas objects into GRASS GIS modules.

The results of the Fire simulation in our Study Area suggest that the modeling raster outputs do not overlap non-vegetated areas. In our case, fire simulations perform better in large areas because our base spatial resolution is 30m which is suitable for studies at a regional scale. Data should be generated to simulate wildfire at a finer scale.

In this tutorial, for the v.surf.idw module we used default parameters. It could be interesting to try out different interpolation parameters and evaluate how they affect simulation outputs.

Finally, we suggest assessing fire spread outputs after running the r.ros and r.sread modules with field measurements of dead fuel moisture, live fuel moisture, wind velocity, and wind direction.

As always, any and all feedback related to this is welcome!


5. References

Petrasova, A., Harmon, B., Petras, V., Tabrizian, P., & Mitasova, H. (2018). Wildfire spread simulation. In Tangible Modeling with Open Source GIS (pp. 155–163). Springer, Cham.

Yang, L., Jin, S., Danielson, P., Homer, C., Gass, L., Case, A., Costello, C., Dewitz, J., Fry, J., Funk, M., Grannemann, B., Rigge, M. and G. Xian. 2018, A New Generation of the United States National Land Cover Database: Requirements, Research Priorities, Design, and Implementation Strategies, p. 108–123.


Related Articles