James Millington
library(terra)
library(rasterVis)
library(tidyverse)
library(googlesheets4)This document presents code and to implement the McArthur Fire model for a region of Cornwall, UK, parameterised using field data collected by King’s College London students in the Dept. of Geography. This is a Quarto document that combines text with code - here we use the R coding language, primarily with the terra package for working with raster data (see a brief tutorial to Raster GIS operations in R with terra here).
Equations for the McArthur fire model we will use come from Noble et al. (1980). To apply the model to a particular location and time we need to provide various numeric inputs.
In the field, students collected data to help estimate key variables for McArthur’s model for different land cover types represented in the UK CEH Land Cover Maps:
-
degree of curing, C, in %
-
fuel load, W, in tons per hectare
By combining these variables with climate variables and the UK CEH Land Cover Map we can estimate fire danger and rate of spread across Cornwall and for different climate scenarios.
Summaries of the data collected by the six groups of students are shown in the Summary tab of the Google spreadsheet used to record the data. In turn, these summaries depend on summary equations in each of the data sheets (see rows at the bottom of input data). Finally, calculations from the combined student field measurements are made in cells A26:E30 (shaded blue) in the Summary tab. Several of these calculations are based on the information presented in Pasalodos-Tato et al. (2015).
We can read these estimated values for C and W for the different land cover types directly our model code here:
gs4_deauth() #run in deauthorised mode as we are accessing a public sheet
sheet_data <-range_read(ss="https://docs.google.com/spreadsheets/d/1HReAFVljlUpjqmHOGYqvDZcUn9lkkPPRQ9ebb5BkQso/edit#gid=937321199",
sheet="Summary",
range="A26:E30")
sheet_data# A tibble: 4 × 5
`Land Cover` Grass Heath Decid Conif
<chr> <dbl> <dbl> <dbl> <dbl>
1 Curing 29.1 32.6 NA NA
2 W-ground 0.716 1.91 0.667 0
3 W-trees 0 6.28 105. 91.8
4 W-total 0.716 8.19 105. 91.8
First we read and display the UK CEH Land Cover Map for our study region. We will use this to determine the location of our fuel types across the study region.
fpath <- "/home/james/OneDrive/Teaching/2022-23/undergrad/fieldwork/cornwall/data/lcm-2020-25m_tiled_4639069/sx/sx06.tif"
sx06 <- rast(fpath, lyrs=1)
plot(sx06)
The colours and their corresponding values in the map above correspond to the UK CEH Land Cover classification.
We shall re-classify these for our fire model as follows:
| Fire Fuel Type | Fire Fuel Type ID | Original CEH Classes |
|---|---|---|
| Deciduous | 1 | 1 |
| Conifer | 2 | 2 |
| Arable | 3 | 3 |
| Grassland | 4 | 4-7 |
| Heathland | 5 | 9-10 |
| Non-Veg | 6 | 8, 11-21 |
Here’s the code for the re-classification:
#first create a table as above to define the re-classification
decid <- c(1,1,1)
conif <- c(2,2,2)
arable <- c(3,3,3)
grass <- c(4,7,4)
nonVeg <- c(8,8,6)
heath <- c(9,10,5)
nonVeg2 <- c(11,21,6)
FMC <- rbind(decid, conif, arable, grass, nonVeg, heath, nonVeg2)
#then run the reclassification on the original raster map
fuelType_int <- classify(sx06, rcl=FMC, right=NA)
#and quickly plot
plot(fuelType_int)
This looks like it has worked, but the map is not very informative, so let’s change colours and add a proper legend:
#first copy the map
fuelType_cat <- fuelType_int
#now define the classes represented in the map
cls <- data.frame(id=1:6,
cover=c("Deciduous", "Conifer", "Arable", "Grassland", "Heathland", "Non-Veg") )
levels(fuelType_cat) <- cls
#set the colours for each class
coltab(fuelType_cat) <- data.frame(values=1:6,
cols=c('forestgreen', 'darkolivegreen','wheat1', 'lightgreen', 'darkorchid', 'gray' ))
#and plot
plot(fuelType_cat, main="Fire Fuel Types")
This looks better. The next line of code will create an interactive map in R-Studio to show the location of the study area in Cornwall, but not in the non-interactive version of the document.
#if running this in RStudio we can visualise in an interactive map
#library(leaflet)
#leaflet() |> addTiles() |> addRasterImage(raster::raster(fuelType_cat), colors = "Spectral", opacity = 0.8)Now we create raster maps for the curing and biomass parameters, applying the estimates from the field to each of our land cover types. First, the curing input (from % dead vegetation) for Grassland and Heathland.
curing <- fuelType_int
curing[curing==1]<-NA
curing[curing==2]<-NA
curing[curing==3]<-NA
curing[curing==4]<-as.numeric(sheet_data[1,2]) #grassland
curing[curing==5]<-as.numeric(sheet_data[1,3]) #heathlad
curing[curing==6]<-NA
levelplot(curing, margin=F, main="Curing (% dead veg)",par.settings = RdBuTheme, at=seq(0, 100, by = 1), contour=T)
Next biomass for Deciduous forest , Conifer forest, Grassland and Heathland.
biomass <- fuelType_int
biomass[biomass==1]<-as.numeric(sheet_data[4,4]) #decid
biomass[biomass==2]<-as.numeric(sheet_data[4,5]) #conifer
biomass[biomass==3]<-NA
biomass[biomass==4]<-as.numeric(sheet_data[4,2]) #grassland
biomass[biomass==5]<-as.numeric(sheet_data[4,3]) #heathland
biomass[biomass==6]<-NA
levelplot(biomass, margin=F, main="Biomass (tons per ha)", par.settings=infernoTheme)
We also need to provide climate inputs - if we assume climate is uniform across the study region we don’t need to make raster maps but can just create individual variables:
temperature <- 20 #degrees C
humidity <- 50 #%
wind <- 10 #km/hr
KDBI <- 50 #KDBI (range is 0 [no drought] to 800 [extreme drought])
daysRain <- 1 #days since rain
pptn <- 10 #last rain in mmLater we can vary the values above to examine how different scenarios of weather influence fire danger and rate of spread of fire.
Now we can calculate Fire Danger using the McArthur grassland fire danger model for Grassland and Heathland fuel types using their respective curing and biomass parameter values, and the forest fire danger model for the Deciduous and Conifer types (based solely on climate inputs).
#see Noble et al. 1980 for equations
drought <- 0.191*(KDBI+104)*(daysRain+1)^1.5/(3.52*(daysRain+1)^1.5+pptn-1)
GrassM <- (((97.7+4.06*humidity)/(temperature+6)-0.00854*humidity))*((100-curing)/100)
GrassFDM <- 3.35*biomass*exp(-0.0897*GrassM+0.0403*wind)
ForestFDM <- 1.25*drought*exp((temperature-humidity)/30+0.0234*wind)And now apply the model calculations to the raster fuel map:
FDM <- GrassFDM
FDM[fuelType_int==1]<-ForestFDM #deciduous
FDM[fuelType_int==2]<-ForestFDM #conifer
FDM[fuelType_int==3]<-max(GrassFDM[fuelType_int==4])/2 #half grassland danger
FDM[fuelType_int==6]<-0 #non-veg
levelplot(FDM, margin=F, main="Fire Danger (McArthur scale)", par.settings=infernoTheme)
So here’s our first Fire Danger map! Note the scale on the right - the McArthur model is unitless but can range from 0 (no fire danger) into the 100s (high danger). Check you can see how the different fuel types have different danger levels. What land cover (fuel type) are are the light coloured areas?
We can also calculate predicted Rate of Spread (RoS) once a fire is alight. For forest fuel types this is dependent on a combination of fire danger and biomass, whereas for fuel types based on the grass model (Grassland, Heathland, Arable) RoS is related only to fire danger.
RoS <- FDM
#again see Noble et al. (1980) for equations
forestMulti<-max(FDM[fuelType_int==1])*0.0012
grassRoS<-max(FDM[fuelType_int==4]*0.13)
heathRoS<-max(FDM[fuelType_int==5]*0.13)
RoS[fuelType_int==1]<-biomass*forestMulti
RoS[fuelType_int==2]<-biomass*forestMulti
RoS[fuelType_int==3]<-grassRoS
RoS[fuelType_int==4]<-heathRoS
RoS[fuelType_int==5]<-grassRoS
RoS[fuelType_int==6]<-0
levelplot(RoS, margin=F, main="Rate of Spread (km per hour)", par.settings=infernoTheme,at=seq(0,5,0.25))
We could now use this model to explore different scenarios, primarily of different climatic conditions but also for different possible land use/cover futures.
To explore the effect of climate, set values we are interested in:
temperature <- 39 #degrees C
humidity <- 10 #%
wind <- 30 #km/hr
KDBI <- 600 #KDBI (range is 0 [no drought] to 800 [extreme drought])
daysRain <- 30 #days since rain
pptn <- 10 #last rain in mmThese values should imply a higher fire danger than our original model application (check you understand why from the values of the variables).
Then we can calculate the fire danger for these new climate values:
#same equations as above but now the climate variables have different values
drought <- 0.191*(KDBI+104)*(daysRain+1)^1.5/(3.52*(daysRain+1)^1.5+pptn-1)
GrassM <- (((97.7+4.06*humidity)/(temperature+6)-0.00854*humidity))*((100-curing)/100)
GrassFDM <- 3.35*biomass*exp(-0.0897*GrassM+0.0403*wind)
ForestFDM <- 1.25*drought*exp((temperature-humidity)/30+0.0234*wind)And finally apply the model to our map of fuel types to estimate fire danger and compare to the original Fire Danger map.
FDM2 <- GrassFDM
FDM2[fuelType_int==1]<-ForestFDM #deciduous
FDM2[fuelType_int==2]<-ForestFDM #conifer
FDM2[fuelType_int==3]<-max(GrassFDM[fuelType_int==4])/2 #half grassland danger
FDM2[fuelType_int==6]<-0 #non-veg
#combine the two fire danger maps together
names(FDM) <- "FDM"
names(FDM2) <- "FDM2"
FDMs <- c(FDM, FDM2)
levelplot(FDMs, margin=F, main="Fire Danger (McArthur scale)", par.settings=infernoTheme)
Again, note the colour scale on the right. Where have the light coloured areas in the original map (on the left) gone?
We can also compare summaries of the values in each of our Fire Danger maps:
summary(FDMs)Warning: [summary] used a sample
FDM FDM2
Min. : 0.000 Min. : 0.00
1st Qu.: 1.769 1st Qu.: 6.65
Median : 1.769 Median : 6.65
Mean : 1.681 Mean : 45.79
3rd Qu.: 1.769 3rd Qu.: 6.65
Max. :20.956 Max. :249.62
We can also calculate RoS for the new climate values and compare:
#copy
RoS2 <- FDM2
#calculate
forestMulti<-max(FDM2[fuelType_int==1])*0.0012
grassRoS<-max(FDM2[fuelType_int==4]*0.13)
heathRoS<-max(FDM2[fuelType_int==5]*0.13)
#apply to appropriate mappes fuel types
RoS2[fuelType_int==1]<-biomass*forestMulti
RoS2[fuelType_int==2]<-biomass*forestMulti
RoS2[fuelType_int==3]<-grassRoS
RoS2[fuelType_int==4]<-heathRoS
RoS2[fuelType_int==5]<-grassRoS
RoS2[fuelType_int==6]<-0
#combine
names(RoS) <- "RoS"
names(RoS2) <- "RoS2"
RoSs <- c(RoS, RoS2)
#plot
levelplot(RoSs, margin=F, main="Rate of Spread (km per hour)", par.settings=infernoTheme)
As above, check you understand how the colour scale as affected the original map. Comparing summaries of the values in the raster maps might help.
summary(RoSs)Warning: [summary] used a sample
RoS RoS2
Min. :0.000 Min. : 0.000
1st Qu.:0.322 1st Qu.: 9.981
Median :2.724 Median : 9.981
Mean :1.805 Mean :11.501
3rd Qu.:2.724 3rd Qu.: 9.981
Max. :2.724 Max. :31.519
Try downloading this document and opening with R-Studio to work with the code interactively. Try changing climate variable values to see the effect on Fire Danger and Rate of Spread across the study area. An interactive version of this (online) will hopefully be ready soon.
Pasalodos-Tato, M., Ruiz-Peinado, R., Del Río, M., & Montero, G. (2015). Shrub biomass accumulation and growth rate models to quantify carbon stocks and fluxes for the Mediterranean region. European Journal of Forest Research, 134, 537-553. https://doi.org/10.1007/s10342-015-0870-6
Noble, I. R., Gill, A. M., & Bary, G. A. V. (1980). McArthur’s fire‐danger meters expressed as equations. Australian Journal of Ecology, 5(2), 201-203. https://doi.org/10.1111/j.1442-9993.1980.tb01243.x
Code and Text: James D.A. Millington, 2023 GNU GENERAL PUBLIC LICENSE
Land Cover Data owned by UK Centre for Ecology & Hydrology © Database Right/Copyright UKCEH