Tornado occurrence in Midwest
Point Model for Log-Gaussian Cox-processes
(Terrain Roughness, Elevation, Distance from Roads, Distance from Cities, Population Density)
Contact: jmh09r@my.fsu.edu
Load Libraries
library("raster")## Warning: package 'raster' was built under R version 3.1.3## Loading required package: splibrary("gridExtra")## Loading required package: gridlibrary("ggplot2")## Warning: package 'ggplot2' was built under R version 3.1.3library("rgdal")## Warning: package 'rgdal' was built under R version 3.1.3## 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/John/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/John/Documents/R/win-library/3.1/rgdal/projlibrary("rasterVis")## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
##
## Attaching package: 'latticeExtra'
##
## The following object is masked from 'package:ggplot2':
##
## layerlibrary("RColorBrewer")
library("wesanderson")
library("mapproj")## Loading required package: mapslibrary("maptools")## Checking rgeos availability: TRUElibrary("rgeos")## rgeos version: 0.3-8, (SVN revision 460)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Polygon checking: TRUElibrary("INLA")## Loading required package: Matrix
## Loading required package: splineslibrary("geostatsp")##
## Attaching package: 'geostatsp'
##
## The following object is masked from 'package:INLA':
##
## inla.modelslibrary("spatstat")##
## spatstat 1.41-1 (nickname: 'Ides of March')
## For an introduction to spatstat, type 'beginner'
##
## Attaching package: 'spatstat'
##
## The following object is masked from 'package:lattice':
##
## panel.histogram
##
## The following objects are masked from 'package:raster':
##
## area, rotate, shiftGet the tornado data
download.file(url = "http://www.spc.noaa.gov/gis/svrgis/zipped/tornado.zip",
destfile = "tornado.zip")
#download.file("http://myweb.fsu.edu/jelsner/data/tornado.zip",
# "tornado.zip", mode = "wb")
unzip("tornado.zip")Read Tornado Data
Torn = readOGR(dsn = "torn", layer = "torn",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "torn", layer: "torn"
## with 58959 features
## It has 21 fieldsTorn$Year = as.integer(Torn$yr)
Torn$Month = as.integer(Torn$mo)
Torn$EF = as.integer(Torn$mag)
Torn$Date = as.Date(Torn$date)
Torn$SLON = as.numeric(Torn$slon)
Torn$SLAT = as.numeric(Torn$slat)
TornE = as.data.frame(Torn)
TornE = subset(TornE, EF >= 0 & Year >= 1994)
TornE$EFf = as.factor(TornE$EF)Domain
Raster
r = raster(xmn = -102, xmx = -95,
ymn = 36, ymx = 42,
resolution = .25)
sp = as(r, 'SpatialPolygons')
Bounds = gUnaryUnion(sp, id = NULL)Point Projection
coordinates(TornE) = ~ SLON + SLAT
ll = "+proj=longlat +ellps=GRS80"
proj4string(TornE) = ll
TornE = spTransform(TornE, CRS(proj4string(Bounds)))Crop TornE to Bounds
TornE.pts = crop(TornE, Bounds)Distance to Cities
Get Cities Data
download.file("http://myweb.fsu.edu/jelsner/data/ci08au12.zip",
"ci08au12.zip", mode = "wb")
unzip("ci08au12.zip")To use the “distmap()” spatstat is needed. Define window and convert cities to .ppp object
Domain.w = as(Bounds, 'owin')
Cities = readOGR(dsn = ".", layer = "ci08au12",
p4s = "+proj=longlat +ellps=GRS80")## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "ci08au12"
## with 43603 features
## It has 23 fieldsCities = spTransform(Cities, CRS(proj4string(Bounds)))
Cities = subset(Cities, Cities$POP_1990 > 1000)
Cities.ppp = as(Cities, "ppp")
Cities.ppp = unmark(Cities.ppp[Domain.w])
unitname(Cities.ppp) = c("meter", "meters")Create a distance raster with distances to nearest town.
Zc = distmap(Cities.ppp)
Zc.r = raster(Zc)
plot(Zc.r)
Distance to Roads
Get Roads Data, specify projection, convert to .psp, and generate “Distance to Roads” raster.
download.file(url = "http://www.mapcruzin.com/fcc-wireless-shapefiles/roads.zip",
destfile = "roads.zip")
unzip("roads.zip")
Roads_US = readOGR(dsn = ".", layer = "roadtrl020",
p4s = "+proj=longlat +ellps=GRS80")## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "roadtrl020"
## with 47014 features
## It has 10 fieldsRoads_US = spTransform(Roads_US, CRS(proj4string(Bounds)))
Roads = crop(Roads_US, Bounds)
Roads.psp = as.psp(Roads)## Warning in as.psp.SpatialLinesDataFrame(Roads): 9 columns of data frame
## discardedZrd = distmap(Roads.psp)
Zrd.r = raster(Zrd)
plot(Zrd.r)
Get Population Data
Fit population density raster to domain
download.file(url = "http://myweb.fsu.edu/jelsner/data/usadens.zip",
destfile = "usadens.zip")
unzip("usadens.zip")
Pop = raster("usadens/usads00g/w001001.adf")
Bounds1 = spTransform(Bounds, CRS(projection(Pop)))
Pop = crop(Pop, Bounds1)
Pop2 = resample(aggregate(Pop, fact = c(nrow(r), ncol(r)), fun = mean), r)
pop = reclassify(Pop2, cbind(Pop2[Pop2<0], 6))
Pd.r = projectRaster(pop, crs = proj4string(Bounds))
plot(Pd.r)
Get DEM
#download.file(url = "http://www.viewfinderpanoramas.org/DEM/TIF15/15-H.ZIP",
# destfile = "15-H.ZIP", mode = "wb")
download.file(url = "http://myweb.fsu.edu/jelsner/data/15-H.tif.zip",
destfile = "15-H.tif.zip", mode = "wb")
unzip("15-H.tif.zip")Create ELV, TR, and TRI rasters
Elev = raster("15-H.tif")
Bounds2 = spTransform(Bounds, CRS(projection(Elev)))
Elev2 = crop(Elev, Bounds2)
ELV = projectRaster(Elev2, crs = proj4string(Bounds))
TR = terrain(ELV, opt = 'roughness')
TRI = terrain(ELV, opt = 'TRI')
plot(ELV)
plot(TR)
plot(TRI)
Model for Log-Gaussian Cox-processes (LGCP)
Via geostatsp and INLA
List potential covariates
covList = list(TR = TR,
TRI = TRI,
ELV = ELV,
Pd.r = Pd.r,
Zc.r = Zc.r,
Zrd.r = Zrd.r)Model0
No Covariates
Model0 = lgcp(TornE.pts ~ 1,
data = TornE.pts,
family = "nbinomial",
grid = 50,
covariates = covList,
control.compute =
list(dic = TRUE,
cpo = TRUE),
shape = 2,
buffer = 15,
border = Bounds)
Model0$parameters$summary## mean sd 0.025quant 0.5quant 0.975quant mode
## (Intercept) 4.1365223 0.0540894 4.0273773 4.1375653 4.2399779 4.1397665
## range 0.7931460 0.1327254 0.5673437 0.7809422 1.0874180 0.7562888
## sd 0.4529506 NA 0.3753045 0.4495422 0.5373912 NA
## kld meanExp meanInvLogit
## (Intercept) 6.071632e-12 63.3073 1.00008
## range NA NA NA
## sd NA NA NADefine State and County Boundaries
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)))Plot the smooth rate
Results are automatically returned as raster objects; but, used rasterVis to control ascetics.
rSR = (exp(Model0$raster$random.mean) - 1)*100
range(values(rSR))## [1] -48.98598 262.92915rng = seq(-50, 300, 50)
cr = wes_palette(name = "Zissou", n = 7,
type = "continuous")
rngL = paste(rng, '%', sep = "")
p0 = levelplot(rSR, margin = FALSE,
sub = "Tornado Reports \n (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))) +
latticeExtra::layer(sp.polygons(bndryCP,col = gray(.85), lwd = 1)) +
latticeExtra::layer(sp.polygons(bndrySP, col = "black", lwd = 2))
p0
Model1
First model using covariates explored in class project Roughness (TR) Elevation (ELV) Population Density (Pd.r) Distance to City (Zc.r)
Model1 = lgcp(TornE.pts ~ TR + ELV + Pd.r + Zc.r,
data = TornE.pts,
family = "nbinomial",
grid = 50,
covariates = covList,
control.compute =
list(dic = TRUE,
cpo = TRUE),
shape = 2,
buffer = 15,
border = Bounds)
Model1$parameters$summary## mean sd 0.025quant 0.5quant
## (Intercept) 4.4786083324 0.1695402788 4.1383109753 4.481139524
## TR -0.0245072965 0.0062792482 -0.0368251648 -0.024512546
## ELV 0.0002902107 0.0002017728 -0.0001034052 0.000289181
## Pd.r 0.0045756873 0.0019523001 0.0007078041 0.004586603
## Zc.r -1.3232798120 0.2899174914 -1.8948805645 -1.322540279
## range 0.8272534981 0.1620082166 0.5684675552 0.806194695
## sd 0.3932440259 NA 0.3223943906 0.390240709
## 0.975quant mode kld meanExp
## (Intercept) 4.8048501075 4.4863680146 2.410178e-12 89.0869129
## TR -0.0121724852 -0.0245229054 7.011889e-14 0.9896042
## ELV 0.0006890577 0.0002870869 8.879654e-13 1.0149214
## Pd.r 0.0083798546 0.0046088394 3.440154e-13 1.0197471
## Zc.r -0.7563048067 -1.3210285628 9.906552e-14 0.2728044
## range 1.2009571198 0.7626769084 NA NA
## sd 0.4763431153 NA NA NA
## meanInvLogit
## (Intercept) 1.0040801
## TR 0.5010287
## ELV 0.5073928
## Pd.r 0.5087616
## Zc.r 0.2119988
## range NA
## sd NAModel2
Dropping Elevation - not significant in Model1
Model2 = lgcp(TornE.pts ~ TR + Pd.r + Zc.r,
data = TornE.pts,
family = "nbinomial",
grid = 50,
covariates = covList,
control.compute =
list(dic = TRUE,
cpo = TRUE),
shape = 2,
buffer = 15,
border = Bounds)
Model2$parameters$summary## mean sd 0.025quant 0.5quant
## (Intercept) 4.66759328 0.103028309 4.4621600952 4.668683863
## TR -0.02530642 0.006297685 -0.0376450264 -0.025316682
## Pd.r 0.00318527 0.001719298 -0.0002176172 0.003193579
## Zc.r -1.20697803 0.279593150 -1.7589882268 -1.206000224
## range 0.76251418 0.140535050 0.5224568302 0.750449941
## sd 0.39054342 NA 0.3198439868 0.387802765
## 0.975quant mode kld meanExp
## (Intercept) 4.866936504 4.670970450 1.255100e-12 107.4225309
## TR -0.012923452 -0.025337447 8.064369e-14 0.9887611
## Pd.r 0.006539322 0.003210653 3.260920e-13 1.0183102
## Zc.r -0.660887458 -1.204010825 9.870561e-14 0.3059223
## range 1.072598530 0.727360838 NA NA
## sd 0.471014722 NA NA NA
## meanInvLogit
## (Intercept) 1.0059211
## TR 0.5007999
## Pd.r 0.5083923
## Zc.r 0.2318229
## range NA
## sd NAModel3
Dropping Population Density - not significant in Model2
Model3 = lgcp(TornE.pts ~ TR + Zc.r,
data = TornE.pts,
family = "nbinomial",
grid = 50,
covariates = covList,
control.compute =
list(dic = TRUE,
cpo = TRUE),
shape = 2,
buffer = 15,
border = Bounds)
Model3$parameters$summary## mean sd 0.025quant 0.5quant 0.975quant
## (Intercept) 4.71668621 0.09808328 4.52151376 4.71755357 4.90695168
## TR -0.02306326 0.00616838 -0.03513662 -0.02307737 -0.01092405
## Zc.r -1.32499540 0.27169157 -1.86175131 -1.32392861 -0.79463605
## range 0.74759331 0.14039584 0.50596445 0.73622956 1.05476269
## sd 0.38763291 NA 0.31864775 0.38470254 0.46895792
## mode kld meanExp meanInvLogit
## (Intercept) 4.71934801 1.015553e-12 112.8270111 1.0063306
## TR -0.02310594 7.337450e-14 0.9909565 0.5013531
## Zc.r -1.32176247 1.060344e-13 0.2715251 0.2115996
## range 0.71509472 NA NA NA
## sd NA NA NA NAChecking for Correlation between distance metrics
cor(values(Zc.r), values(Zrd.r), use = "pairwise.complete.obs")## [1] 0.3265844Model4
Adding Distance to Roads to Model3
Model4 = lgcp(TornE.pts ~ TR + Zc.r + Zrd.r,
data = TornE.pts,
family = "nbinomial",
grid = 50,
covariates = covList,
control.compute =
list(dic = TRUE,
cpo = TRUE),
shape = 2,
buffer = 15,
border = Bounds)
Model4$parameters$summary## mean sd 0.025quant 0.5quant 0.975quant
## (Intercept) 4.81797450 0.101796871 4.61604570 4.81864201 5.01606207
## TR -0.02229298 0.006219583 -0.03446399 -0.02230817 -0.01005018
## Zc.r -0.94433666 0.282196535 -1.50156933 -0.94332874 -0.39319386
## Zrd.r -3.66565113 0.847354411 -5.33870314 -3.66257080 -2.01117496
## range 0.76559146 0.134323066 0.53246298 0.75532017 1.05834465
## sd 0.39601084 NA 0.32586398 0.39295496 0.47879067
## mode kld meanExp meanInvLogit
## (Intercept) 4.82001768 6.856795e-13 124.83103528 1.00706880
## TR -0.02233886 6.936859e-14 0.99170695 0.50154308
## Zc.r -0.94127675 1.083945e-13 0.39798454 0.28134267
## Zrd.r -3.65630942 2.474002e-14 0.03388599 0.03178276
## range 0.73614656 NA NA NA
## sd NA NA NA NAChecking for Correlation between terrain metrics
cor(values(TRI), values(TR), use = "pairwise.complete.obs")## [1] 0.9560921Model5
Swapping TRI for TR
Model5 = lgcp(TornE.pts ~ TRI + Zc.r + Zrd.r,
data = TornE.pts,
family = "nbinomial",
grid = 50,
covariates = covList,
control.compute =
list(dic = TRUE,
cpo = TRUE),
shape = 2,
buffer = 15,
border = Bounds)
Model5$parameters$summary## mean sd 0.025quant 0.5quant 0.975quant
## (Intercept) 4.79741711 0.09979778 4.5995736 4.79803464 4.99171137
## TRI -0.06406776 0.01854116 -0.1003938 -0.06409787 -0.02761197
## Zc.r -0.94102532 0.28319255 -1.5001524 -0.94003868 -0.38787278
## Zrd.r -3.65701819 0.84794878 -5.3312001 -3.65395097 -2.00134160
## range 0.77178663 0.13483647 0.5397188 0.76076118 1.06790401
## sd 0.39492911 NA 0.3238216 0.39213331 0.47575066
## mode kld meanExp meanInvLogit
## (Intercept) 4.79932108 8.917811e-13 122.28881744 1.00685853
## TRI -0.06415922 6.055324e-14 0.95001047 0.49062595
## Zc.r -0.93802942 1.083255e-13 0.39937024 0.28200583
## Zrd.r -3.64771653 2.532298e-14 0.03419534 0.03205374
## range 0.73959551 NA NA NA
## sd NA NA NA NACompare Models
Set-up Metric Tables
Models = c("Model0", "Model1",
"Model2", "Model3",
"Model4", "Model5")
#Deviance information criteria
DICs = c(Model0$inla$dic$dic, Model1$inla$dic$dic,
Model2$inla$dic$dic, Model3$inla$dic$dic,
Model4$inla$dic$dic, Model5$inla$dic$dic)
#Log of conditional predictive ordinances
LCPOs = c(-mean(log(Model0$inla$cpo$cpo)), -mean(log(Model1$inla$cpo$cpo)),
-mean(log(Model2$inla$cpo$cpo)), -mean(log(Model3$inla$cpo$cpo)),
-mean(log(Model4$inla$cpo$cpo)), -mean(log(Model5$inla$cpo$cpo)))
#Failed Observations
Fails = c(sum(Model0$inla$cpo$failure[Model0$inla$cpo$failure==1]),
sum(Model1$inla$cpo$failure[Model1$inla$cpo$failure==1]),
sum(Model2$inla$cpo$failure[Model2$inla$cpo$failure==1]),
sum(Model3$inla$cpo$failure[Model3$inla$cpo$failure==1]),
sum(Model4$inla$cpo$failure[Model4$inla$cpo$failure==1]),
sum(Model5$inla$cpo$failure[Model5$inla$cpo$failure==1]))
#Partial Failed Observations
pFails = c(sum(Model0$inla$cpo$failure[Model0$inla$cpo$failure < 1]),
sum(Model1$inla$cpo$failure[Model1$inla$cpo$failure < 1]),
sum(Model2$inla$cpo$failure[Model2$inla$cpo$failure < 1]),
sum(Model3$inla$cpo$failure[Model3$inla$cpo$failure < 1]),
sum(Model4$inla$cpo$failure[Model4$inla$cpo$failure < 1]),
sum(Model5$inla$cpo$failure[Model5$inla$cpo$failure < 1]))
#probability integral transform (PIT)
PITs = c(mean(Model0$inla$cpo$pit), mean(Model1$inla$cpo$pit),
mean(Model2$inla$cpo$pit), mean(Model3$inla$cpo$pit),
mean(Model4$inla$cpo$pit), mean(Model5$inla$cpo$pit))
Model_mets = as.data.frame(cbind(Model = Models,
DIC = round(DICs, 2),
LCPO = round(LCPOs, 3),
PIT = round(PITs, 3),
Fail = Fails,
pFail = round(pFails, 2)))
Model_mets## Model DIC LCPO PIT Fail pFail
## 1 Model0 6708.43 1.561 0.626 0 0
## 2 Model1 6679.66 1.554 0.628 0 0
## 3 Model2 6680.62 1.554 0.628 0 0
## 4 Model3 6683.87 1.555 0.628 0 0
## 5 Model4 6663.2 1.55 0.629 0 0
## 6 Model5 6663.52 1.55 0.628 0 0ggplotmargin()
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
}Density of the Marginal Terms
(Model4)
results = Model4
results$inla$marginals.fixed$TR[, 1] = (exp(-results$inla$marginals.fixed$TR[, 1]) - 1) * 100
results$inla$marginals.fixed$Zrd.r[, 1] = (exp(-results$inla$marginals.fixed$Zrd.r[, 1]) - 1) * 100
results$inla$marginals.fixed$Zc.r[, 1] = (exp(-results$inla$marginals.fixed$Zc.r[, 1]) - 1) * 100
plotTRP = ggplotmargin(results$inla, type = "fixed",
effect = "TR",
xlab = "Terrain Roughness")
plotRds = ggplotmargin(results$inla, type = "fixed",
effect = "Zrd.r",
xlab = "Distance to Roads")
plotCty = ggplotmargin(results$inla, type = "fixed",
effect = "Zc.r",
xlab = "Distance to City")
grid.arrange(plotTRP, plotRds, plotCty, ncol=2)
Plot Results
Model4 and Model5 perform equally well; Model4 is selected.
Results are automatically returned as raster objects; but, used rasterVis to control ascetics.
rM4 = (exp(Model4$raster$random.mean) - 1)*100
range(values(rM4))## [1] -44.53376 213.25455rng = seq(-45, 215, 26)
cr = wes_palette(name = "Zissou", n = 10,
type = "continuous")
rngL = paste(rng, '%', sep = "")
p1 = levelplot(rM4, margin = FALSE,
sub = "Tornado Reports \n (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))) +
latticeExtra::layer(sp.polygons(bndryCP,col = gray(.85), lwd = 1)) +
latticeExtra::layer(sp.polygons(bndrySP, col = "black", lwd = 2))
p1
Plot Results #2 (change color scheme)
Model3 as shown in previous plot; but, with color scheme more comparable to that in used in class.
cr2 = rev(brewer.pal(11, "RdBu"))
cr2 = cr2[-(6)]
p2 = levelplot(rM4, margin = FALSE,
sub = "Tornado Reports \n (Above/Below Regional Average) ",
xlab = NULL, ylab = NULL,
col.regions = cr2, at = rng,
colorkey = list(at = rng, labels = rngL, col = cr2),
par.settings = list(fontsize = list(text = 13))) +
latticeExtra::layer(sp.polygons(bndryCP,col = gray(.85), lwd = 1)) +
latticeExtra::layer(sp.polygons(bndrySP, col = "black", lwd = 2))
p2