Time Series Prediction Model

A model for the Great Plains.

Get the tornado data. Set the domain and grid resolution.

if(!file.exists("torn")){ 
    download.file(url = "http://www.spc.noaa.gov/gis/svrgis/zipped/tornado.zip",
              destfile = "tornado.zip")
    unzip("tornado.zip")
  }
library("rgdal")

## Loading required package: sp

## rgdal: version: 1.0-4, (SVN revision 548)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.2/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.2/Resources/library/rgdal/proj
##  Linking to sp version: 1.1-1

TornL = readOGR(dsn = "torn", layer = "torn", 
                stringsAsFactors = FALSE)

## OGR data source with driver: ESRI Shapefile 
## Source: "torn", layer: "torn"
## with 58959 features
## It has 21 fields

library("raster")
xmn = -105; xmx = -95; ymn = 33; ymx = 48
r = raster(xmn = xmn, xmx = xmx,
           ymn = ymn, ymx = ymx,
           resolution = 5)
dim(r)

## [1] 3 2 1

nCells = ncell(r)
sp = as(r, 'SpatialPolygons')
spT = spTransform(sp, CRS(proj4string(TornL)))

Create a domain area map.

library("ggmap")

## Loading required package: ggplot2

Map = get_map(location = c(mean(c(xmn, xmx)), mean(c(ymn, ymx))), 
              source = "google",
              zoom = 4,
              color = "bw",
              maptype = "terrain")

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=40.5,-100&zoom=4&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

p1 = ggmap(Map, dev = "extent") +
  geom_segment(aes(x = xmn, xend = xmx, y = ymn, yend = ymn), 
             color = "red") +
  geom_segment(aes(x = xmn, xend = xmx, y = ymx, yend = ymx), 
             color = "red") +
  geom_segment(aes(x = xmn, xend = xmn, y = ymn, yend = ymx), 
             color = "red") +
  geom_segment(aes(x = xmx, xend = xmx, y = ymn, yend = ymx), 
             color = "red") +
  labs(x = expression(paste("Longitude (", degree, "W)")), 
       y = expression(paste("Latitude (", degree, "N)")))
library("maps")

## 
##  # ATTENTION: maps v3.0 has an updated 'world' map.        #
##  # Many country borders and names have changed since 1990. #
##  # Type '?world' or 'news(package="maps")'. See README_v3. #

library("maptools")

## Checking rgeos availability: TRUE

## 
## Attaching package: 'maptools'

## The following object is masked from 'package:sp':
## 
##     nowrapSpatialLines

library("grid")
p1 + theme(panel.grid.minor = element_line(colour = NA), 
           panel.grid.minor = element_line(colour = NA),
           panel.background = element_rect(fill = NA, colour = NA), 
           axis.text.x = element_blank(),
           axis.text.y = element_blank(), 
           axis.ticks.x = element_blank(),
             axis.ticks.y = element_blank(), 
           axis.title = element_blank(),
             rect = element_blank())

Figure: Prediction domain

Subset the tornado data. Overlay tornado paths on the raster. Each element of the list is a cell subset of the tornado track file. The data frame nTor.df contains the per-cell tornado count, area, and rate. The PathsAll.df data frame contains the tornado attributes listed by cell. Information is repeated for each cell in cases where the path intersects more than one cell. Area is square meter. The unadjusted climatological rate (climoU) has units of tornadoes per 100 x 100 km square per year.

library("rgeos")

## rgeos version: 0.3-11, (SVN revision 479)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Linking to sp version: 1.1-0 
##  Polygon checking: TRUE

beginYear = 1954
endYear = 2014
TornL2 = subset(TornL, yr >= beginYear & mag >= 0)
TnT = dim(TornL2)[1]
TornL2$Sid = 1:TnT
ct = over(spT, TornL2, returnList = TRUE)
nT = sapply(ct, function(x) nrow(x))
Nyears = endYear - beginYear + 1
area = rgeos::gArea(spT, byid = TRUE)
nTor.df = data.frame(ID = 1:nCells,
                     Nyears =  Nyears,
                     nT = nT, 
                     area = area, 
                     climoU =  nT/area * 10^6 * 100 * 100 / Nyears
                     )
PathsAll.df = plyr::ldply(ct, data.frame)
ID = rep(nTor.df$ID, nT)
PathsAll.df$ID = ID
TnTdom = length(unique(PathsAll.df$Sid))

Percentage of all tornadoes within the domain by EF rating.

TnTdom/TnT * 100

## [1] 32.2579

length(unique(PathsAll.df$Sid[PathsAll.df$mag >= 1]))/dim(TornL2[TornL2$mag >= 1, ])[1] * 100

## [1] 26.21719

length(unique(PathsAll.df$Sid[PathsAll.df$mag >= 2]))/dim(TornL2[TornL2$mag >= 2, ])[1] * 100

## [1] 27.05633

length(unique(PathsAll.df$Sid[PathsAll.df$mag >= 3]))/dim(TornL2[TornL2$mag >= 3, ])[1] * 100

## [1] 29.82394

length(unique(PathsAll.df$Sid[PathsAll.df$mag >= 4]))/dim(TornL2[TornL2$mag >= 4, ])[1] * 100

## [1] 31.8662

Get annual counts by domain. The data frame PathsAll.df contains the list of tornadoes by grid cell.

library("dplyr")

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:rgeos':
## 
##     intersect, setdiff, union

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

nTor.year = PathsAll.df %>%
  group_by(yr) %>%
  summarize(nT = length(unique(Sid)),
            nTd = length(unique(date)))
nTor.year = as.data.frame(nTor.year)
nTor.year$yr2 = nTor.year$yr
sum(nTor.year$nT)

## [1] 18657

library("ggplot2")
ggplot(nTor.year, aes(x = yr, y = nT)) + 
  geom_line() +
  geom_text(aes(label = nT), vjust = -.5, size = 4, color = "darkblue") +
  xlab("Year") + ylab("Annual Number of EF0+ Tornadoes") +
#  scale_y_continuous(limits = c(0, 2000)) +
#  scale_y_continuous(limits = c(0, 500)) +
  theme(axis.text.x  = element_text(size = 9), 
        legend.position = "none")

Figure: Time series tornado counts over the domain

moYrID.df = PathsAll.df %>%
  group_by(yr, mo, ID) %>%
  summarize(nT = n())
moYr.df = moYrID.df %>%
  group_by(yr, mo) %>%
  summarize(nT = sum(nT))

yr = rep(beginYear:endYear, each = 12)
mo = rep(1:12, Nyears)
moYrNull.df = data.frame(yr, mo, nT = 0)

df = left_join(moYrNull.df, moYr.df, by = c("yr", "mo"))
df[is.na(df)] = 0
df = df %>%
  select(yr, mo, nT = nT.y)

Get GAK and WCA monthly data.

library("tidyr")

## 
## Attaching package: 'tidyr'

## The following object is masked from 'package:raster':
## 
##     extract

GAK = read.table(file = "GAK2.txt")
names(GAK) = c("yr", month.abb)
df2 = gather(GAK, mo, SSTgak, Jan:Dec) %>%
  arrange(yr) %>%
  filter(yr >= beginYear, yr <= endYear)

WCA = read.table(file = "WCA2.txt")
names(WCA) = c("yr", month.abb)
df3 = gather(WCA, mo, SSTwca, Jan:Dec) %>%
  arrange(yr) %>%
  filter(yr >= beginYear, yr <= endYear)

df$SSTgak = df2$SSTgak
df$SSTwca = df3$SSTwca
df$SSTdif = df$SSTwca - df$SSTgak
write.csv(df, file = "ForEthan.csv")

Exploratory analysis

cor(df$SSTgak, df$SSTwca)
cor(df$nT, df$SSTdif)

df2 = df %>%
  filter(mo == 3)

library("INLA")
formula  = nT ~ SSTdif
model = inla(formula, family = "zeroinflatednbinomial0", data = df2)
summary(model)

INLA time series model. After “Drivers” example.

library("INLA")
formula  = nT ~ SSTwca + SSTgak + f(trend, model = "rw2") +
                          f(seasonal, model = "seasonal", season.length = 12, param = c(1, .1))
 
model = inla(formula, family = "zeroinflatednbinomial0", data = df)