And that is how it is documented: the help obtained by
methods?overlay
tells you about the method used:
x = "SpatialPoints", y = "SpatialPolygons" returns a numeric
vector of length equal to the number of points; the number is
the index (number) of the polygon of ‘y’ in which a point
falls; NA denotes the point does not fall in a polygon; if a
point falls in multiple polygons, the last polygon is
recorded.
so
plot(grd[!is.na(clip)])
should give you the pixels inside the polygon(s).
On 10/28/2010 10:23 PM, Peter Larson wrote:
> Hello all,
>
> I have a grid of points and a spatial polygon shape file. I want to
> clip the grid so that I only have the point which lie within the
> shapefile.
>
> I tried the following code:
>
> ## create a grid onto which we will interpolate:
> ## first get the range in data
> x.range <- c(37,50)
> y.range <- c(29,38)
>
> ## now expand to a grid with 500 meter spacing:
> grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=.01),
> y=seq(from=y.range[1], to=y.range[2], by=.01) )
>
> ## convert to SpatialPixel class
> coordinates(grd) <- ~ x+y
> gridded(grd) <- TRUE
>
> ## Read in Shapefile
> ccShapet=readShapePoly('IRQ_adm0') # here cbg00barncnty is an ArcGIS shapefile
>
> ## Clip
> clip <- overlay(grd, ccShapet)
>
> #Plot
> plot(clip)
>
> BUT it does not work. I just get a vector of number and NA's, but not
> the clipped grid.
>
> Any ideas?
>
> Thanks
>
> Pete
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics e.pebe...@wwu.de
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
