Hotspot-climate correlations.

Charts of monthly MODIS hotspots versus rainfall and FFDI for three zones.

Set Zones

opts_knit$set(encoding = "UCS-2LE")
zones = list()
zones[[1]] = c(-9.9, -19.3, 129, 141)
zones[[2]] = c(-19.3, -31.7, 129, 141)
zones[[3]] = c(-31.7, -38.8, 129, 141)

Load hotspot data

hotspots = read.csv("E:\\NASA_NERP\\HotspotFFDI\\data\\ready\\hotspots.csv", 
    fileEncoding = "native.enc")
hotspots$acq_date = as.Date(hotspots$acq_date)
hotspots$year = format(hotspots$acq_date, "%Y")
hotspots$month = format(hotspots$acq_date, "%m")
hotspots$time_idx = as.numeric(hotspots$year) + (as.numeric(hotspots$month) - 
    1)/12

hot_list = list()
for (idx in seq_along(zones)) {
    hot_list[[idx]] = subset(hotspots, latitude < zones[[idx]][1] & latitude > 
        zones[[idx]][2] & longitude > zones[[idx]][3] & longitude < zones[[idx]][4] & 
        time_idx <= 2011.25)
}

st_date = min(hot_list[[1]]$time_idx)
en_date = max(hot_list[[1]]$time_idx)

Load station data

station_root = "G:\\SILO_Data\\complete_stations_min\\"
station_list = read.csv("E:\\NASA_NERP\\HotspotFFDI\\data\\ready\\silo_stations.csv")

stations = list()
for (idx in seq_along(zones)) {
    stations[[idx]] = subset(station_list, lat < zones[[idx]][1] & lat > zones[[idx]][2] & 
        long_ > zones[[idx]][3] & long_ < zones[[idx]][4])$no_
}

Average met data for months

met_list = list()
for (idx in seq_along(zones)) {
    stations_to_get = stations[[idx]]
    if (length(stations_to_get) > 20) {
        stations_to_get = sample(stations_to_get, 20)
    }
    totals = data.frame(station = NA, time_idx = NA, rain = NA, ffdi = NA)
    for (st in seq_along(stations_to_get)) {
        station_data = read.csv(paste(station_root, stations_to_get[st], ".csv", 
            sep = ""))
        station_data$c_date = as.Date(station_data$c_date)
        station_data$year = format(station_data$c_date, "%Y")
        station_data$month = format(station_data$c_date, "%m")
        station_data$time_idx = as.numeric(station_data$year) + (as.numeric(station_data$month) - 
            1)/12
        rain_data = tapply(station_data$rain, station_data$time_idx, sum)
        FFDI_data = tapply(station_data$IDMA, station_data$time_idx, mean)
        this.totals = data.frame(station = rep(st, length(rain_data)), time_idx = names(rain_data), 
            rain = as.numeric(rain_data), ffdi = as.numeric(FFDI_data))
        totals = rbind(totals, this.totals)
    }
    totals = subset(totals, !is.na(station))
    totals$time_idx = as.numeric(totals$time_idx)
    totals = aggregate(totals, by = list(totals$time_idx), mean)
    names(totals) = c("time_idx", "st", "dum", "rain", "ffdi")
    totals = subset(totals, time_idx >= st_date & time_idx <= en_date)
    met_list[[idx]] = totals
}

Rainfall plots

titles = c("Tropics", "Arid", "Temperate")
for (group in 1:3) {
    this.hotspots = hot_list[[group]]
    hot_vec = data.frame(time_idx = this.hotspots, frp = this.hotspots$frp)
    hot_vec = tapply(this.hotspots$frp, this.hotspots$time_idx, length)
    hot_vec = data.frame(count = as.numeric(hot_vec), time_idx = as.numeric(names(hot_vec)))
    op <- par(mar = c(5.1, 4.1, 4.1, 4.1))
    plot(hot_vec$count ~ hot_vec$time_idx, type = "l", xaxt = "n", lwd = 2, 
        xlab = "Date", ylab = "Hotspots", main = titles[group])
    axis(1, at = c(2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011))
    this.met = met_list[[group]]
    par(new = TRUE)
    plot(this.met$rain ~ this.met$time_idx, xaxt = "n", yaxt = "n", xlab = "", 
        ylab = "", type = "h", col = "blue")
    axis(4)
    mtext("Rain (mm)", side = 4, line = 2)
}

FFDI plots

titles = c("Tropics", "Arid", "Temperate")
for (group in 1:3) {
    this.hotspots = hot_list[[group]]
    hot_vec = data.frame(time_idx = this.hotspots, frp = this.hotspots$frp)
    hot_vec = tapply(this.hotspots$frp, this.hotspots$time_idx, length)
    hot_vec = data.frame(count = as.numeric(hot_vec), time_idx = as.numeric(names(hot_vec)))
    op <- par(mar = c(5.1, 4.1, 4.1, 4.1))
    plot(hot_vec$count ~ hot_vec$time_idx, type = "l", xaxt = "n", lwd = 2, 
        xlab = "Date", ylab = "Hotspots", main = titles[group])
    axis(1, at = c(2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011))
    this.met = met_list[[group]]
    par(new = TRUE)
    plot(this.met$ffdi ~ this.met$time_idx, xaxt = "n", yaxt = "n", xlab = "", 
        ylab = "", type = "h", col = "red")
    axis(4)
    mtext("FFDI", side = 4, line = 2)
}