[R-sig-Geo] (no subject); Clipping/Subsetting Sp objects (was no subject)

[Roger Bivand]

On Fri, 7 Apr 2006, Edzer J. Pebesma wrote:

> 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.

Right, and actually also returns which Polygons object the grid centre
belongs to, so that making summary statistics on the pixels by polygons is
a matter of using one of the *apply() functions. The overlay() methods are 
pretty powerful, of course not like GOES for very large data sets, but 
surprisingly fast all the same - certainly faster than implementing robust 
bindings to external software. Another route to this is:

http://starspan.casil.ucdavis.edu/

Roger

> 
> 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
> >>
> >>_______________________________________________
> >>R-sig-Geo mailing list
> >>R-sig-Geo at stat.math.ethz.ch
> >>https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>    
> >>
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no