download.file("http://myweb.fsu.edu/jelsner/data/tornado2011.zip", "tornado.zip",
mode = "wb")
unzip("tornado.zip")
require(rgdal)
## Loading required package: rgdal
## Warning: package 'rgdal' was built under R version 3.0.3
## 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.10.1, released 2013/08/26
## Path to GDAL shared files: C:/Users/Rui_W/Documents/R/win-library/3.0/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/Rui_W/Documents/R/win-library/3.0/rgdal/proj
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)
Torn$EF = factor(Torn$FSCALE)
class(Torn)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
## Specify Iowa state border as the window.
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.0.3
Io.df = map_data("state", "Iowa")
Io.m = cbind(Io.df$long, Io.df$lat)
library(sp)
Io.P = Polygon(Io.m)
Io.PS = Polygons(list(Io.P), 1)
Io.SPS = SpatialPolygons(list(Io.PS))
proj4string(Io.SPS) = CRS("+proj=longlat +ellps=WGS84")
library(rgdal)
Io.SPS = spTransform(Io.SPS, CRS(proj4string(Torn)))
library(rgeos)
## Warning: package 'rgeos' was built under R version 3.0.3
## rgeos version: 0.3-3, (SVN revision 437)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Polygon checking: TRUE
library(spatstat)
## Warning: package 'spatstat' was built under R version 3.0.3
##
## spatstat 1.36-0 (nickname: 'Intense Scrutiny')
## For an introduction to spatstat, type 'beginner'
gIsValid(Io.SPS)
## [1] TRUE
xy = fortify(Io.SPS)
W = owin(poly = list(x = rev(xy$long[-1]), y = rev(xy$lat[-1])))
library(maptools)
## Warning: package 'maptools' was built under R version 3.0.3
## Checking rgeos availability: TRUE
Torn.ppp = as.ppp(Torn["EF"])
TornK.ppp = Torn.ppp[W]
unitname(TornK.ppp) = c("meter", "meters")
summary(TornK.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(split(TornK.ppp))
## Spatial intensity of tornadoes.
den = density(TornK.ppp[TornK.ppp$marks == 1], dimyx = c(256, 256), adjust = 0.8)
plot(den)
str(den)
## List of 10
## $ v : num [1:256, 1:256] NA NA NA NA NA NA NA NA NA NA ...
## $ dim : int [1:2] 256 256
## $ xrange: num [1:2] -49514 483273
## $ yrange: num [1:2] 161758 507522
## $ xstep : num 2081
## $ ystep : num 1351
## $ xcol : num [1:256] -48473 -46392 -44311 -42230 -40148 ...
## $ yrow : num [1:256] 162433 163784 165135 166485 167836 ...
## $ type : chr "real"
## $ units :List of 3
## ..$ singular : chr "meter"
## ..$ plural : chr "meters"
## ..$ multiplier: num 1
## ..- attr(*, "class")= chr "units"
## - attr(*, "class")= chr "im"
## - attr(*, "sigma")= num 34576
## plot G function
plot(Gest(TornK.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(TornK.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
## spatial analysis with nearest city
library(rgdal)
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
summary(US)
## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## coords.x1 -178.88 145.78
## coords.x2 -14.36 71.27
## Is projected: FALSE
## proj4string : [+proj=longlat +ellps=WGS84]
## Number of points: 43603
## Data attributes:
## ST_FIPS SFIPS COUNTY_FIP CFIPS
## 21 : 3185 21 : 3185 003 : 981 003 : 978
## 13 : 1934 13 : 1934 019 : 789 019 : 787
## 48 : 1907 48 : 1910 013 : 782 013 : 782
## 17 : 1771 17 : 1771 031 : 740 031 : 740
## 12 : 1717 12 : 1717 025 : 715 025 : 716
## (Other):33041 (Other):33038 (Other):39406 (Other):39412
## NA's : 48 NA's : 48 NA's : 190 NA's : 188
## PL_FIPS ID NAME ELEVATION
## 50000 : 568 Number : 14 FAIRVIEW : 37 10 : 222
## 0 : 82 2115368: 11 FRANKLIN : 33 0 : 196
## 68232 : 82 0602360: 9 GEORGETOWN: 30 20 : 157
## 39952 : 67 0812399: 8 MARION : 30 5 : 144
## 49656 : 61 2105149: 8 UNION : 30 15 : 128
## (Other):37841 (Other):43325 CLINTON : 29 (Other):26672
## NA's : 4902 NA's : 228 (Other) :43414 NA's :16084
## POP_1990 POPULATION ST STATE
## Min. : 0 0 : 1480 KY : 3185 KENTUCKY: 3185
## 1st Qu.: 188 1604 : 84 GA : 1933 GEORGIA : 1934
## Median : 906 318 : 72 TX : 1902 TEXAS : 1902
## Mean : 7664 2490 : 67 IL : 1770 ILLINOIS: 1770
## 3rd Qu.: 2968 111 : 66 FL : 1717 FLORIDA : 1717
## Max. :7322564 (Other):36028 NY : 1702 NEW YORK: 1702
## NA's : 5806 (Other):31394 (Other) :31393
## WARNGENLEV WARNGENTYP WATCH_WARN ZWATCH_WAR PROG_DISC
## 0 : 1 ` : 1 3 :21419 Min. : 0 Min. : -20
## 1 : 6825 + : 25 2 :15297 1st Qu.: 0 1st Qu.: 400
## 2 :15322 = : 1 1 : 6791 Median : 0 Median : 1450
## 23: 1 0 : 1753 2+ : 15 Mean : 0 Mean : 2390
## 3 :21453 2 : 1 01 : 10 3rd Qu.: 0 3rd Qu.: 3350
## 9 : 1 NA's:41822 (Other): 13 Max. :7980 Max. :606900
## NA's : 58
## ZPROG_DISC COMBOFLAG LAND_WATER RECNUM
## Min. : -20 Min. : 0 0 : 1753 Min. : 0.0
## 1st Qu.: 540 1st Qu.: 0 L : 769 1st Qu.: 0.0
## Median : 1960 Median : 0 NA's:41081 Median : 0.0
## Mean : 2806 Mean : 0 Mean : 43.4
## 3rd Qu.: 4150 3rd Qu.: 0 3rd Qu.: 0.0
## Max. :13470 Max. :140 Max. :1785.0
##
## LON LAT
## Min. :-178.9 Min. :-910000
## 1st Qu.: -97.0 1st Qu.: 35
## Median : -87.7 Median : 39
## Mean : -90.9 Mean : 18
## 3rd Qu.: -81.8 3rd Qu.: 42
## Max. : 145.8 Max. : 71
##
## select city population > 1000
UST = spTransform(US, CRS(proj4string(Torn)))
USt = subset(UST, UST$POP_1990 > 1000)
US.ppp = as.ppp(USt)
W = as.owin(TornK.ppp)
CP.ppp = unmark(US.ppp[W])
plot(CP.ppp, pch = 20)
plot(split(TornK.ppp)$"3", add = TRUE, col = "red")
## here use TornK.ppp$3 as example, but it could be changed to 1.2.4.5 and 0
Zc = distmap(CP.ppp)
plot(Zc)
rhat = rhohat(TornK.ppp[TornK.ppp$marks == 1], Zc, adjust = 2, smoother = "kernel",
method = "transform")
## Warning: data contain duplicated points
library(ggplot2)
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()
## multitype point patten- city and tornadoes
EF3 = split(TornK.ppp)$"3"
plot(Gest(EF3))
## 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)
X = superimpose(CP.ppp, EF3)
marks(X) = as.factor(c(rep("city", CP.ppp$n), rep("tornado", EF3$n)))
library(spatstat)
Kcross(X, "city", "tornado", xlim = c(0, 10000))
## Error: 4 segments do not lie entirely inside the window.
SP = TornK.ppp
plot(SP)
## 0 1 2 3 4 5
## 1 2 3 4 5 6
plot(Gest(SP))
## 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(Kest(SP))
## Error: error in evaluating the argument 'x' in selecting a method for function 'plot': Error: 4 segments do not lie entirely inside the window.
SP.ppm = ppm(SP, trend = ~1, interaction = Strauss(r = 10), rbord = 10)
## Warning: data contain duplicated points
SP.ppm
## Stationary Strauss process
## Possible marks:
## 0 1 2 3 4 5
##
## First order terms:
## beta_0 beta_1 beta_2 beta_3 beta_4 beta_5
## 2.357e-09 2.357e-09 2.357e-09 2.357e-09 2.357e-09 2.357e-09
##
## Interaction: Strauss process
## interaction distance: 10
## Fitted interaction parameter gamma: 1.3448
##
## Relevant coefficients:
## Interaction
## 0.2962
##
## For standard errors, type coef(summary(x))
Tornadoes in Iowa.