polygon - Error for point.in.poly function in the spatialEco package in R stating arguments imply differing number of rows

I am using the point.in.poly function from the spatialEco package. When I ran my original point shapefile and polygon shapefile, I got this error message:

  arguments imply differing number of rows: 36, 0

Here is code to reproduce the error (I cannot share the shapefiles themselves because they are proprietary)

#create point shapefile

points_cor=data.frame(cbind(c(3018, 1069), c(605282.5, 605312.7), c(5022674, 5022730)))

names(points_cor)=c("sta", "long", "lat")
coordinates(points_cor)=~long+lat
proj4string(points_cor)=CRS("+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

#create polygon shapefile

p1=Polygons(list(Polygon(cbind(c(604808.1, 604868.8, 604894.7, 604894.8, 604890.5, 604865.9, 604866.9, 604867.1, 604867.2,604808.1), 
c(5023269, 5023380, 5023366, 5023366, 5023354, 5023239, 5023207, 5023207, 5023206, 5023269)))), '1')
p2=Polygons(list(Polygon(cbind(c(605072.7, 604956.5, 604955.8, 604867.2, 604867.1, 604866.9, 604894.8, 605072.7), 
c(5023274, 5023192, 5023192, 5023206, 5023207, 5023207, 5023366, 5023274)))), '2')

p3=Polygons(list(Polygon(cbind(c(604959.0, 605112.8, 605420.0, 605138.9, 604969.5, 604959.0), 
c(5022650, 5022655, 5022708, 5022550, 5022604, 5022650)))), '3')
p4=Polygons(list(Polygon(cbind(c(605451.4,605515.8, 605452.1, 605451.4), 
c(5023017, 5022954, 5022969, 5023017)))), '4')
p5=Polygons(list(Polygon(cbind(c(605310.3, 605260.8, 605151.7, 605212.8, 605329.4, 
605378.1, 605378.5, 605378.6, 605310.3), 
c( 5022710, 5022691, 5022673, 5022728, 5022785, 5022779, 5022779, 5022779, 5022710)))), '5')
sr=SpatialPolygons(list(p1,p2,p3,p4, p5), proj4string=CRS("+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))
cl_dens=c("C","B", NA, "A", "B")
perturb=rep(NA, 5)

type_cv=rep("F", 5)
vars=data.frame(cbind(cl_dens, perturb, type_cv))
srdf=SpatialPolygonsDataFrame(sr, data.frame(row.names=c('1','2','3','4', '5'), vars))

#run point.in.poly with all variables of the polygon shapefile
pts.poly <- point.in.poly(points_cor, srdf)

This gives the following error message:

Error in data.frame(z@data, stats::na.omit(sp::over(pts, polys))) : 
  arguments imply differing number of rows: 2, 0

If I omit the column/variable with only NA values, it does work with following output:

#run point.in poly with only the variables containing non-empty values
pts.poly <- point.in.poly(points_cor, srdf[,c(1,3)])
> pts.poly
class       : SpatialPointsDataFrame 
features    : 2 
extent      : 605282.5, 605312.7, 5022674, 5022730  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 3

names       :  sta, cl_dens, type_cv 
min values  : 1069,       B,       F 
max values  : 3018,       B,       F 

Answer

I changed the spatialEco::point.in.poly function so that it now retains columns containing all NA values. The behavior of the function was as intended, I just did not anticipate NA columns. I adhere to the philosophy that a function should not unnecessarily modify data. The root of the problem was in how stats::na.omit was being applied. Note that I also added an argument "poly.id" that allows specification of a column in the SpatialPolygonsDataFrame containing unique ID's. If the argument is left undefined (NULL) then unique ID's will be based on the rownames of the SpatialPolygonsDataFrame object.

This change will be in the next release of spatialEco (in a week or so). In the meantime here is the modified function.

#' @title Point and Polygon Intersect
#' @description Intersects point and polygon feature classes and adds polygon attributes to points
#'
#' @param pts      sp SpatialPointsDataFrame or SpatialPoints object

#' @param polys    sp SpatialPolygonsDataFrame object
#' @param poly.id  Name of column contaning unique polygon ID's  
#'
#' @return A SpatialPointsDataFrame with intersected polygon attributes
#'
#' @note Depends: sp
#' @note If poly.id is NULL a column "pids" is created based on the rownames of the polygons
#'
#' @author Jeffrey S. Evans  
#'

#' @examples 
#' require(sp)
#' data(meuse)
#' coordinates(meuse) = ~x+y
#' meuse@data$test.na <- NA
#' 
#' sr1=Polygons(list(Polygon(cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409, 
#'   180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676, 
#'   332618, 332413, 332349)))),'1')
#' sr2=Polygons(list(Polygon(cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437, 

#'   179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683, 
#'   331133, 331623, 332152, 332357, 332373)))),'2')
#' sr3=Polygons(list(Polygon(cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875, 
#'   179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
#'   c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004, 
#'   329783, 329665, 329720, 329933, 330478, 331062, 331086)))),'3')
#' sr4=Polygons(list(Polygon(cbind(c(180304, 180403,179632,179420,180304),
#'   c(332791, 333204, 333635, 333058, 332791)))),'4')
#' sr=SpatialPolygons(list(sr1,sr2,sr3,sr4))
#' srdf=SpatialPolygonsDataFrame(sr, data.frame(row.names=c('1','2','3','4'), PIDS=1:4, y=runif(4)))

#'    srdf@data$pid <- srdf@data$PIDS + 100 
#' 
#' # Intersect points with polygons and add polygon IDS to pts@@data. 
#'   pts.poly <- point.in.poly(meuse, srdf, poly.id = "pid") 
#'     head(pts.poly)
#' 
#' # Point counts for each polygon
#' tapply(pts.poly$lead, pts.poly$PIDS, FUN=length)  
#'
#' @export

point.in.poly <- function(pts, polys, poly.id = NULL) {
    if (!inherits(polys, "SpatialPolygonsDataFrame")) 
        stop("polys must be a SpatialPolygonsDataFrame object")
    if ((inherits(pts, "SpatialPointsDataFrame") | inherits(pts, "SpatialPoints")) == FALSE) 
        stop("pts must be a SpatialPointsDataFrame object")
    if(!is.null(poly.id)) {
      if( length(unique(polys@data[,poly.id])) != nrow(polys) ) stop("poly.id not unique")
    }   
    if(is.null(poly.id)) { 
      polys@data$pids <- rownames(polys@data)

      poly.id = "pids"
    }     
    z <- pts[!is.na(sp::over(pts, sp::geometry(polys))),]
      z@data <- data.frame(sp::over(pts, polys))
        row.ids <- rownames(z@data)
        z@data <- data.frame(z@data[,which(names(polys) %in% poly.id)])
          names(z) <- poly.id
            rownames(z@data) <- row.ids       
          z@data <- data.frame(z@data, merge(z, pts@data, by="row.names", all.x = FALSE))
        names.rm <- c("x", "y", "Row.names", "optional", paste0(poly.id, ".1"))

      z@data <- z@data[,-which(names(z) %in% names.rm)]
    z@proj4string <- pts@proj4string
  return( z )
}