Applied Spatial Statistics: Problem Set # 2

Tyler Fricker

date()

## [1] "Tue Mar 18 16:21:01 2014"

Due Date: March 19, 2014

Total Points: 50

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.

require(rgdal)

## Loading required package: rgdal Loading required package: sp rgdal:
## version: 0.8-14, (SVN revision 496) 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

require(spatstat)

## Loading required package: spatstat
## 
## spatstat 1.36-0 (nickname: 'Intense Scrutiny') For an introduction to
## spatstat, type 'beginner'

require(ggmap)

## Loading required package: ggmap Loading required package: 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-2, (SVN revision 413M) 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])))
library(maptools)

## Checking rgeos availability: TRUE

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

1. 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)

1. 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

1. 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)

1. 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 = "black") + ylab(expression(paste(hat(rho), 
    "(Tornado Reports/km", {
    }^2, ")"))) + xlab("Distance from Nearest Town Center (km)") + theme_bw()

1. 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

1. 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

1. (Optional) Create and interpret a point pattern model of Iowa tornadoes.