GIS5122 Applied Spatial Statistics
Tachanat Bhatrasataponkul (tachanat@coaps.fsu.edu)
Class Project: Tornadoes occur more often over smooth terrain.
Objective(s): To investigate the relationship between tornado occurrence and terrain roughness across the Great Plains during 1954-2014 (61 years).
Different Roughness Indices:
roughness (m)
slope (degrees, radians, tangent)
aspect
TRI: Terrain Ruggedness Index
TPI: Topographic Position Index
According to Wilson et al. (2007), Terrain Ruggedness Index (TRI) is the mean of the absolute differences between the value of a cell and the value of its 8 surrounding cells whereas Topographic Position Index (TPI) is the difference between the value of a cell and the mean value of its 8 surrounding cells. Roughness is the difference between the maximum and the minimum value of a cell and its 8 surrounding cells. Slope and aspect are computed according to Fleming and Hoffer (1979), Horn (1981) and Ritter (1987). The Horn algorithm may be best for rough surfaces, and the Fleming and Hoffer algorithm may be better for smoother surfaces (Jones, 1998; Burrough and McDonnell, 1998).
References:
Burrough, P., and R.A. McDonnell, 1998. Principles of Geographical Information Systems. Oxford University Press.
Fleming, M.D. and Hoffer, R.M., 1979. Machine processing of landsat MSS data and DMA topographic data for forest cover type mapping. LARS Technical Report 062879. Laboratory for Applications of Remote Sensing, Purdue University, West Lafayette, Indiana.
Greenlee, D.D., 1987. Raster and vector processing for scanned linework. Photogrammetric Engineering and Remote Sensing 53:1383-1387
Horn, B.K.P., 1981. Hill shading and the reflectance map. Proceedings of the IEEE 69:14-47
Jones, K.H., 1998. A comparison of algorithms used to compute hill slope as a property of the DEM. Computers & Geosciences 24: 315-323
Ritter, P., 1987. A vector-based slope and aspect generation algorithm. Photogrammetric Engineering and Remote Sensing 53: 1109-1111
Wilson, M.F.J., O’Connell, B., Brown, C., Guinan, J.C., Grehan, A.J., 2007. Multiscale terrain analysis of multibeam bathymetry data for habitat mapping on the continental slope. Marine Geodesy 30: 3-35.
Tornado counts by raster cell
#download.file(url = "http://www.spc.noaa.gov/gis/svrgis/zipped/tornado.zip", destfile = "tornado.zip")
download.file(url = "http://coaps.fsu.edu/~tachanat/data/tornado1950-2014.zip", destfile = "tornado.zip")
unzip("tornado.zip")
library("rgdal")## Loading required package: sp
## rgdal: version: 0.9-2, (SVN revision 526)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.1, released 2014/09/24
## Path to GDAL shared files: C:/Users/Tachanat/Documents/R/win-library/3.1/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: C:/Users/Tachanat/Documents/R/win-library/3.1/rgdal/projTornL = readOGR(dsn = "torn", layer = "torn", stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "torn", layer: "torn"
## with 58882 features
## It has 21 fieldssummary(TornL)## Object of class SpatialLinesDataFrame
## Coordinates:
## min max
## x -2345551 3166079
## y -1819794 1388488
## Is projected: TRUE
## proj4string :
## [+proj=lcc +lat_1=33 +lat_2=45 +lat_0=39 +lon_0=-96 +x_0=0 +y_0=0
## +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0]
## Data attributes:
## om yr mo dy
## Min. : 1.0 Min. :1950 Min. : 1.000 Min. : 1.00
## 1st Qu.: 243.0 1st Qu.:1973 1st Qu.: 4.000 1st Qu.: 8.00
## Median : 499.0 Median :1990 Median : 6.000 Median :16.00
## Mean : 30305.0 Mean :1987 Mean : 5.964 Mean :15.86
## 3rd Qu.: 820.8 3rd Qu.:2003 3rd Qu.: 7.000 3rd Qu.:24.00
## Max. :553549.0 Max. :2014 Max. :12.000 Max. :31.00
## date time tz st
## Length:58882 Length:58882 Min. :0.000 Length:58882
## Class :character Class :character 1st Qu.:3.000 Class :character
## Mode :character Mode :character Median :3.000 Mode :character
## Mean :3.001
## 3rd Qu.:3.000
## Max. :9.000
## stf stn mag inj
## Min. : 1.0 Min. : 0.00 Min. :-9.0000 Min. : 0.000
## 1st Qu.:18.0 1st Qu.: 5.00 1st Qu.: 0.0000 1st Qu.: 0.000
## Median :29.0 Median : 16.00 Median : 1.0000 Median : 0.000
## Mean :29.4 Mean : 27.03 Mean : 0.7933 Mean : 1.579
## 3rd Qu.:45.0 3rd Qu.: 36.00 3rd Qu.: 1.0000 3rd Qu.: 0.000
## Max. :72.0 Max. :232.00 Max. : 5.0000 Max. :1740.000
## fat loss closs
## Min. : 0.00000 Min. : 0.000 Min. : 0.000000
## 1st Qu.: 0.00000 1st Qu.: 0.000 1st Qu.: 0.000000
## Median : 0.00000 Median : 0.100 Median : 0.000000
## Mean : 0.09835 Mean : 2.197 Mean : 0.002117
## 3rd Qu.: 0.00000 3rd Qu.: 4.000 3rd Qu.: 0.000000
## Max. :158.00000 Max. :2800.100 Max. :23.520000
## slat slon elat elon
## Min. :18.33 Min. :-124.47 Min. : 0.00 Min. :-124.46
## 1st Qu.:33.26 1st Qu.: -98.58 1st Qu.: 0.00 1st Qu.: -93.93
## Median :37.10 Median : -93.88 Median :30.85 Median : -81.14
## Mean :37.18 Mean : -92.89 Mean :20.64 Mean : -50.97
## 3rd Qu.:40.98 3rd Qu.: -86.82 3rd Qu.:38.05 3rd Qu.: 0.00
## Max. :49.33 Max. : -67.19 Max. :49.33 Max. : 0.00
## len wid
## Min. : 0.00 Min. : 0.00
## 1st Qu.: 0.10 1st Qu.: 13.00
## Median : 0.53 Median : 40.00
## Mean : 3.49 Mean : 97.58
## 3rd Qu.: 3.00 3rd Qu.: 100.00
## Max. :234.70 Max. :4576.00library("raster")
r = raster(xmn = -102, xmx = -95,
ymn = 36, ymx = 42,
resolution = .25)
sp = as(r, 'SpatialPolygons')
spT = spTransform(sp, CRS(proj4string(TornL)))
library("GISTools")## Loading required package: maptools
## Checking rgeos availability: TRUE
## Loading required package: RColorBrewer
## Loading required package: MASS
##
## Attaching package: 'MASS'
##
## The following objects are masked from 'package:raster':
##
## area, select
##
## Loading required package: rgeos
## rgeos version: 0.3-8, (SVN revision 460)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Polygon checking: TRUElibrary("PBSmapping")##
## -----------------------------------------------------------
## PBS Mapping 2.68.68 -- Copyright (C) 2003-2015 Fisheries and Oceans Canada
##
## PBS Mapping comes with ABSOLUTELY NO WARRANTY;
## for details see the file COPYING.
## This is free software, and you are welcome to redistribute
## it under certain conditions, as outlined in the above file.
##
## A complete user guide 'PBSmapping-UG.pdf' is located at
## C:/Users/Tachanat/Documents/R/win-library/3.1/PBSmapping/doc/PBSmapping-UG.pdf
##
## Packaged on 2015-01-14
## Pacific Biological Station, Nanaimo
##
## All available PBS packages can be found at
## http://code.google.com/p/pbs-software/
##
## To see demos, type '.PBSfigs()'.
## -----------------------------------------------------------data(tornados)
par(mar=c(0,0,0,0))
plot(us_states)
points(TornL$elon, TornL$elat, col="#FB6A4A4C", pch=1, cex=0.4)
Figure 1 Tornado occurrence over the U.S. during 1950-2014
Spatial subset.
library("ggmap")## Loading required package: ggplot2Map = get_map(location = c(-98.5, 39),
source = "google",
zoom = 6,
color = "bw",
maptype = "terrain")## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=39,-98.5&zoom=6&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=falsep1 = ggmap(Map, dev = "extent") +
geom_segment(aes(x = -102, xend = -95, y = 36, yend = 36),
color = "red") +
geom_segment(aes(x = -102, xend = -95, y = 42, yend = 42),
color = "red") +
geom_segment(aes(x = -102, xend = -102, y = 36, yend = 42),
color = "red") +
geom_segment(aes(x = -95, xend = -95, y = 36, yend = 42),
color = "red") +
labs(x = expression(paste("Longitude (", degree, "W)")),
y = expression(paste("Latitude (", degree, "N)")))
library("maps")
library("maptools")
library("grid")
#source('ScaleBarNorth.R')
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 2 Study region.
The CRS of the tornado data is Lambert conic conformal (lcc) with units of meters. We create a raster within a bounding box and at a resolution of .25 degree. The domain is centered on Kansas but includes parts of Nebraska, Iowa, Oklahoma, Missouri, and Texas. Change the raster to spatial polygons. This is used to perform the overlay of the tornado tracks. Convert the polygon grid CRS to the projection of the tornado tracks.
There are 24 cells north to south (rows) and 28 cells east to west (cols) for a total of 672. The CRS is geographic. Overlay the tornado tracks onto the spatial polygons with the over() function from the sp package. It works with various combinations of spatial objects but the rgeos package is needed for it to work with spatial polygons and spatial lines.
Remove duplicate tracks. Overlay tracks onto cells and return a vector indicating either NA (is not inside extent) or the cell number. Subset the tornadoes by the raster extent. Compute the track lengths and store in data frame with start and end lat/lon. Subset by duplicates of this data frame and paths exceeding 100 m.
Temporal subset.
library("rgeos")
tc = over(TornL, spT)
TornL2 = subset(TornL, !is.na(tc))
TornL2 = subset(TornL2, as.numeric(yr) >= 1954)
dim(TornL2@data)## [1] 6990 21PathD = SpatialLinesLengths(TornL2)
df = data.frame(PathD,
SLAT = as.numeric(TornL2$slat),
SLON = as.numeric(TornL2$slon),
ELAT = as.numeric(TornL2$elat),
ELON = as.numeric(TornL2$elon),
DATE = TornL2$date)
dup = duplicated(df)
sum(dup)## [1] 111sum(dup)/dim(TornL2@data)[1] * 100## [1] 1.587983TornL3 = subset(TornL2, !dup)Add a buffer to the tracks based on the width specified in the attribute table to create a spatial polygons data frame. Overlay the path polygons on the cell polygons. Add minimum widths to tornadoes without width.
Width = as.numeric(TornL3$wid) * .9144
sum(Width[Width == 0])## [1] 0#Width[Width == 0] = min(Width[Width != 0])
TornP3 = gBuffer(TornL3, byid = TRUE, width = Width/2)
ct = over(spT, TornP3, returnList = TRUE)
nT = sapply(ct, function(x) length(x))
sumL = sapply(ct, function(x) sum(TornP3$Length[x]))
nt = r
values(nt) = nT
cellStats(nt, stat = "mean")## [1] 14.51042cellStats(nt, stat = "sd")## [1] 7.241823cellStats(nt, stat = "sd")^2/cellStats(nt, stat = "mean")## [1] 3.614231The result is a list. Each element of the list is a vector of row numbers from the attribute table corresponding to paths occurring in the cell. There are 672 cells. The order of the list matches the order of the raster (upper left to lower right in lexigraphical order). From the list we get the length of each vector. This is done with the sapply() (simple apply) where the first argument is the list object and the second argument is a function. These counts are then placed onto the original raster with the values() function. The GIS overlay operation uses projected coordinates. The analysis/display uses geographic coordinates.
Compare using tornadoes as track lines.
ct2 = over(spT, TornL3, returnList = TRUE)
nT2 = sapply(ct2, function(x) nrow(x))
(sum(nT) - sum(nT2))/sum(nT) * 100## [1] 3.066352nt2 = r
values(nt2) = nT2
cellStats(nt2, stat = "mean")## [1] 14.06548cellStats(nt2, stat = "sd")## [1] 7.015112cellStats(nt2, stat = "sd")^2/cellStats(nt2, stat = "mean")## [1] 3.498765There are 3.1% more tornado counts using tornado paths than using tornado tracks.
Create map
Get county & state boundaries. Get the tornado tracks. Convert to the geographic coordinates of the raster.
library("mapproj")
library("maptools")
ext = as.vector(extent(r))
bndryC = map("county", fill = TRUE,
xlim = ext[1:2],
ylim = ext[3:4],
plot = FALSE)
IDs = sapply(strsplit(bndryC$names, ":"), function(x) x[1])
bndryCP = map2SpatialPolygons(bndryC, IDs = IDs,
proj4string = CRS(projection(r)))
bndryS = map("state", fill = TRUE,
xlim = ext[1:2],
ylim = ext[3:4],
plot = FALSE)
IDs = sapply(strsplit(bndryS$names, ":"), function(x) x[1])
bndrySP = map2SpatialPolygons(bndryS, IDs = IDs,
proj4string = CRS(projection(r)))
TornP3T = spTransform(TornP3, CRS(projection(r)))Create the map.
library("rasterVis")## Loading required package: lattice
## Loading required package: latticeExtra
##
## Attaching package: 'latticeExtra'
##
## The following object is masked from 'package:ggplot2':
##
## layerlibrary("wesanderson")
range(values(nt))## [1] 1 48rng = seq(0, 50, 5)
par(mar=c(0,0,0,0))
cr = wes_palette(name = "Zissou", n = 10, type = "continuous")
p2 = levelplot(nt, margin = TRUE,
sub = "Tornado Counts ",
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(space = 'bottom'),
par.settings = list(fontsize = list(text = 15)))
p2 +
latticeExtra::layer(sp.polygons(bndryCP,
col = gray(.85), lwd = 1)) +
latticeExtra::layer(sp.polygons(bndrySP,
col = gray(.15), lwd = 1)) +
latticeExtra::layer(sp.polygons(TornP3T, fill = gray(.4),
col = gray(.5), alpha = .3))
Figure 3 Tornado counts. Paths are shown in gray and the number of tornadoes intersecting each cell is shown with a color ramp. Row and column counts are shown with plots in the left and top margins respectively.
Get elevation raster from a DEM
Digital elevation model data are available from http://www.viewfinderpanoramas.org/DEM/TIF15/ Get the elevation raster and crop it to the extent of the tornado raster r. Compute the elevation roughness with the terrain() function. Here roughness is the difference between the maximum and the minimum elevation in the cell and the eight surrounding cells. Match the resolution and origin of the elevation rasters with those of the tornado raster by degrading the resolution of elevation roughness. Compute the correlation between the number of tornadoes and surface roughness & elevation.
#download.file(url = "http://www.viewfinderpanoramas.org/DEM/TIF15/15-H.ZIP", destfile = "15-H.ZIP", mode = "wb")
#unzip("15-H.ZIP")
download.file(url = "http://coaps.fsu.edu/~tachanat/data/15-H.tif.zip", destfile = "15-H.tif.zip", mode = "wb")
unzip("15-H.tif.zip")
Elev = raster("15-H.tif")
Elev = crop(Elev, nt)
TR = terrain(Elev, opt = 'roughness')
TS = terrain(Elev, opt = 'slope', unit='degrees')
TS_D = terrain(Elev, opt = 'slope', unit='degrees')
TS_R = terrain(Elev, opt = 'slope', unit='radians')
TS_T = terrain(Elev, opt = 'slope', unit='tangent')
TA = terrain(Elev, opt = 'aspect')
TRI = terrain(Elev, opt = 'TRI')
TPI = terrain(Elev, opt = 'TPI')
el = resample(aggregate(Elev, fact = c(nrow(r), ncol(r)), fun = mean), r)
tr = resample(aggregate(TR, fact = c(nrow(r), ncol(r)), fun = mean), r)
ts = resample(aggregate(TS, fact = c(nrow(r), ncol(r)), fun = mean), r)
ta = resample(aggregate(TA, fact = c(nrow(r), ncol(r)), fun = mean), r)
ts_d = resample(aggregate(TS_D, fact = c(nrow(r), ncol(r)), fun = mean), r)
ts_r = resample(aggregate(TS_R, fact = c(nrow(r), ncol(r)), fun = mean), r)
ts_t = resample(aggregate(TS_T, fact = c(nrow(r), ncol(r)), fun = mean), r)
tri = resample(aggregate(TRI, fact = c(nrow(r), ncol(r)), fun = mean), r)
tpi = resample(aggregate(TPI, fact = c(nrow(r), ncol(r)), fun = mean), r)
cellStats(el, stat = "mean"); cellStats(el, stat = "sd")## [1] 592.475## [1] 249.5601cellStats(tr, stat = "mean"); cellStats(tr, stat = "sd")## [1] 14.20804## [1] 6.503738cellStats(ts, stat = "mean"); cellStats(ts, stat = "sd")## [1] 0.6588277## [1] 0.2894508cellStats(ta, stat = "mean"); cellStats(ta, stat = "sd")## [1] 2.768552## [1] 0.4088319cellStats(tri, stat = "mean"); cellStats(tri, stat = "sd")## [1] 4.718679## [1] 2.251977cellStats(tpi, stat = "mean"); cellStats(tpi, stat = "sd")## [1] 0.0007131052## [1] 0.04314302cor(values(el), values(nt)); cor(values(tr), values(nt)); cor(values(ts), values(nt)); cor(values(ta), values(nt)); cor(values(tri), values(nt)); cor(values(tpi), values(nt))## [1] -0.2709387## [1] -0.3647149## [1] -0.3687797## [1] -0.06500166## [1] -0.361772## [1] 0.0001265339Get a population raster
Gridded Population of the World (GPW), v3. http://sedac.ciesin.columbia.edu/data/set/gpw-v3-population-density/data-download Values are persons per square km. Download as grid and read as raster. Crop to extent of tornado raster and match the resolution of the tornado grid.
download.file(url = "http://myweb.fsu.edu/jelsner/data/usadens.zip", destfile = "usadens.zip")
unzip("usadens.zip")
Pop = raster("usadens/usads00g/w001001.adf")
Pop = crop(Pop, r)
pop = resample(aggregate(Pop, fact = c(nrow(r), ncol(r)), fun = mean), r)
cellStats(pop, stat = "mean"); cellStats(pop, stat = "sd")## [1] 12.58063## [1] 17.42387cor(values(pop), values(nt))## [1] 0.223462cor(values(pop), values(el))## [1] -0.5615524Exploratory analysis
Scatter plots
library("ggplot2")
df = data.frame(nT = values(nt),
el = values(el),
tr = values(tr),
ts = values(ts),
ta = values(ta),
ts_d = values(ts_d),
ts_r = values(ts_r),
ts_t = values(ts_t),
tri = values(tri),
tpi = values(tpi),
pop = values(pop))
f3a = ggplot(df, aes(x = log2(pop), y = log(nT))) +
geom_point() +
geom_smooth(method = lm) +
ylab("Number of Tornadoes (log)") +
xlab(expression(paste("2000 Population Density (people per ", km^2, ")"))) +
scale_x_continuous(breaks = c(1, 2, 3, 4, 5),
labels = c(1, 2, 3, 4, 5))
f3b = ggplot(df, aes(x = el, y = log(nT))) +
geom_point() +
geom_smooth(method = lm) +
ylab("Number of Tornadoes (log)") +
xlab("Elevation (m)")
f3c = ggplot(df, aes(x = tr, y = log(nT))) +
geom_point() +
geom_smooth(method = loess) +
ylab("Number of Tornadoes (log)") +
xlab("Terrain Roughness (m)")
f3d = ggplot(df, aes(x = ts_d, y = log(nT))) +
geom_point() +
geom_smooth(method = loess) +
ylab("Number of Tornadoes (log)") +
xlab("Terrain Slope (degrees)")
f3e = ggplot(df, aes(x = ts_r, y = log(nT))) +
geom_point() +
geom_smooth(method = loess) +
ylab("Number of Tornadoes (log)") +
xlab("Terrain Slope (radians)")
f3f = ggplot(df, aes(x = ts_t, y = log(nT))) +
geom_point() +
geom_smooth(method = loess) +
ylab("Number of Tornadoes (log)") +
xlab("Terrain Slope (tangent)")
f3g = ggplot(df, aes(x = ta, y = log(nT))) +
geom_point() +
geom_smooth(method = loess) +
ylab("Number of Tornadoes (log)") +
xlab("Terrain Aspect")
f3h = ggplot(df, aes(x = tri, y = log(nT))) +
geom_point() +
geom_smooth(method = loess) +
ylab("Number of Tornadoes (log)") +
xlab("Terrain Ruggedness Index")
f3i = ggplot(df, aes(x = tpi, y = log(nT))) +
geom_point() +
geom_smooth(method = loess) +
ylab("Number of Tornadoes (log)") +
xlab("Topographic Position Index")
download.file(url = "http://coaps.fsu.edu/~tachanat/codes/multiplot.txt", destfile = "multiplot.txt")
source("http://coaps.fsu.edu/~tachanat/codes/multiplot.txt")
source("multiplot.txt")
mat1 = matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9), nrow = 3, byrow = TRUE)
f3a = f3a + ggtitle("(a)") + theme(plot.title = element_text(hjust = 0))
f3b = f3b + ggtitle("(b)") + theme(plot.title = element_text(hjust = 0))
f3c = f3c + ggtitle("(c)") + theme(plot.title = element_text(hjust = 0))
f3d = f3d + ggtitle("(d)") + theme(plot.title = element_text(hjust = 0))
f3e = f3e + ggtitle("(e)") + theme(plot.title = element_text(hjust = 0))
f3f = f3f + ggtitle("(f)") + theme(plot.title = element_text(hjust = 0))
f3g = f3g + ggtitle("(g)") + theme(plot.title = element_text(hjust = 0))
f3h = f3h + ggtitle("(h)") + theme(plot.title = element_text(hjust = 0))
f3i = f3i + ggtitle("(i)") + theme(plot.title = element_text(hjust = 0))
multiplot(f3a, f3b, f3c, f3d, f3e, f3f, f3g, f3h, f3i, layout = mat1)
Figure 4 Tornadoes versus population, elevation, terrain roughness and terrain slope. The number of tornadoes in each grid cell is given on a log scale. The population density is on a log (base two) scale.
Histogram of the number of tornadoes
p3 = ggplot(df, aes(x = nT)) +
geom_histogram(binwidth = 2, color = "white") +
ylab("Number of Cells") +
xlab("Number of Tornadoes")
p3 + theme_bw()
Figure 5 Histogram of the number of tornadoes by grid cell. The most tornadoes in any cell is 48 and the fewest is one. There are 81 cells with counts of ten or eleven. Cell counts are more dispersed than a Poisson distribution.
Convert raster layers to polygons
Convert raster layers to polygons. Change the attribute name to nT in the spatial polygons data frame and add the elevation roughness as another column.
spdf = as(nt, "SpatialPolygonsDataFrame")
names(spdf) = "nT"
spdf$el = values(el)
spdf$tr = values(tr)
spdf$ts = values(ts)
spdf$tri = values(tri)
spdf$pop = values(pop)
spdfT = spTransform(spdf, CRS(proj4string(TornP3)))
spdf$area = gArea(spdfT, byid = TRUE)
spdf$l2pop = log2(spdf$pop)
spdf$ID = 1:ncell(nt)
cor(spdf$el, spdf$pop)## [1] -0.5615524cor(spdf$el, spdf$tr); cor(spdf$el, spdf$ts); cor(spdf$el, spdf$tri)## [1] -0.07755983## [1] -0.0778609## [1] -0.0728845cor(spdf$pop, spdf$tr); cor(spdf$pop, spdf$ts); cor(spdf$pop, spdf$tri)## [1] 0.1249989## [1] 0.1179515## [1] 0.1269868Spatial model
Some controls for INLA.
source("http://www.math.ntnu.no/inla/givemeINLA.R")##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:base':
##
## crossprod, tcrossprod
##
##
## Attaching package: 'INLA'
##
## The following object is masked _by_ '.GlobalEnv':
##
## inla.update## Warning: package 'Matrix' in library 'C:/Program Files/R/R-3.1.3/library'
## will not be updated##
##
##
## You can later upgrade INLA using: inla.upgrade()library("INLA")Spatial neighborhood definition as an inla graph
library("spdep")
nb = poly2nb(spdf)
nb2INLA("g", nb)
g = inla.read.graph("g")Model for the smoothed relative rate.
formula0 = nT ~ f(ID, model = "besag", graph = g)
model0 = inla(formula0, family = "nbinomial", E = area/10^6,
data = spdf@data)
summary(model0)##
## Call:
## c("inla(formula = formula0, family = \"nbinomial\", data = spdf@data, ", " E = area/10^6)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.5114 2.3090 0.5094 3.3298
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -3.7794 0.0147 -3.808 -3.7795 -3.7504 -3.7796 0
##
## Random effects:
## Name Model
## ID Besags ICAR model
##
## Model hyperparameters:
## mean sd
## size for the nbinomial observations (overdispersion) 17.664 2.6635
## Precision for ID 2.677 0.4032
## 0.025quant 0.5quant
## size for the nbinomial observations (overdispersion) 12.977 17.480
## Precision for ID 1.979 2.644
## 0.975quant mode
## size for the nbinomial observations (overdispersion) 23.427 17.127
## Precision for ID 3.556 2.579
##
## Expected number of effective parameters(std dev): 212.57(23.84)
## Number of equivalent replicates : 3.161
##
## Marginal Likelihood: -2774.53rSR = r
values(rSR) = (exp(model0$summary.random$ID$mean) - 1) * 100Plot the smoothed rate.
library("RColorBrewer")
range(values(rSR))## [1] -63.85859 139.13998rng = seq(-100, 150, 50)
rngL = paste(rng, '%', sep = "")
cr = rev(brewer.pal(8, "RdBu"))
cr = cr[-(1:2)]
p4 = levelplot(rSR, margin = TRUE,
sub = "Tornado Reports (Above/Below Regional Average)",
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at = rng, labels = rngL, col = cr),
par.settings = list(fontsize = list(text = 13)))
p4 +
latticeExtra::layer(sp.polygons(bndryCP, col = gray(.85), lwd = 1)) +
latticeExtra::layer(sp.polygons(bndrySP, col = gray(.15), lwd = 1))
Figure 6 Tornado Reports (Above/Below Regional Average)
Add covariates to the model
formula1tr = nT ~ f(ID, model = "besag", graph = g) + el + l2pop + tr
model1tr = inla(formula = formula1tr, family = "nbinomial", E = area/10^6, data = spdf@data)
summary(model1tr)##
## Call:
## c("inla(formula = formula1tr, family = \"nbinomial\", data = spdf@data, ", " E = area/10^6)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.3277 1.9868 0.0747 2.3892
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -3.5605 0.2608 -4.0782 -3.5587 -3.0528 -3.5551 0
## el -0.0001 0.0003 -0.0008 -0.0001 0.0005 -0.0001 0
## l2pop 0.0849 0.0338 0.0187 0.0848 0.1517 0.0846 0
## tr -0.0267 0.0035 -0.0335 -0.0267 -0.0198 -0.0267 0
##
## Random effects:
## Name Model
## ID Besags ICAR model
##
## Model hyperparameters:
## mean sd
## size for the nbinomial observations (overdispersion) 16.494 2.4218
## Precision for ID 4.336 0.9166
## 0.025quant 0.5quant
## size for the nbinomial observations (overdispersion) 12.172 16.354
## Precision for ID 2.879 4.216
## 0.975quant mode
## size for the nbinomial observations (overdispersion) 21.653 16.104
## Precision for ID 6.455 3.974
##
## Expected number of effective parameters(std dev): 158.06(25.27)
## Number of equivalent replicates : 4.252
##
## Marginal Likelihood: -2770.07formula1ts = nT ~ f(ID, model = "besag", graph = g) + el + l2pop + ts
model1ts = inla(formula = formula1ts, family = "nbinomial", E = area/10^6, data = spdf@data)
summary(model1ts)##
## Call:
## c("inla(formula = formula1ts, family = \"nbinomial\", data = spdf@data, ", " E = area/10^6)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.2154 1.6340 0.0625 1.9119
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -3.5293 0.2613 -4.0481 -3.5275 -3.0208 -3.5239 0
## el -0.0002 0.0003 -0.0008 -0.0002 0.0005 -0.0002 0
## l2pop 0.0836 0.0338 0.0175 0.0835 0.1503 0.0833 0
## ts -0.5878 0.0777 -0.7399 -0.5881 -0.4347 -0.5885 0
##
## Random effects:
## Name Model
## ID Besags ICAR model
##
## Model hyperparameters:
## mean sd
## size for the nbinomial observations (overdispersion) 16.346 2.4009
## Precision for ID 4.389 0.9388
## 0.025quant 0.5quant
## size for the nbinomial observations (overdispersion) 12.070 16.20
## Precision for ID 2.889 4.27
## 0.975quant mode
## size for the nbinomial observations (overdispersion) 21.479 15.940
## Precision for ID 6.548 4.029
##
## Expected number of effective parameters(std dev): 156.62(25.38)
## Number of equivalent replicates : 4.291
##
## Marginal Likelihood: -2767.53formula1tri = nT ~ f(ID, model = "besag", graph = g) + el + l2pop + tri
model1tri = inla(formula = formula1tri, family = "nbinomial", E = area/10^6, data = spdf@data)
summary(model1tri)##
## Call:
## c("inla(formula = formula1tri, family = \"nbinomial\", data = spdf@data, ", " E = area/10^6)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.2167 1.6387 0.0660 1.9214
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -3.5970 0.2601 -4.1135 -3.5952 -3.0908 -3.5916 0
## el -0.0001 0.0003 -0.0008 -0.0001 0.0006 -0.0001 0
## l2pop 0.0873 0.0339 0.0210 0.0872 0.1541 0.0870 0
## tri -0.0780 0.0102 -0.0981 -0.0780 -0.0579 -0.0781 0
##
## Random effects:
## Name Model
## ID Besags ICAR model
##
## Model hyperparameters:
## mean sd
## size for the nbinomial observations (overdispersion) 16.53 2.4288
## Precision for ID 4.31 0.9068
## 0.025quant 0.5quant
## size for the nbinomial observations (overdispersion) 12.197 16.392
## Precision for ID 2.868 4.193
## 0.975quant mode
## size for the nbinomial observations (overdispersion) 21.705 16.142
## Precision for ID 6.406 3.954
##
## Expected number of effective parameters(std dev): 158.70(25.24)
## Number of equivalent replicates : 4.234
##
## Marginal Likelihood: -2769.08Elevation is not significant. Drop it from the model.
formula2tr = nT ~ f(ID, model = "besag", graph = g) + l2pop + tr
model2tr = inla(formula = formula2tr, family = "nbinomial", E = area/10^6, data = spdf@data)
summary(model2tr)##
## Call:
## c("inla(formula = formula2tr, family = \"nbinomial\", data = spdf@data, ", " E = area/10^6)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.3160 1.8002 0.0570 2.1732
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -3.6471 0.1042 -3.8535 -3.6465 -3.4439 -3.6454 0
## l2pop 0.0897 0.0312 0.0285 0.0897 0.1510 0.0896 0
## tr -0.0267 0.0035 -0.0335 -0.0267 -0.0198 -0.0267 0
##
## Random effects:
## Name Model
## ID Besags ICAR model
##
## Model hyperparameters:
## mean sd
## size for the nbinomial observations (overdispersion) 16.460 2.4027
## Precision for ID 4.367 0.9146
## 0.025quant 0.5quant
## size for the nbinomial observations (overdispersion) 12.186 16.311
## Precision for ID 2.895 4.254
## 0.975quant mode
## size for the nbinomial observations (overdispersion) 21.608 16.039
## Precision for ID 6.465 4.026
##
## Expected number of effective parameters(std dev): 156.84(24.97)
## Number of equivalent replicates : 4.285
##
## Marginal Likelihood: -2758.68formula2ts = nT ~ f(ID, model = "besag", graph = g) + l2pop + ts
model2ts = inla(formula = formula2ts, family = "nbinomial", E = area/10^6, data = spdf@data)
summary(model2ts)##
## Call:
## c("inla(formula = formula2ts, family = \"nbinomial\", data = spdf@data, ", " E = area/10^6)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.2447 1.7216 0.0480 2.0143
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -3.6396 0.1050 -3.8477 -3.6390 -3.4350 -3.6378 0
## l2pop 0.0897 0.0311 0.0287 0.0897 0.1510 0.0896 0
## ts -0.5859 0.0774 -0.7374 -0.5861 -0.4335 -0.5866 0
##
## Random effects:
## Name Model
## ID Besags ICAR model
##
## Model hyperparameters:
## mean sd
## size for the nbinomial observations (overdispersion) 16.35 2.388
## Precision for ID 4.40 0.931
## 0.025quant 0.5quant
## size for the nbinomial observations (overdispersion) 12.106 16.208
## Precision for ID 2.905 4.284
## 0.975quant mode
## size for the nbinomial observations (overdispersion) 21.469 15.938
## Precision for ID 6.539 4.051
##
## Expected number of effective parameters(std dev): 155.91(25.09)
## Number of equivalent replicates : 4.31
##
## Marginal Likelihood: -2756.18formula2tri = nT ~ f(ID, model = "besag", graph = g) + l2pop + tri
model2tri = inla(formula = formula2tri, family = "nbinomial", E = area/10^6, data = spdf@data)
summary(model2tri)##
## Call:
## c("inla(formula = formula2tri, family = \"nbinomial\", data = spdf@data, ", " E = area/10^6)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.4896 1.6015 0.0781 2.1692
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -3.6612 0.1033 -3.8659 -3.6607 -3.4597 -3.6596 0
## l2pop 0.0908 0.0312 0.0296 0.0908 0.1522 0.0907 0
## tri -0.0780 0.0102 -0.0980 -0.0780 -0.0579 -0.0780 0
##
## Random effects:
## Name Model
## ID Besags ICAR model
##
## Model hyperparameters:
## mean sd
## size for the nbinomial observations (overdispersion) 16.487 2.4100
## Precision for ID 4.354 0.9108
## 0.025quant 0.5quant
## size for the nbinomial observations (overdispersion) 12.197 16.340
## Precision for ID 2.889 4.242
## 0.975quant mode
## size for the nbinomial observations (overdispersion) 21.644 16.073
## Precision for ID 6.444 4.015
##
## Expected number of effective parameters(std dev): 157.22(25.00)
## Number of equivalent replicates : 4.274
##
## Marginal Likelihood: -2757.65rSR = r
values(rSR) = (exp(model2tr$summary.random$ID$mean) - 1) * 100
range(values(rSR))## [1] -42.58884 74.90679rSS = r
values(rSS) = (exp(model2ts$summary.random$ID$mean) - 1) * 100
range(values(rSS))## [1] -41.92680 74.23533rSI = r
values(rSI) = (exp(model2tri$summary.random$ID$mean) - 1) * 100
range(values(rSI))## [1] -42.72366 75.86747rng = seq(-60, 80, 20)
rngL = paste(rng, '%', sep = "")
cr = rev(brewer.pal(8, "RdBu"))
cr = cr[-(1)]
#terrain roughtness
p5 = levelplot(rSR, margin = TRUE,
sub = "Tornado Reports (Above/Below Regional Average) ",
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at = rng, labels = rngL, col = cr),
par.settings = list(fontsize = list(text = 13)))
p5 +
latticeExtra::layer(sp.polygons(bndryCP, col = gray(.85), lwd = 1)) +
latticeExtra::layer(sp.polygons(bndrySP, col = gray(.15), lwd = 1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(type = 1),
offset = c(-95.5, 41.2),
scale = .75)})
Figure 7 Tornado Reports using Terrain Roughtness (Above/Below Regional Average)
#terrain slope
p6 = levelplot(rSS, margin = TRUE,
sub = "Tornado Reports (Above/Below Regional Average) ",
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at = rng, labels = rngL, col = cr),
par.settings = list(fontsize = list(text = 13)))
p6 +
latticeExtra::layer(sp.polygons(bndryCP, col = gray(.85), lwd = 1)) +
latticeExtra::layer(sp.polygons(bndrySP, col = gray(.15), lwd = 1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(type = 1),
offset = c(-95.5, 41.2),
scale = .75)})
Figure 8 Tornado Reports using Terrain Slope (Above/Below Regional Average)
#TRI
p7 = levelplot(rSI, margin = TRUE,
sub = "Tornado Reports (Above/Below Regional Average) ",
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at = rng, labels = rngL, col = cr),
par.settings = list(fontsize = list(text = 13)))
p7 +
latticeExtra::layer(sp.polygons(bndryCP, col = gray(.85), lwd = 1)) +
latticeExtra::layer(sp.polygons(bndrySP, col = gray(.15), lwd = 1)) +
latticeExtra::layer({SpatialPolygonsRescale(layout.north.arrow(type = 1),
offset = c(-95.5, 41.2),
scale = .75)})
Figure 9 Tornado Reports using Terrain Slope (Above/Below Regional Average)
Function for plotting the density of the marginal term.
ggplotmargin <- function(x, type, effect, xlab, ylab = "Posterior Density",
int.value = c(value = 0, 5, 95),
color = c("red", "gray", "gray")){
xx = as.data.frame(inla.smarginal(x[[paste("marginals", type, sep=".")]][[effect]]))
out = ggplot(xx, aes(x, y)) + geom_line(size = 1) + ylab(ylab) + xlab(xlab)
if(length(int.value) == 0) int.value = 0
int.value = lapply(int.value, function(x) if(is.character(x))
type.convert(x, as.is = TRUE) else x)
int.value = lapply(int.value, function(x) if(is.character(x))
lapply(strsplit(x, "=")[[1]], type.convert, as.is = TRUE) else x)
nx = names(int.value)
if(!is.null(nx))
for(i in which(nx != "")) int.value[[i]] = list(nx[i], int.value[[i]])
int.value = sapply(int.value, function(x) {
if(is.numeric(x)) xx$x[which.max(cumsum(xx$y)/sum(xx$y) >= as.numeric(x/100))]
else switch(x[[1]],
mean = sum(xx$y*xx$x)/sum(xx$y),
median = xx$x[which.max(cumsum(xx$y)/sum(xx$y) >=.5)],
mode = xx$x[which.max(xx$y)],
value = x[[2]],
zero = 0)})
if(length(color) <= length(int.value)) color = rep(color, length = length(int.value))
for(i in 1:length(int.value)) out = out + geom_vline(xintercept = int.value[i], color = color[i])
out
}#terrain roughness
results1 = model2tr
results1$marginals.fixed$tr[, 1] = (exp(-results1$marginals.fixed$tr[, 1]) - 1) * 100
results1$marginals.fixed$l2pop[, 1] = (exp(results1$marginals.fixed$l2pop[, 1]) - 1) * 100
#terrain slope
results2 = model2ts
results2$marginals.fixed$ts[, 1] = (exp(-results2$marginals.fixed$ts[, 1]) - 1) * 100
results2$marginals.fixed$l2pop[, 1] = (exp(results2$marginals.fixed$l2pop[, 1]) - 1) * 100
#TRI
results3 = model2tri
results3$marginals.fixed$tri[, 1] = (exp(-results3$marginals.fixed$tri[, 1]) - 1) * 100
results3$marginals.fixed$l2pop[, 1] = (exp(results3$marginals.fixed$l2pop[, 1]) - 1) * 100
f8a = ggplotmargin(results1, type = "fixed", effect = "l2pop",
xlab = "% increase in tornado reports\n per doubling of population")
f8b = ggplotmargin(results1, type = "fixed", effect = "tr",
xlab = "% decrease in tornado reports\n per meter increase in elevation roughness")
f8c = ggplotmargin(results2, type = "fixed", effect = "ts",
xlab = "% decrease in tornado reports\n per degree increase in elevation slope")
f8d = ggplotmargin(results3, type = "fixed", effect = "tri",
xlab = "% decrease in tornado reports\n per unit increase in TRI")
source("multiplot.txt")
mat2 = matrix(c(1, 2, 3, 4), nrow = 2, byrow = TRUE)
f8a = f8a + ggtitle("(a)") + theme(plot.title = element_text(hjust = 0))
f8b = f8b + ggtitle("(b)") + theme(plot.title = element_text(hjust = 0))
f8c = f8c + ggtitle("(c)") + theme(plot.title = element_text(hjust = 0))
f8d = f8d + ggtitle("(d)") + theme(plot.title = element_text(hjust = 0))
multiplot(f8a, f8b, f8c, f8d, layout = mat2)
Figure 10 Posterior densities of the population and roughness effects.