Download the data and read it into R. Use the rgdal package for the readOGR() function. Print the projection information.
tmp = download.file("http://www.spc.noaa.gov/gis/svrgis/zipped/tornado.zip",
"tornado.zip", mode = "wb")
unzip("tornado.zip")
require(rgdal)
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
Subset on year and F scale and create a factor variable (Category) for the F-scale. Convert length from miles and the width from yards to meters. Compute path area.
Torn2 = subset(Torn, FSCALE > 0 & YEAR >= 1950)
Torn2$Category = factor(Torn2$FSCALE, ordered = TRUE)
Torn2$Length = Torn2$LENGTH * 1609.34
Torn2$Width = Torn2$WIDTH * 0.9144
Torn2$Area = Torn2$Width * Torn2$Length
Select domain: Florida (12). http://www2.census.gov/cgi-bin/shapefiles2009/state-files?state=12 http://www.census.gov/geo/www/cob/st2000.html#shp
tmp = download.file("http://www.census.gov/geo/cob/bdy/st/st00shp/st12_d00_shp.zip",
"st12_d00_shp.zip", mode = "wb")
unzip("st12_d00_shp.zip")
STATE = readOGR(dsn = ".", layer = "st12_d00")
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st12_d00"
## with 14 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
Add projection information and transform the state boundaries.
ll = "+proj=longlat +ellps=GRS80 +datum=NAD83"
proj4string(STATE) = CRS(ll)
STATE = spTransform(STATE, CRS(proj4string(Torn2)))
require(rgeos)
stArea = gArea(STATE)
Subset by state boundary.
st = overlay(Torn2, STATE)
st = !is.na(st)
Torn2 = Torn2[st, ]
Check the subsetting.
plot(STATE)
plot(Torn2, add = TRUE)
table(Torn2$STATE)
##
## AK AL AR AZ CA CO CT DC DE FL GA HI IA ID IL
## 0 0 0 0 0 0 0 0 0 1114 0 0 0 0 0
## IN KS KY LA MA MD ME MI MN MO MS MT NC ND NE
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## NH NJ NM NV NY OH OK OR PA PR RI SC SD TN TX
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## UT VA VT WA WI WV WY
## 0 0 0 0 0 0 0
Compute total area of tornado paths and annual average area.
tArea = sum(Torn2$Area)
areaPerYear = tArea/length(1950:2011)
stArea/areaPerYear # Annual avg recurrence interval for any location
## [1] 15658
areaPerYear/stArea * 100 # Annual probability of getting struck by a tornado
## [1] 0.006387
Next, create and summarize a planar point pattern object.
require(spatstat)
require(maptools)
W = as(STATE, "owin")
## Checking 14 polygons...1, [Checking polygon with 2918 edges...]
## 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14.
## done.
## Checking for cross-intersection between 14 polygons...1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13.
## done.
All.ppp = as.ppp(Torn2["Category"])
All.ppp = All.ppp[W]
unitname(All.ppp) = c("meter", "meters")
summary(All.ppp)
## Marked planar point pattern: 1114 points
## Average intensity 7.25e-09 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
## 1 784 0.704000 5.10e-09
## 2 296 0.266000 1.93e-09
## 3 33 0.029600 2.15e-10
## 4 1 0.000898 6.51e-12
## 5 0 0.000000 0.00e+00
##
## Window: polygonal boundary
## 14 separate polygons (no holes)
## vertices area relative.area
## polygon 1 2918 1.527e+11 9.95e-01
## polygon 2 36 4.759e+07 3.10e-04
## polygon 3 22 5.163e+07 3.36e-04
## polygon 4 26 1.412e+07 9.20e-05
## polygon 5 69 8.935e+07 5.82e-04
## polygon 6 289 1.714e+08 1.12e-03
## polygon 7 25 6.844e+06 4.46e-05
## polygon 8 15 1.089e+07 7.09e-05
## polygon 9 9 7.209e+05 4.69e-06
## polygon 10 106 3.606e+07 2.35e-04
## polygon 11 19 1.975e+07 1.29e-04
## polygon 12 179 3.577e+08 2.33e-03
## polygon 13 11 3.512e+06 2.29e-05
## polygon 14 12 2.759e+07 1.80e-04
## enclosing rectangle: [802437, 1617257] x [-1501045, -800485] meters
## Window area = 1.53583e+11 square meters
## Unit of length: 1 meter
The native spatial units are meters.
Duplicated events:
sum(duplicated(All.ppp))
## [1] 128
Create two ppp objects and join them. This is needed for the relative risk function.
T = All.ppp
plot(split(T))
H = unmark(T[T$marks <= 2 & T$marks >= 1])
I = unmark(T[T$marks >= 3])
T2 = superimpose(H = H, I = I)
## Warning: data contain duplicated points
summary(T2)
## Marked planar point pattern: 1114 points
## Average intensity 7.25e-09 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
## H 1080 0.9690 7.03e-09
## I 34 0.0305 2.21e-10
##
## Window: polygonal boundary
## 14 separate polygons (no holes)
## vertices area relative.area
## polygon 1 2918 1.527e+11 9.95e-01
## polygon 2 36 4.759e+07 3.10e-04
## polygon 3 22 5.163e+07 3.36e-04
## polygon 4 26 1.412e+07 9.20e-05
## polygon 5 69 8.935e+07 5.82e-04
## polygon 6 289 1.714e+08 1.12e-03
## polygon 7 25 6.844e+06 4.46e-05
## polygon 8 15 1.089e+07 7.09e-05
## polygon 9 9 7.209e+05 4.69e-06
## polygon 10 106 3.606e+07 2.35e-04
## polygon 11 19 1.975e+07 1.29e-04
## polygon 12 179 3.577e+08 2.33e-03
## polygon 13 11 3.512e+06 2.29e-05
## polygon 14 12 2.759e+07 1.80e-04
## enclosing rectangle: [802437, 1617257] x [-1501045, -800485] meters
## Window area = 1.53583e+11 square meters
## Unit of length: 1 meter
Plot the intensy events.
plot(I, pch = 25, cols = "red", main = "")
Compute spatial densities and relative risk. calculate bw first, then do densities/rr.
eps = 500
bw = bw.relrisk(T2, hmin = 500, hmax = 1e+10)
dAll = density(unmark(T2), dimyx = c(eps, eps), sigma = bw)
mean(dAll$v, na.rm = TRUE)
## [1] 7.246e-09
dH = density(H, dimyx = c(eps, eps), signma = bw)
dI = density(I, dimyx = c(eps, eps), sigma = bw)
rr = relrisk(T2, dimyx = c(eps, eps), sigma = bw)
Statewide mean relative risk.
mean(rr$v, na.rm = TRUE)
## [1] 0.03039
sd(as.vector(rr$v), na.rm = TRUE)
## [1] 1.413e-05
range(rr$v, na.rm = TRUE)
## [1] 0.03032 0.03042
Map them on a 4-panel plot.
par(mfrow = c(2, 2))
plot(dAll, col = terrain.colors(10), ribbon = FALSE, main = "")
contour(dAll, add = TRUE)
## NULL
plot(dH, col = terrain.colors(10), ribbon = FALSE, main = "")
contour(dH, add = TRUE)
## NULL
plot(dI, col = terrain.colors(10), ribbon = FALSE, main = "")
contour(dI, add = TRUE)
## NULL
plot(rr, ribbon = FALSE, main = "")
contour(rr, add = TRUE)
## NULL
Create a SpatialGridDataFrame (S4 class).
rr.sgdf = as(rr, "SpatialGridDataFrame")
proj4string(rr.sgdf) = CRS(proj4string(Torn2))
Map the relative risk.
require(maps)
## Loading required package: maps
brd = map("county", "florida", plot = FALSE)
brd_ll = map2SpatialLines(brd, proj4string = CRS(ll))
brd_lcc = spTransform(brd_ll, CRS(proj4string(Torn2)))
require(RColorBrewer)
## Loading required package: RColorBrewer
rr.sgdf$v = rr.sgdf$v * 100
range(rr.sgdf$v, na.rm = TRUE)
## [1] 3.032 3.042
rng = seq(3.032, 3.044, 0.002)
cls = brewer.pal(6, "YlOrRd")
l3 = list("sp.lines", brd_lcc, lwd = 0.5, first = FALSE)
p1 = spplot(rr.sgdf, "v", col.regions = cls, at = rng, sp.layout = l3, colorkey = list(space = "bottom",
labels = list(cex = 1.2)), sub = list("Relative probability (as %) that a tornado will be intense (EF3 or stronger)",
cex = 1.2, font = 1))
p1
Download the county boundaries (all US).
tmp = download.file("http://www.nws.noaa.gov/geodata/catalog/county/data/c_02ap13.zip",
"c_02ap13.zip", mode = "wb")
unzip("c_02ap13.zip")
COUNTY = readOGR(dsn = ".", layer = "c_02ap13")
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "c_02ap13"
## with 3312 features and 11 fields
## Feature type: wkbPolygon with 2 dimensions
COUNTY = COUNTY[COUNTY$STATE == "FL", ]
COUNTY = spTransform(COUNTY, CRS(proj4string(Torn2)))
Extract the mean relative risk per county.
COUNTY$v = over(COUNTY, rr.sgdf, fn = mean)
spplot(COUNTY, "v")
county.df = as.data.frame(COUNTY)
county.df[rev(order(county.df$v))[1:6], ]
## OBJECTID_1 STATE CWA COUNTYNAME FIPS TIME_ZONE FE_AREA LON LAT
## 3168 3180 FL JAX Nassau 12089 E ne -81.81 30.61
## 3167 3179 FL JAX Duval 12031 E ne -81.67 30.33
## 2478 2483 FL JAX St. Johns 12109 E ne -81.44 29.90
## 2576 2581 FL JAX Flagler 12035 E ne -81.31 29.46
## 3212 3224 FL MLB Brevard 12009 E <NA> -80.74 28.26
## 3213 3225 FL MLB Volusia 12127 E ec -81.19 29.06
## Shape_Leng Shape_Area v
## 3168 4.374 0.1597 3.042
## 3167 4.051 0.2033 3.041
## 2478 5.149 0.1589 3.041
## 2576 3.252 0.1213 3.041
## 3212 9.249 0.2496 3.041
## 3213 8.093 0.2965 3.041
Subset by county boundary.
CTY = subset(COUNTY, COUNTY$COUNTYNAME == "Escambia")
ct = overlay(Torn2, CTY)
ct = !is.na(ct)
Torn3 = Torn2[ct, ]
table(Torn3$Category)
##
## 1 2 3 4 5
## 20 7 2 0 0
Convert coordinates to units of kilometers.
# T = affine(All.ppp, mat = diag(c(1/1, 1/1))) unitname(T) = c('km',
# 'kms')
summary(T)
Get data on the location of cities.
tmp = download.file("http://www.nws.noaa.gov/geodata/catalog/national/data/ci08au12.zip",
"ci08au12.zip", mode = "wb")
## Error: cannot open destfile 'ci08au12.zip', reason 'Permission denied'
unzip("ci08au12.zip")
US = readOGR(dsn = ".", layer = "ci08au12", p4s = ll)
## 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(Torn2)))
st = overlay(UST, STATE)
st = !is.na(st)
Cities = UST[st, ]
plot(unmark(T), pch = 16, cex = 0.4)
plot(Cities, pch = 15, col = "green", cex = 0.4, add = TRUE)
C = as.ppp(Cities)
window = as.owin(T)
C = C[window]
Create a map of distance from nearest city centers. Distance from nearest city is a covariate in the point pattern model (ppm).
Zc = distmap(C)
plot(Zc)
Examine the smoothed estimate of spatial density as a function of the covariate.
H = unmark(T[T$marks <= 2 & T$marks >= 1])
I = unmark(T[T$marks >= 3])
plot(rhohat(H, Zc, smoother = "kernel", method = "transform"), legend = FALSE,
main = "")
plot(rhohat(I, Zc, smoother = "kernel", method = "transform"), legend = FALSE,
main = "")
Spatial density as function of distance to nearest city
rhat = rhohat(unmark(T), Zc, smoother = "kernel", method = "transform") # max report density as fcn of distance
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## Warning: 1 out of 192 pixel values were outside the pixel image domain and
## were estimated by convolution
rptsPerYear = max(rhat$rho)/length(1950:2011) * stArea # Reports per year
areaPerYear = rptsPerYear * mean(Torn2$Area)
stArea/areaPerYear # Annual avg recurrence interval for any location
## [1] 6125
areaPerYear/stArea * 100 # Annual probability of a tornado at any location
## [1] 0.01633
Annual average recurrence interval for any location (bias corrected): 6125 Annual probability of a tornado at any location: 0.01163
Model the spatial density of the intense tornadoes using distance from city center as a covariate.
fit1 = ppm(I, ~Zc, covariates = list(Zc = Zc), correction = "border", rbord = 5)
summary(fit1)
## Point process model
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm(Q = I, trend = ~Zc, covariates = list(Zc = Zc), correction = "border",
## rbord = 5)
## Edge correction: "border"
## [border correction distance r = 5 ]
##
## ----------------------------------------------------
## Quadrature scheme = data + dummy + weights
## Data pattern:
## Planar point pattern: 34 points
## Average intensity 2.21e-10 points per square meter
## Window: polygonal boundary
## 14 separate polygons (no holes)
## vertices area relative.area
## polygon 1 2918 1.527e+11 9.95e-01
## polygon 2 36 4.759e+07 3.10e-04
## polygon 3 22 5.163e+07 3.36e-04
## polygon 4 26 1.412e+07 9.20e-05
## polygon 5 69 8.935e+07 5.82e-04
## polygon 6 289 1.714e+08 1.12e-03
## polygon 7 25 6.844e+06 4.46e-05
## polygon 8 15 1.089e+07 7.09e-05
## polygon 9 9 7.209e+05 4.69e-06
## polygon 10 106 3.606e+07 2.35e-04
## polygon 11 19 1.975e+07 1.29e-04
## polygon 12 179 3.577e+08 2.33e-03
## polygon 13 11 3.512e+06 2.29e-05
## polygon 14 12 2.759e+07 1.80e-04
## enclosing rectangle: [802437, 1617257] x [-1501045, -800485] meters
## Window area = 1.53583e+11 square meters
## Unit of length: 1 meter
##
##
## Dummy quadrature points:
## ( 32 x 32 grid, plus 4 corner points)
## Planar point pattern: 336 points
## Average intensity 2.19e-09 points per square meter
## Window: polygonal boundary
## 14 separate polygons (no holes)
## vertices area relative.area
## polygon 1 2918 1.527e+11 9.95e-01
## polygon 2 36 4.759e+07 3.10e-04
## polygon 3 22 5.163e+07 3.36e-04
## polygon 4 26 1.412e+07 9.20e-05
## polygon 5 69 8.935e+07 5.82e-04
## polygon 6 289 1.714e+08 1.12e-03
## polygon 7 25 6.844e+06 4.46e-05
## polygon 8 15 1.089e+07 7.09e-05
## polygon 9 9 7.209e+05 4.69e-06
## polygon 10 106 3.606e+07 2.35e-04
## polygon 11 19 1.975e+07 1.29e-04
## polygon 12 179 3.577e+08 2.33e-03
## polygon 13 11 3.512e+06 2.29e-05
## polygon 14 12 2.759e+07 1.80e-04
## enclosing rectangle: [802437, 1617257] x [-1501045, -800485] meters
## Window area = 1.53583e+11 square meters
## Unit of length: 1 meter
##
##
## Quadrature weights:
## (counting weights based on 32 x 32 array of rectangular tiles)
## All weights:
## range: [34800000, 5.57e+08] total: 1.54e+11
## Weights on data points:
## range: [34800000, 2.79e+08] total: 7.75e+09
## Weights on dummy points:
## range: [34800000, 5.57e+08] total: 1.46e+11
##
## ----------------------------------------------------
## FITTED MODEL:
##
## Nonstationary Poisson process
##
##
## ---- Intensity: ----
##
## Trend formula: ~Zc
## Model depends on external covariate 'Zc'
##
## Covariates provided:
## Zc: im
##
## Fitted coefficients for trend formula:
## (Intercept) Zc
## -2.126e+01 -2.016e-04
##
## Estimate S.E. Ztest CI95.lo CI95.hi
## (Intercept) -2.126e+01 0.3022897 na -2.185e+01 -2.067e+01
## Zc -2.016e-04 0.0000643 ** -3.276e-04 -7.562e-05
##
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept) Zc
## -2.126e+01 -2.016e-04
##
## Fitted exp(theta):
## (Intercept) Zc
## 5.852e-10 9.998e-01
Create another covariate (distance from nearest weak tornado).
Zh = distmap(H)
Model the spatial density using two covariates.
fit2 = ppm(I, ~Zc + Zh, covariates = list(Zc = Zc, Zh = Zh), correction = "border",
rbord = 5)
## Warning: glm.fit: algorithm did not converge
summary(fit2)
## Point process model
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## *** Algorithm did not converge ***
## Call:
## ppm(Q = I, trend = ~Zc + Zh, covariates = list(Zc = Zc, Zh = Zh),
## correction = "border", rbord = 5)
## Edge correction: "border"
## [border correction distance r = 5 ]
##
## ----------------------------------------------------
## Quadrature scheme = data + dummy + weights
## Data pattern:
## Planar point pattern: 34 points
## Average intensity 2.21e-10 points per square meter
## Window: polygonal boundary
## 14 separate polygons (no holes)
## vertices area relative.area
## polygon 1 2918 1.527e+11 9.95e-01
## polygon 2 36 4.759e+07 3.10e-04
## polygon 3 22 5.163e+07 3.36e-04
## polygon 4 26 1.412e+07 9.20e-05
## polygon 5 69 8.935e+07 5.82e-04
## polygon 6 289 1.714e+08 1.12e-03
## polygon 7 25 6.844e+06 4.46e-05
## polygon 8 15 1.089e+07 7.09e-05
## polygon 9 9 7.209e+05 4.69e-06
## polygon 10 106 3.606e+07 2.35e-04
## polygon 11 19 1.975e+07 1.29e-04
## polygon 12 179 3.577e+08 2.33e-03
## polygon 13 11 3.512e+06 2.29e-05
## polygon 14 12 2.759e+07 1.80e-04
## enclosing rectangle: [802437, 1617257] x [-1501045, -800485] meters
## Window area = 1.53583e+11 square meters
## Unit of length: 1 meter
##
##
## Dummy quadrature points:
## ( 32 x 32 grid, plus 4 corner points)
## Planar point pattern: 336 points
## Average intensity 2.19e-09 points per square meter
## Window: polygonal boundary
## 14 separate polygons (no holes)
## vertices area relative.area
## polygon 1 2918 1.527e+11 9.95e-01
## polygon 2 36 4.759e+07 3.10e-04
## polygon 3 22 5.163e+07 3.36e-04
## polygon 4 26 1.412e+07 9.20e-05
## polygon 5 69 8.935e+07 5.82e-04
## polygon 6 289 1.714e+08 1.12e-03
## polygon 7 25 6.844e+06 4.46e-05
## polygon 8 15 1.089e+07 7.09e-05
## polygon 9 9 7.209e+05 4.69e-06
## polygon 10 106 3.606e+07 2.35e-04
## polygon 11 19 1.975e+07 1.29e-04
## polygon 12 179 3.577e+08 2.33e-03
## polygon 13 11 3.512e+06 2.29e-05
## polygon 14 12 2.759e+07 1.80e-04
## enclosing rectangle: [802437, 1617257] x [-1501045, -800485] meters
## Window area = 1.53583e+11 square meters
## Unit of length: 1 meter
##
##
## Quadrature weights:
## (counting weights based on 32 x 32 array of rectangular tiles)
## All weights:
## range: [34800000, 5.57e+08] total: 1.54e+11
## Weights on data points:
## range: [34800000, 2.79e+08] total: 7.75e+09
## Weights on dummy points:
## range: [34800000, 5.57e+08] total: 1.46e+11
##
## ----------------------------------------------------
## FITTED MODEL:
##
## Nonstationary Poisson process
##
##
## ---- Intensity: ----
##
## Trend formula: ~Zc + Zh
## Model depends on external covariates 'Zc' and 'Zh'
##
## Covariates provided:
## Zc: im
## Zh: im
##
## Fitted coefficients for trend formula:
## (Intercept) Zc Zh
## -2.046e+01 -6.061e-05 -2.429e-04
##
## Estimate S.E. Ztest CI95.lo CI95.hi
## (Intercept) -2.046e+01 3.285e-01 na -2.111e+01 -1.982e+01
## Zc -6.061e-05 6.993e-05 -1.977e-04 7.645e-05
## Zh -2.429e-04 6.201e-05 *** -3.644e-04 -1.213e-04
##
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept) Zc Zh
## -2.046e+01 -6.061e-05 -2.429e-04
##
## Fitted exp(theta):
## (Intercept) Zc Zh
## 1.296e-09 9.999e-01 9.998e-01
Can you drop one of the covariates?
fit3 = ppm(I, ~Zh, covariates = list(Zh = Zh), correction = "border", rbord = 5)
## Warning: glm.fit: algorithm did not converge
AIC(fit1)
## [1] 1571
AIC(fit2)
## [1] 1549
AIC(fit3)
## [1] 1548
2 * (as.numeric(logLik(fit2)) - as.numeric(logLik(fit3)))/I$n
## [1] 0.0234
Does adding a term for clustering improve things?
fit4 = kppm(I, ~Zc + Zh, covariates = list(Zc = Zc, Zh = Zh), interaction = "Thomas",
correction = "border", rbord = 5)
## Warning: glm.fit: algorithm did not converge
fit4
## Inhomogeneous cluster point process model
## Fitted to point pattern dataset 'I'
## Fitted by minimum contrast
## Summary statistic: inhomogeneous K-function
##
## Trend formula: ~Zc + Zh
##
## Fitted coefficients for trend formula:
## (Intercept) Zc Zh
## -2.046e+01 -6.061e-05 -2.429e-04
##
## Cluster model: Thomas process
## Fitted cluster parameters:
## kappa sigma2
## 1.621e-10 1.575e+09
## Mean cluster size: [pixel image]
Prediction correlation.
M = quadrat.test(fit2, nx = 8, ny = 8)
## Warning: Some expected counts are small; chi^2 approximation may be
## inaccurate
cor.test(M$expected, M$observed)
##
## Pearson's product-moment correlation
##
## data: M$expected and M$observed
## t = 3.913, df = 30, p-value = 0.0004849
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2917 0.7733
## sample estimates:
## cor
## 0.5813
There are only 37 “I” events when F3 is used. Here I try again by defining F2 and higher as “I” events. This time there are 330 “I” events.
Create two ppp objects and join them. This is needed for the relative risk function.
H2 = unmark(T[T$marks == 1])
I2 = unmark(T[T$marks >= 2])
T3 = superimpose(H2 = H2, I2 = I2)
## Warning: data contain duplicated points
summary(T3)
## Marked planar point pattern: 1114 points
## Average intensity 7.25e-09 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
## H2 784 0.704 5.10e-09
## I2 330 0.296 2.15e-09
##
## Window: polygonal boundary
## 14 separate polygons (no holes)
## vertices area relative.area
## polygon 1 2918 1.527e+11 9.95e-01
## polygon 2 36 4.759e+07 3.10e-04
## polygon 3 22 5.163e+07 3.36e-04
## polygon 4 26 1.412e+07 9.20e-05
## polygon 5 69 8.935e+07 5.82e-04
## polygon 6 289 1.714e+08 1.12e-03
## polygon 7 25 6.844e+06 4.46e-05
## polygon 8 15 1.089e+07 7.09e-05
## polygon 9 9 7.209e+05 4.69e-06
## polygon 10 106 3.606e+07 2.35e-04
## polygon 11 19 1.975e+07 1.29e-04
## polygon 12 179 3.577e+08 2.33e-03
## polygon 13 11 3.512e+06 2.29e-05
## polygon 14 12 2.759e+07 1.80e-04
## enclosing rectangle: [802437, 1617257] x [-1501045, -800485] meters
## Window area = 1.53583e+11 square meters
## Unit of length: 1 meter
Plot the intensy events.
plot(I2, pch = 25, cols = "red", main = "")
Compute spatial densities and relative risk. calculate bw first, then do densities/rr.
eps = 500
bw = bw.relrisk(T3, hmax = 1e+10)
dAll2 = density(unmark(T3), dimyx = c(eps, eps), sigma = bw)
mean(dAll2$v, na.rm = TRUE)
## [1] 7.239e-09
dH2 = density(H2, dimyx = c(eps, eps), signma = bw)
dI2 = density(I2, dimyx = c(eps, eps), sigma = bw)
rr2 = relrisk(T3, dimyx = c(eps, eps), sigma = bw)
Statewide mean relative risk.
mean(rr2$v, na.rm = TRUE)
## [1] 0.2955
sd(as.vector(rr2$v), na.rm = TRUE)
## [1] 0.002866
range(rr2$v, na.rm = TRUE)
## [1] 0.2898 0.3013
Map them on a 4-panel plot.
par(mfrow = c(2, 2))
plot(dAll2, col = terrain.colors(10), ribbon = FALSE, main = "")
contour(dAll2, add = TRUE)
## NULL
plot(dH2, col = terrain.colors(10), ribbon = FALSE, main = "")
contour(dH2, add = TRUE)
## NULL
plot(dI2, col = terrain.colors(10), ribbon = FALSE, main = "")
contour(dI2, add = TRUE)
## NULL
plot(rr2, ribbon = FALSE, main = "")
contour(rr2, add = TRUE)
## NULL
Create a SpatialGridDataFrame (S4 class).
rr2.sgdf = as(rr2, "SpatialGridDataFrame")
proj4string(rr2.sgdf) = CRS(proj4string(Torn2))
Map the relative risk.
require(maps)
brd = map("county", "florida", plot = FALSE)
brd_ll = map2SpatialLines(brd, proj4string = CRS(ll))
brd_lcc = spTransform(brd_ll, CRS(proj4string(Torn2)))
require(RColorBrewer)
rr2.sgdf$v = rr2.sgdf$v * 100
range(rr2.sgdf$v, na.rm = TRUE)
## [1] 28.98 30.13
rng = seq(28.8, 30.2, 0.2)
cls = brewer.pal(7, "YlOrRd")
l3 = list("sp.lines", brd_lcc, lwd = 0.5, first = FALSE)
p1 = spplot(rr2.sgdf, "v", col.regions = cls, at = rng, sp.layout = l3, colorkey = list(space = "bottom",
labels = list(cex = 1.2)), sub = list("Relative probability (as %) that a tornado will be intense (EF2 or stronger)",
cex = 1.2, font = 1))
p1
The relative risk flips from the previous version.