(no subject); Clipping/Subsetting Sp objects (was no subject)

Andy, Tim, the magic phrases (following your example) are:

clip.sp = SpatialPolygons(list(Polygons(list(Polygon(clip)), ID="clip")))
fullgrid(x) = FALSE
x.clip = x[!is.na(overlay(x, clip.sp)),]
image(x.clip,col="blue",add=T)

the first command creates a SpatialPolygons object from
the coordinates. Yes, this seems complex, but the class
does allow for real-world polygons, having multiple parts,
holes, islands, etc.

The second converts x from SpatialGridDataFrame into SpatialPixelsDataFrame,
which is basically an object with points and their coordinates that 
happen to lie
on a grid layout (instead of a full 2D matrix with NA's on the cells 
outside the
study region).

The third uses the overlay() method, which basically does a point-in-polygon
and returns NA for points outside any of the polygons.

I'm thinking, right now, how to do this without the fullgrid(x)=F statement,
but selection for a SpatialGridDataFrame needs rows/cols, as
in x[row_sel, col_sel, attribute_sel], so can't be done with a list of
pixel indexes. Maybe it should, but I need to be pushed for this. And
that has it's own use. Perhaps here overlay should return a matrix in this
case, and [ should accept (and understand) a matrix as first argument.

Tim Keitt wrote:

We need R wrappers for:

http://www.vividsolutions.com/jts/jtshome.htm

or

http://geos.refractions.net/

My preference would be the latter.

THK
 

"We need", as in "I would like to see", yes, I agree, I would also like to
see this. In the sense of "I would put my own resources on this", I doubt
whether I would do this, or work on a (better?) interface to PostGIS.
After all, GEOS was primarily written to provide the spatial 
functionality there.

Currently we have in sp overlay methods (which mainly do point-in-polygons,
and are in development, comments more than welcome!) to combine Spatial
objects of different classes (I'm pretty sure that overlay can, or 
should be able to,
deal with Ulrich's problem posted last week to statsgrass; haven't had the
opportunity to look into it because of computer crashes) -- we also have
polygon/polygon overlay (polygon clipping) using gpclib, interfaced to sp
classes. (Is it still in spgcplib, Roger?)

What exacty do you believe, Tim, would a direct interface to GEOS add to
what we now have?

I believe that a next challenge is to make R work with massive spatial 
data sets;
rgdal does a lot for this, but for polygon and point data (think of laser
altimetry data, yielding e.g. 1e7-1e8 point observations) we may need
an out-of-memory processor (PostGIS?) that can fast retrieve selections
in local search neigbhourhoods fast (GiST?).
--
Edzer

On Thu, 2006-04-06 at 16:38 -0400, Andy Bunn wrote:
 

List,

I want to clip (subset) an sp object SpatialGridDataFrame using a polygon.
For instance in this example I want to create a new object of class
SpatialGridDataFrame that is clipped to the area inside the polygon on the
map.

 data(meuse.grid)
 coordinates(meuse.grid) <- ~x+y
 gridded(meuse.grid) <- T
 x = as(meuse.grid, "SpatialGridDataFrame")
 clip <- cbind(x=c(180000,180000,180500,180500,180000),
               y=c(330500,331000,331000,330500,330500))
 image(x, "dist")
 polygon(clip)
 # y = x[...clip?]


How can I go about it?

Thanks in advance, Andy