Fire Index Shift Plot

In this plot, we represent FFDI values for a location on two axes, the y-axis representing year, and the x-axis representing day of year. FFDI value is represented by color. We might be able to see a shift in “season” through time.

We will perform this analysis for a single given cluster, in this case Perth (Zone 21).

Set up our libraries and data sources.

library(fields)  # Provides the quilt.plot function, nicest way to plot x-y data
library(reshape2)  # Provides the melt function
library(mgcv)  # Provides the gam function

station_root = "G:\\SILO_Data\\complete_stations\\"  # Where the met station data is stored
station_list = read.csv("E:\\NASA_NERP\\HotspotFFDI\\data\\ready\\silo_stations.csv")  # Station geographic locations

Slice up the data into a matrix for a given cluster.

zone = 21  # Let's do Perth
stations_to_get = subset(station_list, station_list$GRIDCODE == zone)$no_
if (length(stations_to_get) > 20) {
    # If more than 20 stations, pick 20 randomly to use
    stations_to_get = sample(stations_to_get, 20)
}
all_matrices = list()  # We're going to put all the matrices of weather variables in a list
for (idx in seq_along(stations_to_get)) {
    station = stations_to_get[idx]
    station_data = read.csv(paste(station_root, station, ".csv", sep = ""), 
        stringsAsFactors = F)
    station_data$c_date = as.Date(station_data$c_date)
    year_list = 1901:2010
    x_len = 365
    y_len = length(year_list)
    this_mat = matrix(0, nrow = y_len, ncol = x_len)
    for (the_year_idx in seq_along(year_list)) {
        the_year = year_list[the_year_idx]
        year_sub = subset(station_data, year == the_year)
        year_vals = as.numeric(year_sub$IDMA[1:365])
        this_mat[the_year_idx, ] = year_vals
    }
    rownames(this_mat) = year_list
    colnames(this_mat) = 1:365
    all_matrices[[idx]] = this_mat
}
matrix_sum = Reduce("+", all_matrices)  # Adds all the matrices in the list together, to make a single matrix
matrix_mean = matrix_sum/length(year_list)  # And now a mean - if only we could do this in one step
as_table = melt(matrix_mean)  # Now turn it into a data-frame

So here's a raw plot

quilt.plot(as_table, nx = length(year_list), ny = 365)  # Plot of raw data

And a smooth:

fit = gam(value ~ s(Var1, Var2, k = 5, bs = "tp"), data = as_table)  # GAM is the quickest, nicest way to smooth. Couldn't get Tps to work.
pred = predict(fit)
quilt.plot(x = as_table$Var1, y = as_table$Var2, z = pred, nx = length(year_list), 
    ny = 365)  # Plot of smoothed model

So fire seasons in Perth appear to be getting longer, and more extreme.