Hello,
I am using the following code to create ppp files from csv data and map shape
files, but I am getting some errors which I have been unable to fix by
searching them online:
library(spatstat)
library(maps)
library(maptools)
NYC2<-readShapePoly("nybb.shp") # this is a map of the NYC boroughs without
waterways and no census tract divisions (but it does include lines separating
the 5 boroughs)
plot(NYC2)
NYBorough<-readShapePoly("NYBoroughsShapeMerge.shp") # this is a map of the NYC
boroughs with census tract divisions
plot(NYBorough)
myfile<-read.csv("myfile.csv")
# my csv file contains 4 columns (ID, variable_type, X, Y)
# my goal is to determine whether individuals of variable_type=1 are spatially
distributed more than what would be expected due to chance, if individuals of
variable_type=0 are spatially distributed more than what would be expected due
to chance, and finally if the spatial distribution of variable_type=1 differs
from that of variable_type=0.
x<-NYBorough@polygons[[1]]@Polygons
# ppp(x.coordinates, y.coordinates, x.range, y.range)
coords<-slot(x[[1]],"coords")
coords
coords<-coords[19:1,]
coords<-coords[-1,]
coords
u<-myfile$type==0
Type_0<-ppp(myfile$X[u],myfile$Y[u],
window=owin(poly=list(x=coords[,1],y=coords[,2])))
Type_1<-ppp(myfile$X[!u],myfile$Y[!u],
window=owin(poly=list(x=coords[,1],y=coords[,2])))
ERROR Message:
Warning message:
In ppp(myfile$X[u], myfile$Y[u], window = owin(poly = list(x = coords[, :
217 points were rejected as lying outside the specified window
Warning message:
In ppp(myfile$X[!u], myfile$Y[!u], window = owin(poly = list(x = coords[, :
435 points were rejected as lying outside the specified window
# I only have 652 points, so all of my points are rejected as lying outside the
specified window
Can someone please tell me what I am doing wrong. The code I was using as a
template was only looking at one county so it had a defined border which was
somewhat rectangular. Could it be because the shapefile I am using for NYC has
too many borders? Is it the coords code that is incorrect?
Any advice or suggestions would be greatly appreciated.
In case it is helpful, the code that follows is:
qtype_1<-quadratcount(type_1,10,10)
plot(type_1,pch="+",main=Title")
plot(type_1,add=TRUE,col="red",ltw=2,cex=1.2)
polygon(coords,lwd=2)
points(type_1,pch=".")
Thanks,
Abby
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
