Untitled
Gregory Burris
Thursday, April 02, 2015
Load up some libraries.
library(maptools)## Loading required package: sp
## Checking rgeos availability: FALSE
## Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
## which has a restricted licence. It is disabled by default;
## to enable gpclib, type gpclibPermit()library(rgdal)## rgdal: version: 0.9-2, (SVN revision 526)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
## Path to GDAL shared files: C:/Users/Franklin/Documents/R/win-library/3.2/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: C:/Users/Franklin/Documents/R/win-library/3.2/rgdal/projlibrary(spatstat)##
## spatstat 1.41-1 (nickname: 'Ides of March')
## For an introduction to spatstat, type 'beginner'setwd("C:/Project3/point_data")
USshape = readOGR(".","US_simple_dissolve")## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "US_simple_dissolve"
## with 1 features
## It has 1 fieldsplot(USshape)
USshape2 = as(USshape, "SpatialPolygons")
USshape3 = as(USshape2, "owin")
stations = readOGR(".","stations")## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "stations"
## with 1920 features
## It has 11 fieldsstations.sp = as(stations, "SpatialPoints")
stations.ppp = as(stations.sp, 'ppp')
stations.ppp$window = USshape3
plot(stations.ppp)
```
Setup the boundary. Recall that point pattern statistics are sensitive to the boundary chosen. In this case we want the boundary to be the contiguous United States. As seen above, the R default is a bounding box around the points. My version of Arc is won’t let me simplify the US polygon. If you are able to do it, please share it with the rest of us.
K function (this will be slow)
K_stations = Kest(stations.ppp, correction="Ripley")
plot(K_stations)
## lty col key label
## iso 1 1 iso italic(hat(K)[iso](r))
## theo 2 2 theo italic(K[pois](r))
## meaning
## iso Ripley isotropic correction estimate of K(r)
## theo theoretical Poisson K(r)L function for all the data (using the “translation” edge correction to speed things up).
L_stations_9 <- envelope(stations.ppp, Lest, correction="translation", nsim = 9, nrank = 1, verbose=TRUE)## Generating 9 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9.
##
## Done.L_stations_9## Pointwise critical envelopes for L(r)
## and observed value for 'stations.ppp'
## Edge correction: "trans"
## Obtained from 9 simulations of CSR
## Alternative: two.sided
## Significance level of pointwise Monte Carlo test: 2/10 = 0.2
## .....................................................................
## Math.label Description
## r r distance argument r
## obs hat(L)[obs](r) observed value of L(r) for data pattern
## theo L[theo](r) theoretical value of L(r) for CSR
## lo hat(L)[lo](r) lower pointwise envelope of L(r) from simulations
## hi hat(L)[hi](r) upper pointwise envelope of L(r) from simulations
## .....................................................................
## Default plot formula: .~r
## where "." stands for 'obs', 'theo', 'hi', 'lo'
## Columns 'lo' and 'hi' will be plotted as shading (by default)
## Recommended range of argument r: [0, 729590]
## Available range of argument r: [0, 729590]plot(L_stations_9, .-r~r) # the second argument ensures that CSR is a horizontal line at zero
## lty col key label
## obs 1 1 obs italic(hat(L)[obs](r)-r)
## theo 2 2 theo italic(L[theo](r)-r)
## hi 1 8 hi italic(hat(L)[hi](r)-r)
## lo 1 8 lo italic(hat(L)[lo](r)-r)
## meaning
## obs observed value of L(r) for data pattern
## theo theoretical value of L(r) for CSR
## hi upper pointwise envelope of L(r) from simulations
## lo lower pointwise envelope of L(r) from simulationsFull L function
alpha = 0.05
sims = 999
rank = trunc((alpha*(sims+1))/2) # nrank must be an integer; use trunc (truncate) instead of round to be conservative
L_stations_999 = envelope(stations.ppp, Lest, correction="translation", nsim=sims, nrank=rank, verbose=TRUE)## Generating 999 simulations of CSR ...
## 1, 2, 3, ......10 [etd 12:57] .........20 [etd 12:33] .........30 [etd 12:22] .........40 [etd 12:12] .........50
## [etd 11:59] .........60 [etd 11:52] .........70 [etd 11:42] .........80 [etd 11:31] .........90 [etd 11:24] .........100
## [etd 11:17] .........110 [etd 11:13] .........120 [etd 11:05] .........130 [etd 10:58] .........140 [etd 10:52] .........150
## [etd 10:45] .........160 [etd 10:37] .........170 [etd 10:31] .........180 [etd 10:24] .........190 [etd 10:16] .........200
## [etd 10:09] .........210 [etd 10:00] .........220 [etd 9:52] .........230 [etd 9:44] .........240 [etd 9:38] .........250
## [etd 9:33] .........260 [etd 9:25] .........270 [etd 9:20] .........280 [etd 9:12] .........290 [etd 9:05] .........300
## [etd 8:57] .........310 [etd 8:49] .........320 [etd 8:42] .........330 [etd 8:34] .........340 [etd 8:26] .........350
## [etd 8:18] .........360 [etd 8:10] .........370 [etd 8:03] .........380 [etd 7:55] .........390 [etd 7:47] .........400
## [etd 7:40] .........410 [etd 7:32] .........420 [etd 7:24] .........430 [etd 7:17] .........440 [etd 7:09] .........450
## [etd 7:02] .........460 [etd 6:54] .........470 [etd 6:46] .........480 [etd 6:38] .........490 [etd 6:30] .........500
## [etd 6:23] .........510 [etd 6:15] .........520 [etd 6:07] .........530 [etd 6:00] .........540 [etd 5:52] .........550
## [etd 5:44] .........560 [etd 5:37] .........570 [etd 5:29] .........580 [etd 5:21] .........590 [etd 5:13] .........600
## [etd 5:06] .........610 [etd 4:58] .........620 [etd 4:50] .........630 [etd 4:42] .........640 [etd 4:34] .........650
## [etd 4:27] .........660 [etd 4:19] .........670 [etd 4:12] .........680 [etd 4:04] .........690 [etd 3:57] .........700
## [etd 3:49] .........710 [etd 3:41] .........720 [etd 3:34] .........730 [etd 3:26] .........740 [etd 3:18] .........750
## [etd 3:11] .........760 [etd 3:03] .........770 [etd 2:55] .........780 [etd 2:48] .........790 [etd 2:40] .........800
## [etd 2:32] .........810 [etd 2:25] .........820 [etd 2:17] .........830 [etd 2:09] .........840 [etd 2:02] .........850
## [etd 1:54] .........860 [etd 1:46] .........870 [etd 1:39] .........880 [etd 1:31] .........890 [etd 1:23] .........900
## [etd 1:16] .........910 [etd 1:08] .........920 [etd 1:00] .........930 [etd 53 sec] .........940 [etd 45 sec] .........950
## [etd 38 sec] .........960 [etd 30 sec] .........970 [etd 22 sec] .........980 [etd 15 sec] .........990 [etd 7 sec] ........ 999.
##
## Done.L_stations_999 # note that the reported "Significance level" may not match alpha due to truncation of nrank## Pointwise critical envelopes for L(r)
## and observed value for 'stations.ppp'
## Edge correction: "trans"
## Obtained from 999 simulations of CSR
## Alternative: two.sided
## Significance level of pointwise Monte Carlo test: 50/1000 = 0.05
## .....................................................................
## Math.label Description
## r r distance argument r
## obs hat(L)[obs](r) observed value of L(r) for data pattern
## theo L[theo](r) theoretical value of L(r) for CSR
## lo hat(L)[lo](r) lower pointwise envelope of L(r) from simulations
## hi hat(L)[hi](r) upper pointwise envelope of L(r) from simulations
## .....................................................................
## Default plot formula: .~r
## where "." stands for 'obs', 'theo', 'hi', 'lo'
## Columns 'lo' and 'hi' will be plotted as shading (by default)
## Recommended range of argument r: [0, 729590]
## Available range of argument r: [0, 729590]plot(L_stations_999, .-r~r)
## lty col key label
## obs 1 1 obs italic(hat(L)[obs](r)-r)
## theo 2 2 theo italic(L[theo](r)-r)
## hi 1 8 hi italic(hat(L)[hi](r)-r)
## lo 1 8 lo italic(hat(L)[lo](r)-r)
## meaning
## obs observed value of L(r) for data pattern
## theo theoretical value of L(r) for CSR
## hi upper pointwise envelope of L(r) from simulations
## lo lower pointwise envelope of L(r) from simulationsCompare the weather stations and USCSR
# extract the two types of point data
stations_ws_shp = stations[stations$type %in% c('ws','USCR'),]
# the cross-K function expects a factor to mark the different types of points
stations_ws_shp$mark = factor(stations_ws_shp$type)
# print out some descriptives
summary(stations_ws_shp$mark)## USCR ws
## 116 1804Convert the data into a point pattern
stations_ws_sp = as(stations_ws_shp, "SpatialPoints")
stations_ws = as(stations_ws_sp, 'ppp')
stations_marked_ws = ppp(stations_ws$x, stations_ws$y,
marks=stations_ws_shp$mark,
window=USshape3)## Warning in ppp(stations_ws$x, stations_ws$y, marks = stations_ws_shp$mark,
## : data contain duplicated pointsplot(stations_marked_ws)
Compute the K-function and L-function. This requires a “marked” point pattern, which was created above (robberies_marked_tw). We need to communicate more information to the envelope function than earlier. The parameters i and j relate to the cases and controls we discussed in class. The simulate parameter communicates HOW the significance testing should be done; similar to Goreaud and Pelissier (2003) we consider two cases: random labeling (rlabel) and random shift (rshift).
# random shift
L_stations_ws_rs99 = envelope(stations_marked_ws, Lcross,
i="USCR", j="ws",
correction="translation", nsim=sims, nrank=rank,
simulate=expression(rshift(stations_marked_ws,
which='ws', radius=2000,
edge='none')),
verbose=TRUE)## Generating 999 simulations by evaluating expression ...
## 1, 2, 3, ......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## 10 [etd 4:55] .........20 [etd 4:39] .........30 [etd 4:32] ...## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ...## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ...40 [etd 4:36] .....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ....50
## [etd 4:34] ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ...## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..60 [etd 4:38]## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ...70 [etd 4:40] .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .80 [etd 4:41] .........90 [etd 4:34] .........100
## [etd 4:29] ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .....110 [etd 4:25] .........120 [etd 4:20] ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .....130 [etd 4:17] .........140 [etd 4:13] .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ...150
## [etd 4:12] ......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ...## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## 160 [etd 4:09]## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..170 [etd 4:08] ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .....180 [etd 4:04] .....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..190 [etd 4:02] ......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..200
## [etd 4:00] .....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ....210 [etd 3:56] .....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .220 [etd 3:54]## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ........## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .230 [etd 3:52] .........240 [etd 3:48] .........250
## [etd 3:44] .........260 [etd 3:40]## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .........270 [etd 3:37] .........280 [etd 3:33] ......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ...290 [etd 3:30] ........## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .300
## [etd 3:27]## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..310 [etd 3:25]## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .........320 [etd 3:22] ......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ...330 [etd 3:19] .........340 [etd 3:15] ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .350
## [etd 3:12] ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .......360 [etd 3:09] ...## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ......370 [etd 3:06] .......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..380 [etd 3:04] ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .......390 [etd 3:01] ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .....400
## [etd 2:57]## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .........410 [etd 2:54] ........## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .420 [etd 2:51]## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ....430 [etd 2:49] .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ........## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## 440 [etd 2:46]## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .........450
## [etd 2:43] ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ......460 [etd 2:40]## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..470 [etd 2:37] ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .480 [etd 2:35] ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ...## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .490 [etd 2:32] .........500
## [etd 2:29] .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .......510 [etd 2:26]## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .........520 [etd 2:23] .........530 [etd 2:20] .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ....540 [etd 2:17] .........550
## [etd 2:14] ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .....560 [etd 2:11] .........570 [etd 2:08] .........## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## 580 [etd 2:05] ......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ...590 [etd 2:02] .........600
## [etd 1:58] ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ...610 [etd 1:56] .....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ....620 [etd 1:53] .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ........630 [etd 1:50] ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .640 [etd 1:47] ........## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .650
## [etd 1:44] ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .....660 [etd 1:41] ...## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ......670 [etd 1:38] .........680 [etd 1:35] .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..690 [etd 1:32] .........700
## [etd 1:29]## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .......710 [etd 1:26] .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ....720 [etd 1:23] .......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..730 [etd 1:20] ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .....740 [etd 1:17] ...## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ....750
## [etd 1:14] ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .......760 [etd 1:11] .........770 [etd 1:08] .......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..780 [etd 1:05] ...## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ...## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ...790 [etd 1:02] .........800
## [etd 59 sec] ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .......810 [etd 56 sec] ...## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## 820 [etd 53 sec] ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .....830 [etd 50 sec] ........## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .840 [etd 47 sec] .........850
## [etd 44 sec] .....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..860 [etd 41 sec]## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .....870 [etd 38 sec] ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .....880 [etd 36 sec] .....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..890 [etd 33 sec] ......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ...900
## [etd 30 sec] ...## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ..910 [etd 27 sec] ......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ...920 [etd 24 sec] .........930 [etd 21 sec] .........940 [etd 18 sec] .........950
## [etd 15 sec] .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .960 [etd 12 sec] .........970 [etd 9 sec]## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .......## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## .980 [etd 6 sec] .## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ........## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## 990 [etd 3 sec] .....## Warning in ppp(XY$x, XY$y, window = W, marks = marx, check = check &
## needcheck): data contain duplicated points## ... 999.
##
## Done.plot(L_stations_ws_rs99, .-mmean~r, legend=FALSE)
## lty col key
## obs 1 1 obs
## mmean 2 2 mmean
## hi 1 8 hi
## lo 1 8 lo
## label
## obs italic({hat(L)[list(USCR,ws)]^{obs}}(r)-{bar(L)[list(USCR,ws)]}(r))
## mmean italic({bar(L)[list(USCR,ws)]}(r)-{bar(L)[list(USCR,ws)]}(r))
## hi italic({hat(L)[list(USCR,ws)]^{hi}}(r)-{bar(L)[list(USCR,ws)]}(r))
## lo italic({hat(L)[list(USCR,ws)]^{lo}}(r)-{bar(L)[list(USCR,ws)]}(r))
## meaning
## obs observed value of L["USCR","ws"](r) for data pattern
## mmean sample mean of L["USCR","ws"](r) from simulations
## hi upper pointwise envelope of L["USCR","ws"](r) from simulations
## lo lower pointwise envelope of L["USCR","ws"](r) from simulationsOne more thing. Recall that the null hypothesis for the cross-K function is that K_11 = K_22 = K_12. This means that the point pattern of any one set of points alone (or of the points combined) cannot be distinguished. Let’s plot these for USCR and ws.
stations_USCR = ppp(stations_ws$x[stations_ws_shp$type=="USCR"],
stations_ws$y[stations_ws_shp$type=="USCR"],
window=USshape3)
K_stations_USCR = Kest(stations_USCR, correction="translation")
stations_ws = ppp(stations_ws$x[stations_ws_shp$type=="ws"],
stations_ws$y[stations_ws_shp$type=="ws"],
window=USshape3)## Warning in ppp(stations_ws$x[stations_ws_shp$type == "ws"],
## stations_ws$y[stations_ws_shp$type == : data contain duplicated pointsK_stations_ws = Kest(stations_ws, correction="translation")
plot(K_stations_USCR$r, K_stations_USCR$trans - K_stations_ws$trans, type='l')
plot(K_stations_USCR$r, K_stations_USCR$trans - K_stations_ws$trans, type='l')
plot(K_stations_ws$r, K_stations_ws$trans - K_stations_ws$trans, type='l')