date()
## [1] "Wed Mar 19 21:53:45 2014"
Tornadoes in Iowa.
(1.) Create a ppp object for Iowa tornado touchdown reports (all) using the period 1954-2011. Use the touchdown locations in http://myweb.fsu.edu/jelsner/data/tornado2011.zip and the Iowa state boundary as the window.
library(rgdal)
## Loading required package: sp
## rgdal: version: 0.8-16, (SVN revision 498)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.0/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.0/Resources/library/rgdal/proj
library(spatstat)
##
## spatstat 1.36-0 (nickname: 'Intense Scrutiny')
## For an introduction to spatstat, type 'beginner'
library(maptools)
## Checking rgeos availability: TRUE
library(maps)
library(ggplot2)
download.file("http://myweb.fsu.edu/jelsner/data/tornado2011.zip", "tornado.zip",
mode = "wb")
unzip("tornado.zip")
Torn = readOGR(dsn = ".", layer = "tornado")
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "tornado"
## with 56221 features and 28 fields
## Feature type: wkbPoint with 2 dimensions
Torn = subset(Torn, FSCALE >= 0 & YEAR >= 1954 & YEAR <= 2011)
Torn$EF = factor(Torn$FSCALE)
class(Torn)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
ia.df = map_data("state", "iowa")
ia.m = cbind(ia.df$long, ia.df$lat)
library(sp)
ia.P = Polygon(ia.m)
ia.PS = Polygons(list(ia.P), 1)
ia.SPS = SpatialPolygons(list(ia.PS))
proj4string(ia.SPS) = CRS("+proj=longlat +ellps=WGS84")
library(rgdal)
ia.SPS = spTransform(ia.SPS, CRS(proj4string(Torn)))
library(rgeos)
## rgeos version: 0.3-3, (SVN revision 437)
## GEOS runtime version: 3.3.3-CAPI-1.7.4
## Polygon checking: TRUE
gIsValid(ia.SPS)
## [1] TRUE
xy = fortify(ia.SPS)
W = owin(poly = list(x = rev(xy$long[-1]), y = rev(xy$lat[-1])))
Torn.ppp = as.ppp(Torn["EF"])
TornI.ppp = Torn.ppp[W]
unitname(TornI.ppp) = c("Meter", "Meters")
summary(TornI.ppp)
## Marked planar point pattern: 2114 points
## Average intensity 1.47e-08 points per square Meter
##
## *Pattern contains duplicated points*
##
## Coordinates are given to 2 decimal places
## i.e. rounded to the nearest multiple of 0.01 Meters
##
## Multitype:
## frequency proportion intensity
## 0 925 0.43800 6.42e-09
## 1 660 0.31200 4.58e-09
## 2 399 0.18900 2.77e-09
## 3 90 0.04260 6.25e-10
## 4 35 0.01660 2.43e-10
## 5 5 0.00237 3.47e-11
##
## Window: polygonal boundary
## single connected closed polygon with 255 vertices
## enclosing rectangle: [-49514, 483273] x [161758, 507522] Meters
## Window area = 1.43999e+11 square Meters
## Unit of length: 1 Meter
plot(TornI.ppp)
## 0 1 2 3 4 5
## 1 2 3 4 5 6
(2.) Create a spatial density map with the events plotted as points.
den = density(TornI.ppp[TornI.ppp$marks == 1], dimyx = c(256, 256), adjust = 0.8)
plot(den)
(3.) Plot the G function (nearest neighbor distance function) and test for clustering using a simulation envelope.
G = Gest(TornI.ppp)
G
## Function value object (class 'fv')
## for the function r -> G(r)
## Entries:
## id label description
## -- ----- -----------
## r r distance argument r
## theo G[pois](r) theoretical Poisson G(r)
## han hat(G)[han](r) Hanisch estimate of G(r)
## rs hat(G)[bord](r) border corrected estimate of G(r)
## km hat(G)[km](r) Kaplan-Meier estimate of G(r)
## hazard hat(lambda)[km](r) Kaplan-Meier estimate of hazard function lambda(r)
## --------------------------------------
##
## Default plot formula: . ~ r
##
## Recommended range of argument r: [0, 6727.2]
## Available range of argument r: [0, 15800]
## Unit of length: 1 Meter
plot(G)
## lty col key label meaning
## km 1 1 km hat(G)[km](r) Kaplan-Meier estimate of G(r)
## rs 2 2 rs hat(G)[bord](r) border corrected estimate of G(r)
## han 3 3 han hat(G)[han](r) Hanisch estimate of G(r)
## theo 4 4 theo G[pois](r) theoretical Poisson G(r)
plot(envelope(TornI.ppp, fun = Gest, nsim = 99))
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
## 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
## 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
## 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
## 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
## lty col key label
## obs 1 1 obs hat(G)[obs](r)
## theo 2 2 theo G[theo](r)
## hi 1 8 hi hat(G)[hi](r)
## lo 1 8 lo hat(G)[lo](r)
## meaning
## obs observed value of G(r) for data pattern
## theo theoretical value of G(r) for CSR
## hi upper pointwise envelope of G(r) from simulations
## lo lower pointwise envelope of G(r) from simulations
(4.) Create a ppp object of the cities and towns in the state with population exceeding 1000.
download.file("http://myweb.fsu.edu/jelsner/data/ci08au12.zip", "ci08au12.zip",
mode = "wb")
unzip("ci08au12.zip")
US = readOGR(dsn = ".", layer = "ci08au12", p4s = "+proj=longlat +ellps=WGS84")
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "ci08au12"
## with 43603 features and 23 fields
## Feature type: wkbPoint with 2 dimensions
UST = spTransform(US, CRS(proj4string(Torn)))
USt = subset(UST, UST$POP_1990 > 1000)
US.ppp = as.ppp(USt)
W = as.owin(TornI.ppp)
CP.ppp = unmark(US.ppp[W])
plot(CP.ppp, pch = 20)
(5.) Plot an estimate of the spatial density of tornado reports as a function of distance from nearest town.
Zc = distmap(CP.ppp)
plot(Zc)
rhat = rhohat(TornI.ppp[TornI.ppp$marks == 1], Zc, adjust = 2, smoother = "kernel",
method = "transform")
## Warning: data contain duplicated points
m2km = 0.001
km2 = 1e+06
dist = rhat$Zc * m2km
rho = rhat$rho * km2
hi = rhat$hi * km2
lo = rhat$lo * km2
Rho.df = data.frame(dist = dist, rho = rho, hi = hi, lo = lo)
ggplot(Rho.df) + geom_ribbon(aes(x = dist, ymin = lo, ymax = hi), alpha = 0.3) +
geom_line(aes(x = dist, y = rho), col = "red") + ylab(expression(paste(hat(rho),
"(Tornado Reports/km", {
}^2, ")"))) + xlab("Distance From Nearest Town Center (km)") + theme_bw()
(6.) Create a multitype point pattern object by combining town locations and tornado touchdown locations.
marks(CP.ppp) = NULL
marks(TornI.ppp) = NULL
X = superimpose(CP.ppp, TornI.ppp)
## Warning: data contain duplicated points
marks(X) = as.factor(c(rep("city", CP.ppp$n), rep("tornado", TornI.ppp$n)))
plot(X)
## city tornado
## 1 2
(7.) Plot the Kcross function and test for 'interaction' between town and tornado reports for distances less than 10 km. Use the envelope function to create 99 simulations of the Kross function
class(X)
## [1] "ppp"
Kcross(X, "city", "tornado", correction = "border")
## Function value object (class 'fv')
## for the function r -> Kcross["city", "tornado"](r)
## Entries:
## id label description
## -- ----- -----------
## r r distance argument r
## theo {K[list(city,tornado)]^{pois}}(r) theoretical Poisson Kcross["city", "tornado"](r)
## border {hat(K)[list(city,tornado)]^{bord}}(r) border-corrected estimate of Kcross["city", "tornado"](r)
## --------------------------------------
##
## Default plot formula: . ~ r
##
## Recommended range of argument r: [0, 86441]
## Available range of argument r: [0, 86441]
## Unit of length: 1 Meter
plot(Kcross(X, "city", "tornado", correction = "border"), xlim = c(0, 10000))
## lty col key label
## border 1 1 border {hat(K)[list(city,tornado)]^{bord}}(r)
## theo 2 2 theo {K[list(city,tornado)]^{pois}}(r)
## meaning
## border border-corrected estimate of Kcross["city", "tornado"](r)
## theo theoretical Poisson Kcross["city", "tornado"](r)
plot(envelope(X, correction = "border", fun = Kcross, nsim = 99))
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
## 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
## 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
## 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
## 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
## lty col key label
## obs 1 1 obs {hat(K)[list(city,tornado)]^{obs}}(r)
## theo 2 2 theo {K[list(city,tornado)]^{theo}}(r)
## hi 1 8 hi {hat(K)[list(city,tornado)]^{hi}}(r)
## lo 1 8 lo {hat(K)[list(city,tornado)]^{lo}}(r)
## meaning
## obs observed value of Kcross["city", "tornado"](r) for data pattern
## theo theoretical value of Kcross["city", "tornado"](r) for CSR
## hi upper pointwise envelope of Kcross["city", "tornado"](r) from simulations
## lo lower pointwise envelope of Kcross["city", "tornado"](r) from simulations
(8.) (Optional) Create and interpret a point pattern model of Iowa tornadoes.