I have a  dataframe of points that I would like to "crop"  or classify
according to
a test of whether then are inside a polygon or outside a polygon.

I use "sp"  and  "rgdal"

The polygon is supplied as a shapefile and comes with its own  projection.
my datapoints are all in the projection listed below. ("+proj=longlat
+datum=WGS84")

what I'd like to do is create a function  "crop" where I feed the function,
the datapoints and the shapefile
and I get back a dataframe of only those points that lie inside the
polygon. polygons will never have holes.
I'm thinking that points.in.polygon() is going to useful, but wondering if
there is a better way or if I am missing something


 library(sp)
 library(rgdal)
#   Generate some random data that is nearby the shapefile I will use to
crop
dataPoints <- matrix(NA,nrow=100, ncol = 2)

dataPoints[ ,1] <- sample( x = seq(-80, -45, .1), size = 100, replace =
FALSE)
dataPoints[ ,2] <- sample( x = seq(40, 70, .1),   size = 100, replace =
FALSE)
dataPoints      <- SpatialPoints(dataPoints)

df              <- data.frame(Id = as.character(1:100))

pointsDF        <- SpatialPointsDataFrame(dataPoints,df,proj4string =
CRS("+proj=longlat +datum=WGS84" ))

#  read in the shapefile.. data listed for information pruposes
 X              <- readOGR(dsn = getwd(), layer = "New_Shapefile")

ogrInfo(dsn=getwd(),layer = "New_Shapefile")
Source: "G:/ShapTest", layer: "New_Shapefile"


  Driver: ESRI Shapefile number of rows 1
  Feature type: wkbPolygon with 2 dimensions
  +proj=utm +zone=20 +datum=NAD83 +units=m +no_defs
  Number of fields: 1
  name type length typeName
  1   Id    0      6  Integer


  dput(X)


new("SpatialPolygonsDataFrame"
    , data = structure(list(Id = 1L), .Names = "Id", row.names = 0L,
class = "data.frame")
    , polygons = list(<S4 object of class structure("Polygons",
package = "sp")>)
    , plotOrder = 1L
    , bbox = structure(c(227198.913276668, 5878401.92368424, 968033.728288047,
6714486.92918889), .Dim = c(2L, 2L), .Dimnames = list(c("x",
"y"), c("min", "max")))
    , proj4string = new("CRS"
    , projargs = " +proj=utm +zone=20 +datum=NAD83 +units=m +no_defs
+ellps=GRS80 +towgs84=0,0,0"
)
)

	[[alternative HTML version deleted]]

