L11 Spatial Methods: Point Patterns

Plan: 1) Show some basic plots of restaurant data 2) Convert restaurants to PPP 3) Count how many Italian restaurants within given radii 4) Plot count against radii 5) Compare to randomly drawing restaurants Note: there are much better ways to do this using specific functions from the SpatStat package (and others). However, I think it’s helpful to start by doing everything “by hand” before relying on packages.

Load necessary libraries:

library(maptools)
library(rgdal)
library(sp)
library(spatstat)

Loading required package: nlme
Loading required package: rpart

spatstat 1.51-0       (nickname: <U+393C><U+3E31>Poetic Licence<U+393C><U+3E32>) 
For an introduction to spatstat, type <U+393C><U+3E31>beginner<U+393C><U+3E32> 

Bring in point data and plot

rest_data<-read.table("data/menupages_Wb126.csv", header=T, sep=",",stringsAsFactors =FALSE)
rest_spoints<-SpatialPointsDataFrame(cbind(rest_data$x,rest_data$y), rest_data, coords.nrs=c(4,5),match.ID=TRUE)
#Read shapefiles with readOGR  
nyc_tracts<-readOGR("data","manhattan_tracts")

OGR data source with driver: ESRI Shapefile 
Source: "data", layer: "manhattan_tracts"
with 221 features
It has 46 fields

plot(nyc_tracts,main="Manhattan Restaurants 2008")
points(rest_spoints)

Subset Plots

plot(nyc_tracts,main="Italian Restaurants")
points(rest_spoints[rest_spoints$cuisName1=="Italian",],col="red")

Now look at statistics on point patterns. In order to take advantage of spatstat distance methods we use a projection (NAD1983 Long Island, feet). This allows simple Euclidean distances.

min_x<-min(rest_spoints$x_feet)
max_x<-max(rest_spoints$x_feet)
min_y<-min(rest_spoints$y_feet)
max_y<-max(rest_spoints$y_feet)
W <- owin(xrange=c(min_x,max_x),yrange=c(min_y,max_y)) #a window defined by limits of data
allrests.ppp<-ppp(rest_spoints$x_feet,rest_spoints$y_feet,window=W)

data contain duplicated points

plot(allrests.ppp)

italian.ppp<-allrests.ppp[rest_spoints$cuisName1=="Italian",]
chinese.ppp<-allrests.ppp[rest_spoints$cuisName1=="Chinese",]
pizza.ppp<-allrests.ppp[rest_spoints$cuisName1=="Pizza",]
#plot(allrests.ppp,main="Chinese Restaurants")
print(paste("There are ",length(chinese.ppp$x)," Chinese restaurants and ",length(pizza.ppp$x)," pizza restaurants."))

[1] "There are  338  Chinese restaurants and  316  pizza restaurants."

plot(chinese.ppp)

plot(italian.ppp)

plot(pizza.ppp)

Look at count of Chinese restaurants within given radii

neighDist<-pairdist(chinese.ppp)
radii_set<-seq(200,3600,200)
chinese_mat<-array(,c(length(radii_set),2))
for (i in 1:length(radii_set)) {
  chinese_mat[i,1]<-radii_set[i]
  chinese_mat[i,2]<-mean(rowSums((neighDist<radii_set[i]), na.rm=TRUE))-1 #remove own distance
}
plot(chinese_mat,xlab="radius",ylab="Count",main="Count Chinese Restaurants within Radius")

neighDist<-pairdist(pizza.ppp)
pizza_mat<-array(,c(length(radii_set),1))
pizza_mat<-cbind(chinese_mat[,1],pizza_mat)
for (i in 1:length(radii_set)) {
  pizza_mat[i,2]<-mean(rowSums((neighDist<radii_set[i]), na.rm=TRUE))-1 #remove own distance
}
plot(chinese_mat,xlab="radius",ylab="Count",main="Chinese (red) and Pizza (green)",
     ylim=range(c(chinese_mat[,2],pizza_mat[,2])),col="red")
par(new=T)

plot(pizza_mat,xlab="",ylab="",axes=FALSE,col="green")

Now compare Chinese to simulated measure

numSims<-1000
sim_mat<-array(,c(length(chinese_mat[,1]),numSims))
dim(sim_mat)

[1]   18 1000

for (i in 1:numSims) {
  sim.ppp<-allrests.ppp[sample(1:length(allrests.ppp$x),length(chinese.ppp$x),replace=F)]
  neighDist_sim<-pairdist(sim.ppp)
  for (j in 1:length(radii_set)) {
    sim_mat[j,i]<-mean(rowSums(neighDist_sim<radii_set[j], na.rm=TRUE))-1
  }
}
conf_mat<-array(,c(length(chinese_mat[,1]),2))
for (j in 1:length(radii_set)) {
  conf_mat[j,1]<-quantile(sim_mat[j,],probs=.025,na.rm=TRUE)
  conf_mat[j,2]<-quantile(sim_mat[j,],probs=.975,na.rm=TRUE)
}
chinese_mat<-cbind(chinese_mat,conf_mat)
plot(chinese_mat[,1],chinese_mat[,2],ylim=range(c(chinese_mat[,2],chinese_mat[,3],chinese_mat[,4])),
     xlab="Radius",ylab="Count",pch=20,col="red")
par(new=T)

plot(chinese_mat[,1],chinese_mat[,3],ylim=range(c(chinese_mat[,2],chinese_mat[,3],chinese_mat[,4])),
     axes=FALSE,xlab="",ylab="",main="Count Chinese Restaurants within Radius",pch=25,cex=0.5)
par(new=T)

plot(chinese_mat[,1],chinese_mat[,4],ylim=range(c(chinese_mat[,2],chinese_mat[,3],chinese_mat[,4])),
     axes=FALSE,xlab="",ylab="",pch=24,cex=0.5)

Do same for pizza restaurants

sim_mat<-array(,c(length(pizza_mat[,1]),numSims))
dim(sim_mat)

[1]   18 1000

for (i in 1:numSims) {
  sim.ppp<-allrests.ppp[sample(1:length(allrests.ppp$x),length(pizza.ppp$x),replace=F)]
  neighDist_sim<-pairdist(sim.ppp)
  for (j in 1:length(radii_set)) {
    sim_mat[j,i]<-mean(rowSums(neighDist_sim<radii_set[j], na.rm=TRUE))-1
  }
}
conf_mat<-array(,c(length(pizza_mat[,1]),2))
for (j in 1:length(radii_set)) {
  conf_mat[j,1]<-quantile(sim_mat[j,],probs=.025,na.rm=TRUE)
  conf_mat[j,2]<-quantile(sim_mat[j,],probs=.975,na.rm=TRUE)
}
pizza_mat<-cbind(pizza_mat,conf_mat)
plot(pizza_mat[,1],pizza_mat[,2],ylim=range(c(pizza_mat[,2],pizza_mat[,3],pizza_mat[,4])),
     xlab="Radius",ylab="Count",col="red",cex=1,pch=20)
par(new=T)

plot(pizza_mat[,1],pizza_mat[,3],ylim=range(c(pizza_mat[,2],pizza_mat[,3],pizza_mat[,4])),
     axes=FALSE,xlab="",ylab="",main="Count Pizza Restaurants within Radius",pch=25,cex=0.5)
par(new=T)

plot(pizza_mat[,1],pizza_mat[,4],ylim=range(c(pizza_mat[,2],pizza_mat[,3],pizza_mat[,4])),
     axes=FALSE,xlab="",ylab="",pch=24,cex=0.5)