Tasmania McArthur Forest Fire Danger Index

This script produces an interpolated map of current Forest Fire Danger Index for Tasmania. It relies on data from the Tasmanian Fire Service, and the Bureau of Meteorology, and uses Generalized Additive Models to generate a surface based on data from a few dozen stations. Interpolation does include elevation in the model.

First, we load the libraries we require. sp, raster, maptools for the Shapefile and Raster processing, mgcv for the GAM model.

# Load libraries
library(sp)
library(maptools)
library(raster)
library(mgcv)
# library(fields)

The source data involves:
1. A table of stations, including their latitude, longitude, WMO code and look up name for the TFS drought factor data.
2. A ArcGIS shapefile of Tasmania.
3. A file made available from the TFS with daily drought factor values for Tasmanian locations.
4. Recent meteorological conditions tables for each station from the Bureau of Meteorology.
5. A raster of elevation in Tasmania, derived from the NASA SRTM.

# Load our table of stations to use
ffdi_stats = read.csv("ffdi_stats.csv")

# Load elevation raster
ele = raster("tas_dem.tif")

## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.9.0, released 2011/12/29 Path to GDAL shared
## files:
## /Library/Frameworks/R.framework/Versions/2.15/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480] Path to
## PROJ.4 shared files:
## /Library/Frameworks/R.framework/Versions/2.15/Resources/library/rgdal/proj

# Load Shapefile of Tasmania
tas = readShapePoly("tas_ol.shp")

# Get Drought Index file from TFS
f = url("http://www.fire.tas.gov.au/mysite/Show?pageId=colWeatherfwxsdi")
drought_file = readLines(f)

The following two functions retrieve the drought factor for a given station name, and temperature, humidity and wind speed for a given WMO code.

# Function to retrieve the drought factor for a given station
get_df = function(df.file, st.name) {
    for (te in 1:(length(df.file) - 1)) {
        if (length(grep(st.name, df.file[te], value = TRUE)) != 0 & df.file[te] != 
            "") {
            lin = te
        }
    }

    theline = df.file[lin]
    tlcomp = strsplit(theline, " ")
    tlcomp = unlist(tlcomp)
    tlcomp = tlcomp[tlcomp != ""]
    df = as.numeric(tlcomp[5])
    if (is.na(df)) {
        df = as.numeric(tlcomp[4])
    }
    df
}

# Function to retrieve latest RH, TEMP, WS from BOM
get_met = function(wmo_code) {
    to_get = paste0("http://www.bom.gov.au/fwo/IDT60801/IDT60801.", wmo_code, 
        ".axf")
    g = read.csv(to_get, skip = 19)
    g = subset(g, !is.na(wmo))
    g = subset(g, air_temp != -9999)

    temp = g$air_temp[1]
    rh = g$rel_hum[1]
    ws = g$wind_spd_kmh[1]
    c(temp, rh, ws)
}

Here is the meat of the code - we go through each station in the table, retrieve the drought factor and other variables, calculate the FFDI and store it back in the data table. We join the elevation data, then we then create a GLM model and interpolate it to a raster for display.

# Make a new column for the FFDI
ffdi_stats$ffdi = 0

# Calculate FFDI for each station, add to table
for (idx in seq_along(ffdi_stats$Station)) {
    wmo = ffdi_stats$File[idx]
    df_line = ffdi_stats$Line[idx]
    df = get_df(drought_file, df_line)
    met_vars = get_met(wmo)
    rh = met_vars[2]
    temp = met_vars[1]
    ws = met_vars[3]
    F <- 2 * exp(-0.45 + 0.987 * log(df + 0.001) - 0.0345 * rh + 0.0338 * temp + 
        0.0234 * ws)
    ffdi_stats$ffdi[idx] = F
}

# Make a new raster object
# f_rast=raster(nrows=100,ncols=100,xmn=144.6,xmx=148.5,ymn=-43.8,ymx=-40.6,crs='+proj=longlat
# +datum=WGS84')

# Convert station table to SpatialPoints object
coordinates(ffdi_stats) = c("Lon", "Lat")

# Extract elevation raster to points. Some islands don't have elevations,
# so make them 5m.
e_vec = extract(ele, ffdi_stats)
e_vec[is.na(e_vec)] = 5
ffdi_stats$ele = e_vec

# For the GAM model, we need coordinates in the frame.
proj4string(ffdi_stats) = "+proj=longlat +datum=WGS84"
ffdi_stats$x = coordinates(ffdi_stats)[, 1]
ffdi_stats$y = coordinates(ffdi_stats)[, 2]

# Run the GAM model, predict back to a raster.
gam.mod = gam(ffdi ~ s(x, y, k = 5) + s(ele, k = 5), data = ffdi_stats)
xy <- data.frame(xyFromCell(ele, 1:ncell(ele)))
xy$ele <- values(ele)
gam.p = predict(gam.mod, newdata = xy)
op = ele
values(op) = as.vector(gam.p)

# Mask the raster to Tasmania
op_m = mask(op, tas)

And finally, produce our plot.

# Plot raster predictions, stations, and state outline
plot(op_m, main = "FFDI", col = rev(rainbow(20)))
points(ffdi_stats)
plot(tas, add = T)