[R-sig-Geo] Warning message: points were rejected as lying outside the specified window

Rolf Turner r.turner at auckland.ac.nz
Sat May 24 02:05:56 CEST 2014


On 24/05/14 05:55, Hossain, Md wrote:

> With sincere appreciation for all the suggestions so far I got, and thought it's better to keep you all updated.
> To my datasets (coded as "points"), I did the point-in-polygon test with the R code:
>
> library("maps")
> library ("sp")
> library("spatstat")
> library(SDMTools)
> usmap <- map("usa", fill=TRUE)
> polypnts = cbind(x=usmap$x, y=usmap$y)
> points = cbind(x=my_usdat$x, y=my_usdat$y)
> plot(rbind(polypnts, points))
> out = pnt.in.poly(points,polypnts)
> sum(out==0)
>> [1] 1746
> sum(out==1)
>> [1] 0
>
> It shows that in my dataset there is no points outside of USA map. So the problem I was getting, mainly because of not using the correct coordinate reference system (CRS). It prompted me to apply another strategy that does not require to specify any CRS, posted by Rolf Turner as:
>
> us.df <- with(usmap,data.frame(x=rev(x),y=rev(y)))
> us.df <- unique(us.df)
> B <- list(list(x=us.df$x, y=us.df$y))
> owin(poly=B)
>
> But, got the error message as:
> Error in owin(poly = B) :
>    poly must be either a list(x,y) or a list of list(x,y)
>
> After removing all the "NA" cases by the code:
>
> us.df <- us.df[is.na(us.df$x)==FALSE, ]
> B <- list(list(x=us.df$x, y=us.df$y))
> owin(poly=B)
>
> Got another error message as:
> Error in owin(poly = B) : Area of window is negative;
>   check that all polygons were traversed in the right direction
>
> It seems I am very close, any idea please!

Not really sure what you are up to, but I think a major mistake is 
chucking out the NAs from us.df.

These NAs are there to separate different, distinct, polygons (the 
overall outline of the USA plus various lakes and islands).  This is a 
commonly used convention in GIS data it would seem.  Hellishly kludgy 
and inconvenient, but there's nothing we can do about that.

To make these data into an owin object you need to create a *list*, (say 
"L") each entry of which is a data frame consisting of one of the 
sub-polygons (the parts of us.df "sandwiched" between NA rows).

Then use owin(poly=L).

cheers,

Rolf Turner



More information about the R-sig-Geo mailing list