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

[Andy Bunn]

> -----Original Message-----
> From: Edzer J. Pebesma [mailto:e.pebesma at geo.uu.nl]
> Sent: Friday, April 07, 2006 5:30 AM
> To: Tim Keitt
> Cc: Andy Bunn; r-sig-geo at stat.math.ethz.ch; uleopold at science.uva.nl
> Subject: Re: [R-sig-Geo] (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)


Edzer, Tim, and Roger: Thanks so much for the magic. That solution does
work. The only sticky bit is setting fullgrid(x) to FALSE. In my application
this takes about 25 seconds as the SpatialGridDataFrame is 1405x621. Keep up
the good work fellers.

Thanks again, Andy




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