2

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 
  • There is no limitation to the number of columns and i imagine that there is something in the data that is causing a data mismatch in the NA handling. Can you make you data avalible so I can troubleshoot the issue? – Jeffrey Evans Sep 29 '17 at 22:56
  • Do NA values in a given column cause errors? There are some columns/variables with no values for some features. – user3386170 Sep 30 '17 at 13:54
  • I have updated the post having tested out different columns following on your lead that NA may be problematic. My dataset has columns without any values at all (the polygon is a portion of a larger dataset), i.e. all NAs. When I omit those columns specifically, the function works (and I've keep some columns where there are some NA values). Was that the problem? – user3386170 Oct 1 '17 at 1:41
  • @JeffreyEvans : given what I have added into my answer and my comment, could you confirm that I found the root of the problem I was encountering? – user3386170 Oct 18 '17 at 19:23
3

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  <jeffrey_evans@@tnc.org>
#'
#' @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 )
}
  • Please update to the latest version of the package! As of 3/12/2018, it is at version 1.1-0. I completely rewrote the function. – Jeffrey Evans Mar 19 '18 at 16:21
  • Do you automatically update the citation that goes with the package? I have updated the package today and am putting together a bibliography but noticed that citation(spatialEco) refers to v. 0.0.1-7 and year 2017. Just letting you know in case it matters to you. – user3386170 Mar 20 '18 at 15:35
  • @user3386170 thanks for pointing this out. I was in a flurry to get this version out the door and forgot to update the citation. – Jeffrey Evans Mar 20 '18 at 18:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.