R help - Jan 2013 - points rejected as lying outside the specified window

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 at 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

This is a question about the "spatstat" package and as such would have
been more appropriately directed to the R-sig-Geo list.

(1) Your code seems unnecessarily complicated; I haven't followed it
through in detail however.

(2) You did not get an "error", you got a "warning".  This
warning means
what it says:   A number of the points that you are trying to make into
a planar point pattern lie out side of the window that you specified.

Since you say that *all* of your points lie outside the specified window
one might guess that the coordinate system of the points in
"myfile.csv"
(God how I hate that "mythis" and "mythat" terminology that
Micro$soft
has introduced!) is different from the coordinate system that is used in the
"NYBorough" data set.

Look at

     range(coords[,1])
     range(coords[,2])
     range(myfile$X)
     range(myfile$Y)

to get some insight into the nature of the coordinates.

If the coordinates seem to be at least roughly commensurate
your could try:

WR <- ripras(x=myfile$X[u],y=myfile$Y[u],shape="rectangle")
Type_0 <- ppp(myfile$X[u],myfile$Y[u],window=WR)
WP <- owin(poly=list(x=coords[,1],y=coords[,2]))
plot(Type_0)
plot(WP,add=TRUE)

which might be revealing.

Other than that it's difficult to say what your problem is without having
access to your actual data.

     cheers,

         Rolf Turner

On 01/30/2013 09:29 AM, Abby Rudolph wrote:
> 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 at 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=".")