point.in.poly (spatialEco) returns non-intersecting points

Question

I just updated the spatialEco package for R, and now the point.in.poly function behaves differently from previous behavior. Previously, it returned only points that intersected at least one of the polygons in the polygon layer. Now, it returns all points, including those that do not intersect with any polygons (the polygon fields are NA for the non-intesecting points). It seems that an extra step is now necessary to retain only the intersecting points, which doesn't seem desirable to me. It is not too burdensome to add the extra step, but it does break old code that uses point.in.poly. Is there a way to retain only the intersecting points in the output from the point.in.poly function?

Below is some example code that illustrates what I mean.

library(rgeos)
library(sp)
library(spatialEco)
library(raster)

coords <- matrix(c(0,0,0,10,10,10,10,0,0,0),5,2,byrow=TRUE)
SA <- SpatialPolygonsDataFrame(SpatialPolygons(list(Polygons(list(Polygon(coords)),1))),data.frame(pol=1))

pts <- spsample(SA,50,"regular")
pts <- SpatialPointsDataFrame(pts,data.frame(ptid=1:length(pts)))
smallSA <- gBuffer(SA,width= - 1,byid=TRUE)
grd <- raster(extent(smallSA),resolution=2)
grd <- as(grd,"SpatialPolygons")
grd <- SpatialPolygonsDataFrame(grd,data.frame(pol=1:length(grd)))
plot(SA)
plot(grd,add=TRUE)
plot(pts,add=TRUE,col="blue",pch=20)
pts2 <- point.in.poly(pts,smallSA) # includes non-intersecting points!
plot(pts2,add=TRUE,col="red",pch=20)
length(pts)
length(pts2)

# remove non-intersecting points
pts3 <- pts2[-which(is.na(pts2$pol)==TRUE),] 
length(pts3)
plot(pts3,add=TRUE,col="blue",pch=20)

Answer 1

You can use st_intersection() from sf, which retains the data from both data sets and keeps only those points actually intersecting your polygon geometry. The sole drawback here is that your Spatial* objects must be made available as sf objects, which can be handled through st_as_sf(). The resulting object can be transformed back to Spatial* using as() as shown below. See also this related answer.

library(sf)

tmp = st_intersection(st_as_sf(pts), st_as_sf(smallSA))
pts4 = as(tmp, "Spatial")

plot(SA)
plot(grd, add = TRUE)
plot(pts4, col = "blue", pch = 20, add = TRUE)

identical(pts3, pts4)
# TRUE