[R-sig-Geo] Cropping a set of points with a polygon

[Tom Gottfried]

You might be interested in over() from sp and the rgeos package.

Tom

Am 01.05.2012 05:24, schrieb steven mosher:
> 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]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

-- 
Technische Universität München
Department für Pflanzenwissenschaften
Lehrstuhl für Grünlandlehre
Alte Akademie 12
85350 Freising / Germany
Phone: ++49 (0)8161 715324
Fax:   ++49 (0)8161 713243
email: tom.gottfried at wzw.tum.de
http://www.wzw.tum.de/gruenland