date()
## [1] "Tue Mar 18 18:32:05 2014"
Tornadoes in Iowa.
download.file("http://myweb.fsu.edu/jelsner/data/tornado2011.zip", "tornado2011.zip",
mode = "wb")
unzip("tornado2011.zip")
library(spatstat)
##
## spatstat 1.36-0 (nickname: 'Intense Scrutiny')
## For an introduction to spatstat, type 'beginner'
library(sp)
library(rgdal)
## 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(rgeos)
## rgeos version: 0.3-3, (SVN revision 437)
## GEOS runtime version: 3.3.3-CAPI-1.7.4
## Polygon checking: TRUE
library(maptools)
## Checking rgeos availability: TRUE
library(ggplot)
## Error: there is no package called 'ggplot'
require(maps)
## Loading required package: maps
require(ggmap)
## Loading required package: ggmap
## Loading required package: ggplot2
Torn = readOGR(dsn = ".", layer = "tornado", stringsAsFactors = FALSE)
## 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)
Torn$EF = factor(Torn$FSCALE)
ia.df = map_data("state", "iowa")
ia.m = cbind(ia.df$long, ia.df$lat)
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")
ia.SPS = spTransform(ia.SPS, CRS(proj4string(Torn)))
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
den = density(TornI.ppp[TornI.ppp$marks == 1], dimyx = c(256, 256), adjust = 0.8)
plot(den, main = "Spatial Density of Recorded Tornado Touchdowns \nin Iowa, 1954-2011")
points(TornI.ppp, pch = 16, cex = 0.5)
plot(Gest(TornI.ppp))
## 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 = 250), main = "Nearest Neighbor Distance Function for \n1954-2011 Iowa Tornados, with Simulated Clustering Envelope")
## Generating 250 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24
## .26.28.30.32.34.36.38.40.42.44.46.48
## .50.52.54.56.58.60.62.64.66.68.70.72
## .74.76.78.80.82.84.86.88.90.92.94.96
## .98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144
## .146.148.150.152.154.156.158.160.162.164.166.168
## .170.172.174.176.178.180.182.184.186.188.190.192
## .194.196.198.200.202.204.206.208.210.212.214.216
## .218.220.222.224.226.228.230.232.234.236.238.240
## .242.244.246.248. 250.
##
## 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
print("As the simulation envelope does not overlap with the observed nearest neighbor distribution of tornadoes, there is evidence to support the existence of clustering.")
## [1] "As the simulation envelope does not overlap with the observed nearest neighbor distribution of tornadoes, there is evidence to support the existence of clustering."
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)
towns.ppp = unmark(US.ppp[W])
plot(towns.ppp, pch = 16, cex = 0.5, main = "Location of Towns in Iowa \nwith Populations Exceeding 1000")
Zc = distmap(towns.ppp)
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 = "black") + ylab(expression(paste(hat(rho),
"(Tornado Reports/km", {
}^2, ")"))) + xlab("Distance from Nearest Town Center (km)") + ggtitle("Spatial Density of Iowa Tornado Reports \nas a Function of Distance from Town Centers") +
theme_bw()
marks(towns.ppp) = NULL
marks(TornI.ppp) = NULL
TornTown.ppp = superimpose(towns.ppp, TornI.ppp)
## Warning: data contain duplicated points
marks(TornTown.ppp) = as.factor(c(rep("city", towns.ppp$n), rep("tornado", TornI.ppp$n)))
plot(TornTown.ppp, cex = 0.5, main = "Iowa Towns (Circles) and \nTornado Touchdowns (Triangles)")
## city tornado
## 1 2
plot(Kcross(TornTown.ppp, "city", "tornado", correction = "border"), xlim = c(0,
10000), main = "Kcross Function, Iowa Torns and Tornadoes")
## 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(TornTown.ppp, 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
print("From the lack of overlap between the simulation envelope and the observed K-Cross function, it appears that much of the clustering seen in tornado touchdowns in Iowa is explained by city locations.")
## [1] "From the lack of overlap between the simulation envelope and the observed K-Cross function, it appears that much of the clustering seen in tornado touchdowns in Iowa is explained by city locations."