Hobart McArthur Forest Fire Danger Index

This script produces a graph of the current McArthur FFDI for Hobart, Tasmania, based on the drought factor (obtained from the Tasmanian Fire Service), and antecedent meteorological conditions.

First, we download the Tasmanian Fire Service's daily fire weather summary page, search for the line containing Hobart City, and grab the drought factor from it.

f = url("http://www.fire.tas.gov.au/mysite/Show?pageId=colWeatherfwxsdi")
ff = readLines(f)
for (te in 1:(length(ff) - 1)) {
    if (length(grep("Hobart City", ff[te], value = TRUE)) != 0 & ff[te] != "") {
        lin = te
    }
}

theline = ff[lin]
tlcomp = strsplit(theline, " ")
tlcomp = unlist(tlcomp)
tlcomp = tlcomp[tlcomp != ""]
df = as.numeric(tlcomp[5])

Now we access the Bureay of Meteorology's current observations table for Hobart, and extract Temperature, Relative Humidity and Windspeed into vectors.

g = url("http://www.bom.gov.au/fwo/IDT60901/IDT60901.94970.axf")
g = read.csv("http://www.bom.gov.au/fwo/IDT60901/IDT60901.94970.axf", skip = 19)
g = subset(g, !is.na(wmo))
g = subset(g, air_temp != -9999)

temp = g$air_temp
rh = g$rel_hum
ws = g$wind_spd_kmh

Here is the core Forest Fire Danger Index equation.

F <- 2 * exp(-0.45 + 0.987 * log(df + 0.001) - 0.0345 * rh + 0.0338 * temp + 
    0.0234 * ws)

We want to add a time field to our data. The way the Bureau of Meteorology stores date/time is a bit weird, so we're just going to do it roughly, and count backwards from the current system time in half-hour increments.

time_now = as.POSIXct(Sys.time())
subvec = seq(from = 0, by = 60 * 30, length.out = length(F))
time_vec = time_now - subvec
stt = time_vec[length(time_vec)]
ent = time_vec[1]

And here's the plot, with fire danger categories indicated in colours, and night time indicated by the grey vertical bars.

plot(F ~ time_vec, main = paste("Hobart McArthur FFDI:", round(F[1], 3), sep = ""), 
    lwd = 2, type = "l", xlab = "", ylab = "FFDI", ylim = c(0, 130))
rect(stt, 0, ent, 7.5, col = "#88FF0088")
rect(stt, 7.5, ent, 15, col = "#FFFF0088")
rect(stt, 15, ent, 25, col = "#FF880088")
rect(stt, 25, ent, 75, col = "#FF000088")
rect(stt, 75, ent, 100, col = "#CC004488")
rect(stt, 100, ent, 130, col = "#99008888")
for (cc in 1:length(time_vec)) {
    hour = as.numeric(format(time_vec[cc], "%H"))
    if (hour >= 18 | hour < 6) {
        rect(time_vec[cc], 0, time_vec[cc] + (30 * 60), 130, col = "#44444488", 
            border = "#444444FF")
    }
}