MODIS fire data in GRASS GIS
In this post, I’ll show how to calculate the yearly fire frequency based on the Moderate Resolution Imaging Spectroradiometer (MODIS) MYD14A2 and MOD14A2 products1–3. As input, the 8-day composite products, available from 2002 to February 2022, will be used.
Calculating the fire frequencies is easy enough. The tricky bit, as is often the case, is to get all the data in order. You’ll first have to import 900+ MYD14A2 and MOD14A2 images into a GRASS GIS database4, and then register them in a time-space data-set5,6. Mind you, you could do without the last step. However, it makes life a lot easier, especially if at a later stage you want to do more temporal analyses.
Importing and registring all these data layers by hand is a bit much. Luckily, with just a bit of Python code, this can be automated. I used the code in this post to create a fire frequency map for the TAL, a region covering the lowlands of Nepal and an area in the north of India. But you can use any region you wish. Just be aware that it may take a bit of time to download the data.
Preparations
Get the data
You can get the data from the AρρEEARS portal from the NASA7. The site offers a simple and efficient way to access geospatial data from a variety of federal data archives, including MODIS satellite data. The first step is to select your area of interest and date range through a simple menu.
You can opt to download the layers in the native projection (MODIS Sinusoidal projection) or in several other projections. In addition, you have the option to the layers in Geotifs or NetCDF-4 format. In the example below, I assume you have downloaded the data as Geotifs. I’ll leave the choice of the projection to you.
After selecting your region of interest and the product, you can batch download the selected layers directly. See this help page for other download options. Note that for this tutorial, you only need the FireMask layers. If you want to import the QA layers as well, you’ll need to adapt the code in this post.
Create a location/mapset
To work in GRASS GIS you first need to create a GRASS GIS database that includes a location with the same projection as the data layers you downloaded earlier. You can use the Location Wizard to create one. To open the wizard, use the ![[Add Location]](location-add.png)
button. Alternatively, in the menu, go to settings → grass working environment → create new location.
The wizard guides you through a series of dialogs to browse and select or define the appropriate projections. These steps are fairly self-explanatory. Just make sure that in the second step (select the CRS), you select the option to read the CRS from a georeferenced data file. On the next screen, you can select one of the Geotifs you just downloaded. Make sure to select No when asked if you want to import the file.
Import the data
Libraries
GRASS GIS has advanced temporal data processing functions. To use them, you need to register the layers in a grass gis space-time raster dataset (strds). If you are not familiar with the temporal data processing framework, I suggest you first read this overview.
To run a Python script, open GRASS GIS’s own Python editor, or start up your favorite editor from the GRASS GIS command line. Now, the first step is to import the required libraries. Next, define the working directory with the downloaded Geotif files.
The Firemask layers have 9 categories, including a number of not processed, non-fire and fire categories. These categories are described in the metadata. The code below defines the labels of each of these categories. These will later be used to define the labels of the imported layers using the r.category function. Obviously, this step is only useful if you plan to keep the imported layers.
# Categories satellite data
categories = """\
0|Not processed (missing input data)
1|Not processed (obsolete; not used since Collection 1)
2|Not processed (other reason)
3|Non-fire water pixel
4|Cloud (land or water)
5|Non-fire land pixel
6|Unknown (land or water)
7|Fire (low confidence, land or water)
8|Fire (nominal confidence, land or water)
9|Fire (high confidence, land or water)
"""Optionally, you can prepare the color definition rules. These will later be used to define the colors scheme of the imported layers. Again, this is only useful if you plan to keep the imported layers.
Import MOD14A data
To import the MOD14A data layers, you need to create a list with all the file names of the MOD14A layers. A convenient way to get such a list is by using the pathname pattern expansion offered by the glob function.
You may have noticed that the file names include the string doy followed by 7 digits. The number is an ordinal date (also called the Julian date). The first 4 digits represent the year, and the last 3 digits represent the day of the year. I am actually not sure whether the date represents the start, middle or end date of the 8-day period, so I am assuming here that it represents the end date.
The first part of the for loop uses the strptime function in the Python datetime library to extract the end date (image_end) of the 8-day period represented by the image. It uses this to calculate the start date (image_start) and time period (period). You’ll need the latter as input in the r.timestamp function to set the timestamp of the raster layer.
The second part of the for loop imports the Geotif in GRASS GIS, defines the categories, writes the metadata, and defines the color table. In addition, the name of the layer is written to the list_of_layers list. This list is later used to define the layers that need to be registered in the space-time raster dataset (strds).
# Import layers in grass
list_of_layers = []
for i, input_layer in enumerate(list_files):
# time period image
doy = (list_files[i]
.replace("./MOD14A2.061_FireMask_doy", "")
.replace("_aid0001.tif", ""))
image_end = datetime.datetime.strptime(doy, '%Y%j')
image_start = (image_end - relativedelta(days=7)).strftime("%-d %b %Y")
image_end = ((image_end + relativedelta(hours=23, minutes=59))\
.strftime("%-d %b %Y"))
period = f"{image_start}/{image_end} 23:59:59"
# import
layer_name = "{}_{}".format("MOD14A", doy)
list_of_layers.append(layer_name)
Module("r.in.gdal", flags=["r"], input=input_layer, output=layer_name,
overwrite=True)
Module("r.category", map=layer_name, separator="pipe", rules="-",
stdin_=categories)
layer_title=("Confidence of fire {}:{}"
.format(str(image_start), str(image_end)))
Module("r.support", map=layer_name, title = layer_title)
Module("r.colors", map=layer_name, rules="-", stdin_=color_rules)
Module("r.timestamp", map=layer_name, date=period)After the layers are imported (this may take some time), you define the space-time raster dataset MOD14A.
Now you can register the layers in the strds using the t.register function. It uses the list with the layers you created in the loop above as input.
Note that this function has the option to directly define the start and end time. So, you could have made this step part of the loop above. However, I have the (not properly tested) impression that doing this step separately is faster. But let me know in the comment section below if I am wrong.
Import MYD14A data
Repeat the steps above to import the MYD14A layers and register them in a strds.
# Get a list with the MOD14a files
list_files = glob.glob(pathname = "./MYD*FireMask*", recursive=True)
list_files.sort()
# Import layers in grass
list_of_layers = []
for i, input_layer in enumerate(list_files):
# dates
doy = (list_files[i]
.replace("./MYD14A2.061_FireMask_doy", "")
.replace("_aid0001.tif", ""))
image_end = datetime.datetime.strptime(doy, '%Y%j')
image_start = (image_end - relativedelta(days=7)).strftime("%-d %b %Y")
image_end = ((image_end + relativedelta(hours=23, minutes=59))\
.strftime("%-d %b %Y"))
period = f"{image_start}/{image_end} 23:59:59"
# import
layer_name = "{}_{}".format("MOD14A", doy)
list_of_layers.append(layer_name)
Module("r.in.gdal", flags=["r"], input=input_layer,
output=layer_name, overwrite=True)
Module("r.category", map=layer_name, separator="pipe", rules="-",
stdin_=categories)
layer_title=("Confidence of fire {}:{}"
.format(str(image_start), str(image_end)))
Module("r.support", map=layer_name, title = layer_title)
Module("r.colors", map=layer_name, rules="-",
stdin_=color_rules)
Module("r.timestamp", map=layer_name, date=period)
# Create the timespace data-set
Module("t.create", output="MYD14A", semantictype="sum",
title='8-day fire occurrences',
description='Confidence of fire 8-day interval',
temporaltype="absolute")
# Register the layers in the time-space data set
Module("t.register",
type = "raster",
input=timespacedata,
maps = list_of_layers,
start = image_start,
end = image_end) Reclassify the layers
Now that you have all the data layers, you need to reclassify them to -99 (no data), 0 (no fire) and 1 (fire). In the code below, the low, nominal, and high confidence fires are all three reclassified as fire. If you want to go for a more conservative estimate, you can opt to reclassify category 7 or 7 and 8 as no fire.
You’ll reclassify the layers MOD14A and MYD14A according to the rules defined above. But first, define the region to match that of the imported layers.
Reclassify MOD14A layers
Get a list of the MOD14A layers using g.list. As an alternative, you can also use the t.rast.list to get the list of raster layers in the strds.
Now, you can reclassify the raster layers according to the rules you defined earlier. In the example, the reclassified layers get the name of the original layers with the suffix _recl. Note that in the last two lines, the timestamp of the original layers is copied to the reclassified layer.
Create a new strds MOD14A_reclass79 and register the layers.
# Create spacetime for MOD14A reclass
Module("t.create", output="MOD14A_reclass79",
semantictype="sum",
title='fire/no fire 8-day interval',
description='fire (cats 7-9) or no fire',
temporaltype="absolute")
# Add reclass layers to new strds
Module("t.register",
type = "raster",
input="MOD14A_reclass79",
maps = newlays)Reclassify MYD14A layers
Now, repeat the steps above to reclassify the MYD14A layers
# Get list of MYD14A layers
lays = Module("g.list", type="raster",
pattern="MYD14A_*",
stdout_=PIPE).outputs.stdout.split("\n")
lays = [_s for _s in lays if _s]
# Reclassify the layers and copy timestamp
for _,layname in enumerate(lays):
Module("r.reclass",
input=layname,
output="{}_recl".format(layname),
rules="-", stdin_=reclassrules)
timst = Module("r.timestamp", map=layname, stdout_=PIPE).outputs.stdout
Module("r.timestamp", map="{}_recl".format(layname), date=timst)
# Create spacetime for MYD14A reclass
Module("t.create", output="MYD14A_reclass79",
semantictype="sum",
title='fire/no fire 8-day interval',
description='fire (cats 7-9) or no fire (MYD14A)',
temporaltype="absolute")
# Add reclass layers to new strds
Module("t.register",
type = "raster",
input="MYD14A_reclass79",
maps = newlays)Fire frequencies
Select layers
You can visualize the temporal coverage of both data sets using the timeline tool. The results will show you that the MOD14A data set covers a longer period than the MYD14A dataset. It starts at an earlier date, but coverage of this earlier period is incomplete.
To have a uniform approach for the entire period, select the subset of MOD14A data layers covering the same time period as the MYD14A dataset. This can be done using the t.rast.extract function. The result is a new strds with the extracted data layers.
# Start time MYD14A
strt = Module("t.info", flags="g", input="MYD14A", stdout_=PIPE).outputs
strt = strt.stdout.split("\n")[8].split("=")[1]
# Extract all layers after the MYD14A start date
Module("t.rast.extract", input="MOD14A_reclass79",
where="start_time >= {}".format(strt),
output="terra", basename="terra")Combine the layers
Now combine for each 8-day period the mod14a and myd14a layers by computing for each grid cell the maximum value of the corresponding cells in the input grid layers. This means that in the output layers, a fire occurrence will be recorded if the category value of at least one of the two corresponding input layers is 7 or higher.
If you want to be more cautious, and only accept a fire occurrence if recorded by both the MOD14A and MYD14A layers, you can use the minimum instead of the maximum value of the two layers.
Next, convert all raster cells with a value of -99 to NULL (no data). This can be done for all layers in once using the t.rast.null function.
This step ensures that no-data values are ignored when computing the average fire frequency numbers. This may result in a biased estimate if fires are more (or less) likely to occur in periods with an high likelihood of cloud coverage. To identify areas where this may be an issue, you can calculate the number of 8-day periods with fire/no-fire observations as well.
Compute fire frequencies
You use the t.rast.series function to compute the average number of fires per 8-day period. In the code below, the number of layers with a value is counted as well. You can use this layer to check the number of 8-day periods without value (nodata).
The last step is to compute the estimated average number of fire occurrences per year. Note that the underlying assumption is that there will never be more than one fire per 8-day period. Or in other words, within a grid cell, multiple fires within an 8-day period will always be registered as one fire occurrence.
As a very last, optional, step, you can remove all the intermediate layers and strds.
Printing the map
To print out the map presented at the top of this post, I used the code below. And yes, a lot of code to create a map. But in case you are going to create many maps of the same area, a bit of coding will save you time in the longer run.
# Input raster
inputraster="Fire_YearlyAverage@ModisFire"
# Output settings
outputfile = "TAL_fireyearly.png"
width_image = 3600
title = "Average fire occurrences per year"
# Get width/height ratio of image
Module("g.region", raster=inputraster)
region_settings = gs.parse_command("g.region", flags="g")
height = float(region_settings['n']) - float(region_settings['s'])
width = float(region_settings['e']) - float(region_settings['w'])
height_image = width_image / (width/height)
# Set environmental variables
os.environ['GRASS_RENDER_IMMEDIATE'] = "cairo"
os.environ['GRASS_RENDER_FILE'] = outputfile
os.environ['GRASS_RENDER_HEIGHT'] = str(height_image)
os.environ['GRASS_RENDER_WIDTH'] = str(width_image)
os.environ['GRASS_RENDER_FILE_READ'] = "TRUE"
os.environ['GRASS_FONT'] = "DejaVuSansCondensed"
# Render image
Module("d.rast", map=inputraster, values=value_range)
Module("d.grid", size=1, color="190:190:190",
text_color="100:100:100", border_color="160:160:160",
fontsize=32, width=1)
Module("d.vect", map="ne_10m_admin_0", type="area",
color="0:0:0", fill_color="none", width=8)
Module("d.vect", map="ne_10m_admin_0", type="area",
color="255:255:255", fill_color="none", width=4)
Module("d.legend", flags=['t', 's', 'b'],
raster=inputraster, bgcolor="255:255:255",
title=title, title_fontsize=42, fontsize=42,
digits=digits, font="arial", at=[90,92,70,90],
range=legend_range, border_color="none")