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 ordinary Krieging to generate a surface based on data from a few dozen stations. Interpolation as yet does not include elevation in the model, resulting in some fairly inrealistic representations in alpine areas.

First, we load the libraries we require. sp and maptools for the Shapefile and Raster processing, and gstat for the Krieg model.

# Load libraries
library(sp)
library(maptools)
library(raster)
library(gstat)

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.

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

# 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 then create a Krieg 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")
proj4string(ffdi_stats) = "+proj=longlat +datum=WGS84"

# Fit an ordinary Krieging model
v = variogram(ffdi ~ 1, ffdi_stats)
m = fit.variogram(v, vgm(1, "Sph", 200, 1))
ffdi.kr = gstat(NULL, "ffdi", ffdi ~ 1, ffdi_stats, model = m)

# Interpolate Krieg to the raster
op = interpolate(f_rast, ffdi.kr)

# 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, cex = 0.2)
plot(tas, add = T)