Find points that are %in% a SpatialPolygonsDataFrame

spatial_in.r

library(sp)
library(raster)

setClass("coords", contains="matrix")

# x should be a matrix. The first column holds longitudes, the
# second holds latitudes.
coords <- function(x) {
  if (!inherits(x, "matrix") && !inherits(x, "data.frame")) {
    stop(paste0("Don't know how to make coordinates from an object ",
                "of class ", class(x), "."))
  }
  if (inherits(x, "data.frame")) {
    x <- as.matrix(x)
  }
  if (ncol(x) != 2) {
    stop("x should have 2 columns - 1st is lon, 2nd is lat.")
  }
  colnames(x) <- c("lon", "lat")
  rownames(x) <- NULL
  new("coords", x)
}

setMethod("%in%", signature(x="coords", table="SpatialPolygonsDataFrame"),
  function(x, table) {
    ret <- rep(FALSE, nrow(x))
    if (length(ret) > 0) {
      # Get the points in the bounding box.
      in_bb <- as.vector(which( (x[, "lon", drop=FALSE] >= xmin(table)) &
                           (x[, "lon", drop=FALSE] <= xmax(table)) &
                           (x[, "lat", drop=FALSE] <= ymax(table)) &
                           (x[, "lat", drop=FALSE] >= ymin(table)) ))

      # Assume the projection. It would be nice to detect this and project or error.
      sp_coords <- SpatialPointsDataFrame(
        coords=x[in_bb,,drop=FALSE],
        proj4string=CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"),
        data=as.data.frame(in_bb))
      ret[raster::intersect(sp_coords, table)$in_bb] <- TRUE
    }
    ret
  })