Problem Set #4
Andy Green.
11/20/2020
setwd("C:/R")##Library
library(sf)## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tmap)
library(USAboundaries)## Warning: package 'USAboundaries' was built under R version 4.0.3library(ggplot2)
library(spatstat)## Warning: package 'spatstat' was built under R version 4.0.3## Loading required package: spatstat.data## Warning: package 'spatstat.data' was built under R version 4.0.3## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapse## Loading required package: rpart##
## spatstat 1.64-1 (nickname: 'Help you I can, yes!')
## For an introduction to spatstat, type 'beginner'##
## Note: spatstat version 1.64-1 is out of date by more than 6 months; we recommend upgrading to the latest version.library(hexbin)## Warning: package 'hexbin' was built under R version 4.0.3library(maptools)## Warning: package 'maptools' was built under R version 4.0.3## Loading required package: sp## Checking rgeos availability: TRUElibrary(USAboundariesData)Problem Set #4
- Consider tornado reports across Iowa over the period 1950–2017. Subset the reports by the state border and use the border as the analysis domain (window). Use Albers equal area conic projection centered on Iowa (EPSG:26975). Rescale the units to kilometers. Estimate the annual tornado rate (intensity) per 100 square kilometers. (20)
if (!file.exists("1950-2018-torn-initpoint")) {
download.file("http://www.spc.noaa.gov/gis/svrgis/zipped/1950-2018-torn-initpoint.zip",
destfile = "temporary.zip")
unzip("temporary.zip")
}sfdfP <- st_read(dsn = "1950-2018-torn-initpoint") %>%
st_transform(crs = 26975) %>% ### Albers Equal Area Conic projection, Iowa
filter(st == "IA", mag >= 0) %>%
mutate(EF = as.factor(mag)) %>%
select(EF)## Reading layer `1950-2018-torn-initpoint' from data source `C:\R\1950-2018-torn-initpoint' using driver `ESRI Shapefile'
## Simple feature collection with 63645 features and 22 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -163.53 ymin: 18.13 xmax: -64.9 ymax: 61.02
## geographic CRS: WGS 84IA.sf <- us_states(states = "IOWA")
tm_shape(IA.sf) +
tm_borders(col = "black") +
tm_shape(sfdfP) +
tm_bubbles(size = .05,
col = "red",
alpha = .4)
IA.sf <- us_states(states = "Iowa") %>%
st_transform(crs = st_crs(sfdfP)$proj4string)
IA.sp <- as(IA.sf, "Spatial")
IA.win <- as(IA.sp, 'owin')spdfP <- as(sfdfP, "Spatial")
T.ppp <- as(spdfP["EF"], "ppp")
plot(T.ppp)
T.ppp <- T.ppp[IA.win]
plot(T.ppp)
T.ppp <- rescale(T.ppp,
s = 1000,
unitname = "km")
summary(T.ppp)## Marked planar point pattern: 2545 points
## Average intensity 0.01746782 points per square km
##
## *Pattern contains duplicated points*
##
## Coordinates are given to 5 decimal places
##
## Multitype:
## frequency proportion intensity
## 0 1138 0.447151300 0.0078107570
## 1 820 0.322200400 0.0056281380
## 2 435 0.170923400 0.0029856590
## 3 106 0.041650290 0.0007275398
## 4 40 0.015717090 0.0002745433
## 5 6 0.002357564 0.0000411815
##
## Window: polygonal boundary
## single connected closed polygon with 253 vertices
## enclosing rectangle: [1244.1494, 1778.2757] x [877.5527, 1226.8102] km
## (534.1 x 349.3 km)
## Window area = 145696 square km
## Unit of length: 1 km
## Fraction of frame area: 0.781- Create a spatial density map with the events plotted as points. (10)
plot(density(T.ppp))
plot(T.ppp, add = TRUE)
- Create a
pppobject using the locations of cities and towns that have a population exceeding 500. Make a map showing the locations (20)
unzip("ci08au12.zip")
C.sf <- st_read(dsn = ".",
layer = "ci08au12",
quiet = TRUE) %>%
st_set_crs(4326) %>% ## mercator
st_transform(crs = st_crs(sfdfP)) %>%
filter(POP_1990 > 500)C.spdf <- as(C.sf, "Spatial")
C.ppp = as(C.spdf, "ppp")
C.ppp = unmark(C.ppp[IA.win])
C.ppp = rescale(C.ppp, s = 1000, unitname = "km")
plot(C.ppp)
Zc <- distmap(C.ppp)
plot(Zc)
- Plot the spatial density of tornado reports as a function of distance from nearest town. Use an adjustment to the bandwidth of 1.5. Make a plot. (10)
Y <- unmark(T.ppp)
summary(Y)## Planar point pattern: 2545 points
## Average intensity 0.01746782 points per square km
##
## *Pattern contains duplicated points*
##
## Coordinates are given to 5 decimal places
##
## Window: polygonal boundary
## single connected closed polygon with 253 vertices
## enclosing rectangle: [1244.1494, 1778.2757] x [877.5527, 1226.8102] km
## (534.1 x 349.3 km)
## Window area = 145696 square km
## Unit of length: 1 km
## Fraction of frame area: 0.781plot(Y)
X <- rpoispp(lambda = summary(Y)$intensity,
win = window(Y))
summary(X)$intensity## [1] 0.01653437plot(X)
Z <- superimpose(Y = Y, X = X)## Warning: data contain duplicated pointsplot(density(split(Z)))
rhat <- rhohat(Y,
covariate = Zc,
adjust = 1.5,
method = "transform")##Thank you for your willingness to help!K.df <- as.data.frame(Kest(T.ppp))K.df <- as.data.frame(Kest(T.ppp))
ggplot(K.df, aes(x = r, y = iso)) +
geom_line() +
geom_line(aes(y = theo), color = "blue") +
geom_vline(xintercept = 60, lty = 'dashed') +
geom_hline(yintercept = 12500, lty = 'dashed') +
geom_hline(yintercept = 11000, lty = 'dashed') +
xlab("Distance (km)") +
ylab("Tornadoes Near Towns")
- Plot the nearest neighbor distance function and test for clustering using a simulation envelope. Use 99 simulations. Describe the results. (10)
G.df <- as.data.frame(Gest(T.ppp)) %>%
filter(r < 8)
ggplot(G.df, aes(x = r, y = km)) +
geom_line() +
geom_line(aes(y = theo), color = "blue") +
geom_hline(yintercept = .4, lty = 'dashed') +
geom_vline(xintercept = c(3.2, 4), lty = 'dashed') +
xlab("Distance in Meters") +
ylab("Percent Distances within a distance of different tornado")
ST.ppp <- unmark(T.ppp[T.ppp$marks == 2 |
T.ppp$marks == 3 |
T.ppp$marks == 4 |
T.ppp$marks == 5])
env <- envelope(ST.ppp,
fun = Kest,
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.env.df <- as.data.frame(env)
ggplot(env.df, aes(x = r, y = obs)) +
geom_ribbon(aes(ymin = lo, ymax = hi), fill = "gray50") +
geom_line() +
geom_line(aes(y = theo), color = "blue", lty = 'dashed') +
xlab("Distance (km)") +
ylab("Possible number of additional tornadoes \nwithin a distance of another tornado")
- Create a multitype point pattern object by combining city locations and tornado touchdown locations. Note: remove the marks from the tornado point pattern data. Use “city” and “tornado” as factor marks. Plot the Ripley’s K cross function and test for ‘interaction’ between towns and tornado reports for distances less than 10 km. Use the
envelope()function to create 99 simulations. Interpret the results. (30)
T2.ppp <- unmark(T.ppp)
M.ppp <- superimpose(C.ppp, T2.ppp)## Warning: data contain duplicated pointsmarks(M.ppp) <- as.factor(c(rep("city", C.ppp$n),
rep("tornado", T2.ppp$n)))
plot(Kcross(M.ppp))
Kc= Kcross(M.ppp)modelM <- ppm(M.ppp, trend = ~ 1,
interaction = Strauss(r = 10),
rbord = 10)
modelM## Stationary Strauss process
##
## Possible marks: 'city' and 'tornado'
##
## First order terms:
## beta_city beta_tornado
## 0.006955742 0.006955742
##
## Interaction distance: 10
## Fitted interaction parameter gamma: 1.0641396
##
## Relevant coefficients:
## Interaction
## 0.06216661
##
## For standard errors, type coef(summary(x))
##
## *** Model is not valid ***
## *** Interaction parameters are outside valid range ***plot(envelope(modelM, Lest,
nsim = 99,
global = TRUE,
correction = 'border'),
legend = FALSE)## Generating 198 simulated realisations of fitted Gibbs model (99 to estimate the
## mean and 99 to calculate envelopes) ...
## 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.
##
## Done.
summary(modelM)## Point process model
## Fitting method: maximum pseudolikelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = M.ppp, trend = ~1, interaction = Strauss(r = 10),
## rbord = 10)
## Edge correction: "border"
## [border correction distance r = 10 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Marked planar point pattern: 3169 points
## Average intensity 0.0218 points per square km
## Multitype:
## frequency proportion intensity
## city 624 0.197 0.00428
## tornado 2540 0.803 0.01750
##
## Window: polygonal boundary
## single connected closed polygon with 253 vertices
## enclosing rectangle: [1244.1494, 1778.2757] x [877.5527, 1226.8102] km
## (534.1 x 349.3 km)
## Window area = 145696 square km
## Unit of length: 1 km
## Fraction of frame area: 0.781
##
## Dummy quadrature points:
## 120 x 120 grid of dummy points, plus 4 corner points
## dummy spacing: 4.451052 x 2.910479 km
##
## Original dummy parameters: =
## Marked planar point pattern: 25919 points
## Average intensity 0.178 points per square km
## Multitype:
## frequency proportion intensity
## city 13900 0.537 0.0955
## tornado 12000 0.463 0.0824
##
## Window: polygonal boundary
## single connected closed polygon with 253 vertices
## enclosing rectangle: [1244.1494, 1778.2757] x [877.5527, 1226.8102] km
## (534.1 x 349.3 km)
## Window area = 145696 square km
## Unit of length: 1 km
## Fraction of frame area: 0.781
## Quadrature weights:
## (counting weights based on 120 x 120 array of rectangular tiles)
## All weights:
## range: [0.559, 13] total: 291000
## Weights on data points:
## range: [0.559, 6.48] total: 17100
## Weights on dummy points:
## range: [0.559, 13] total: 274000
## --------------------------------------------------------------------------------
## FITTED MODEL:
##
## Stationary Strauss process
## Possible marks:
## city tornado
## ---- Trend: ----
##
##
## First order terms:
## beta_city beta_tornado
## 0.006955742 0.006955742
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -4.96818780 0.040564668 -5.04769309 -4.8886825 *** -122.47574
## Interaction 0.06216661 0.004712277 0.05293072 0.0714025 *** 13.19248
##
## ---- Interaction: -----
##
## Interaction: Strauss process
## Interaction distance: 10
## Fitted interaction parameter gamma: 1.0641396
##
## Relevant coefficients:
## Interaction
## 0.06216661
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept) Interaction
## -4.96818780 0.06216661
##
## Fitted exp(theta):
## (Intercept) Interaction
## 0.006955742 1.064139627
##
## *** Model is not valid ***
## *** Interaction parameters are outside valid range ***plot(modelM)
###Those in a city where a tornado has been observed are less likely to see another tornado.