Function Window of spatstat lose spatialpoints

Question

I have a code for my test, with a short data set (20 points) and a small spatial area. It's working fine.

After the test, I wanted to increase the spatial area, I defined a bigger polygon, plotted the 20 points and they are within the area.

However the function Window of spatstat didn't find the points, there appears 0 points.

What's happening here?

day<-c(22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22)
longitude<- c(617237.1,495004.2,618067.6,617088.9,616980.4,358600.5,617759.7,617188.4,617029.6,649704.0,284809.5,358567.5,711395.7,617637.2,617069.2,358272.2,575173.3,726735.2,684353.3,636383.3)
latitude<-c(4177364,4284973,4177484,4177504,4177411,4266954,4177658,4177359,4177455,4130436,4400643,4266136,4185628,4177310,4177367,4266181,4126490,4109880,3989603,4074784)
data<-cbind(day,longitude,latitude)
dataDF<-as.data.frame(data)
coordinates(dataDF)<-  ~longitude+latitude
sp::proj4string(dataDF) <-CRS("+proj=utm +zone=26 +ellps=intl +towgs84=-104,167,-38,0,0,0,0 +units=m +no_defs")
# small polygon 
y_coord <- c(3900000,3900000,4500000 ,4500000   )
x_coord <- c(2e+05,  8e+05, 8e+05,2e+05)
xym1 <- cbind(x_coord, y_coord)
xym1
library(sp)
p = Polygon(xym1)
ps = Polygons(list(p),1)
quadrat1 = SpatialPolygons(list(ps))
sp::proj4string(quadrat1) <-CRS("+proj=utm +zone=26 +ellps=intl +towgs84=-104,167,-38,0,0,0,0 +units=m +no_defs")
plot(quadrat1,axes=TRUE,col="green")
plot(dataDF,add=TRUE,col="red")# points within extent
p.ppp <- as(dataDF, "ppp")      
class(p.ppp)
ma1    <- as(quadrat1, "owin")
#Window works fine
marks(p.ppp)  <- NULL
Window(p.ppp) <- ma1#
plot(p.ppp, main=NULL, cols=rgb(0,0,0,.2), pch=20)
#big polygon 
y_coord <- c(2900000,2900000,5900000 ,5900000   )
x_coord <- c(-15e+05,  2e+06, 2e+06,-15e+05)
xym2 <- cbind(x_coord, y_coord)
xym2
library(sp)
p2 = Polygon(xym2)
ps2 = Polygons(list(p2),1)
quadrat2 = SpatialPolygons(list(ps2))
sp::proj4string(quadrat2) <-CRS("+proj=utm +zone=26 +ellps=intl +towgs84=-104,167,-38,0,0,0,0 +units=m +no_defs")

plot(quadrat2,axes=TRUE,col="green")
plot(dataDF,add=TRUE,col="red")# 20 points within
#Window didn-t work
p.ppp <- as(dataDF, "ppp")      
ma2    <- as(quadrat2, "owin")#big polygon
marks(p.ppp)  <- NULL
Window(p.ppp) <- ma2# 0 points?
plot(p.ppp, main=NULL, cols=rgb(0,0,0,.2), pch=20)# 0 points

Ege Rubak Ege Rubak · Accepted Answer · 2019-02-08T21:01:08

Edit: This was a bug in the handling of a special case in spatstat (rectangle encoded as a general polygon) which has now been fixed. So from version 1.58-2.026 onward this should no longer occur. It was not related to computer precision at all. The answer below still contains the valid point of encoding a observation window (owin) explicitly as a rectangle rather than using a general polygon which introduces unnecessary complexity.

This error appears to be due to a problem with computer precision when using numerically large coordinates in the polygon. Polygon computations and point-in-polygon tests are difficult to implement with high precision on a computer, so everything polygon related in spatstat is scaled to integers and handled by some C++ routines. Apparently something is problematic here. We will investigate.

Since your polygon is actually a simple rectangle you can avoid the problem by specifying it directly in spatstat rather than converting a spatial polygon to a polygonal window in spatstat.

library(spatstat)
x <- c(617237.1,495004.2,618067.6,617088.9,616980.4,358600.5,617759.7,617188.4,617029.6,649704.0,284809.5,358567.5,711395.7,617637.2,617069.2,358272.2,575173.3,726735.2,684353.3,636383.3)
y <- c(4177364,4284973,4177484,4177504,4177411,4266954,4177658,4177359,4177455,4130436,4400643,4266136,4185628,4177310,4177367,4266181,4126490,4109880,3989603,4074784)
y_coord <- c(2900000,2900000,5900000,5900000)
x_coord <- c(-15e+05,  2e+06, 2e+06,-15e+05)

Problematic example:

W1 <- owin(poly=cbind(x_coord, y_coord))
W1
#> window: polygonal boundary
#> enclosing rectangle: [-1500000, 2e+06] x [2900000, 5900000] units
X1 <- ppp(x, y, window = W1)
#> Warning: 20 points were rejected as lying outside the specified window
X1
#> Planar point pattern: 0 points
#> window: polygonal boundary
#> enclosing rectangle: [-1500000, 2e+06] x [2900000, 5900000] units
#> *** 20 illegal points stored in attr(,"rejects") ***

Working example:

W2 <- owin(x_coord[1:2], y_coord[2:3])
X2 <- ppp(x, y, window = W2)
X2
#> Planar point pattern: 20 points
#> window: rectangle = [-1500000, 2e+06] x [2900000, 5900000] units

^{Created on 2019-02-08 by the reprex package (v0.2.1)}

Function Window of spatstat lose spatialpoints

1 Answers