[R-sig-Geo] SpatialPolygons to SpatialGrid

Roger Bivand Roger.Bivand at nhh.no
Mon Aug 18 18:27:14 CEST 2008


On Sun, 17 Aug 2008, Murray Richardson wrote:

> Thanks again Roger - I see about the slopes assignment.  I'm just learning 
> how to work with these spatial objects.  That should make them very nice to 
> work with actually, glad to realize this.
>
> I don't suppose there is any way to go from SpatialPolygons to a SpatialGrid 
> is there?  I'm having a problem using the shapes to  grid function in SAGA 
> via RSAGA geoprocessor, although most functions are working very well for me 
> otherwise.  There is a problem with the SAGA source code it seems. 
> Essentially I need a relatively rapid conversion from a shapefile to a 
> SpatialGridDataFrame (i.e. something not too computationally intensive) as I 
> am having to do some expensive formatting to make things work via RSAGA.

Did your subsequent message indicate that you had found a work-around?

If not, this reply copied from a recent posting to the statsgrass list may 
be relevant, using overlay() to assign grid cell centres to polygons:

library(sp)
library(maptools)
polyg <- readShapeSpatial(system.file("shapes/sids.shp",
  package="maptools")[1],  IDvar="FIPSNO",
  proj4string=CRS("+proj=longlat +ellps=clrk66"))
# the SpatialGrid is created here - if you already have one, just use it
grd <- Sobj_SpatialGrid(polyg)$SG
d <- overlay(grd, polyg)
# d contains the indices of the Polygons object into which each cell
# centre falls, or NA
summary(d)
names(polyg)
v <- polyg$SID74[d]
# use d to index the variable of interest
# if grd had been a SpatialGridDataFrame, just assign "[[<-" or "$<-",
# grd[["SID74"]] <- polyg[["SID74"]][d]
# or:
grdd <- SpatialGridDataFrame(grid=slot(grd, "grid"),
  data=data.frame(SID74=v), proj4string=CRS(proj4string(grd)))
summary(grdd)
image(grdd)
spplot(polyg, "SID74")

Roger

>
> Thanks again,
>
> Murray
>
>
> Roger Bivand wrote:
>>  On Fri, 15 Aug 2008, Murray Richardson wrote:
>> 
>> >  Oh maybe I've misunderstood the function.
>> > 
>> >  I am using unionSpatialPolygons to dissolve boundaries between adjacent 
>> >  polygons that have a similar attribute (I use cut to classify the ID 
>> >  into say, 10 different categories), i.e.:
>> > 
>> >  x<-readShapePoly("up1polys", IDvar="ID")
>> >  slopes<-slot(x[(1:length(slot(x, "polygons"))),3], "data")
>> >  breaks<-c(0,5,10,20,30,40,50,60,70,100)
>> >  ID <- cut(slopes[,1], breaks)
>> >  sptmp <- unionSpatialPolygons(x, ID)
>> >  newID<-data.frame(c(1:length(slot(sptmp, "polygons"))))
>> >  merged<-SpatialPolygonsDataFrame(sptmp, data=newID, match.ID=FALSE)
>> >  writePolyShape(merged, "up1merged", factor2char = TRUE, max_nchar=254)
>> > 
>> >  So I am just using unionSpatialPolygons as a dissolve tool but I would 
>> >  like non-adjacent polys within the same slope class to be separate 
>> >  polygons when I'm done.
>>
>>  The ID argument groups all, both contiguous and so merged, and
>>  non-contiguous, Polygon objects in the same Polygons object. So if you
>>  don't want them merged, only set the same ID values for the slivers to be
>>  joined, merge the polygons, then go on and do the other things.
>>
>>  Your assignment to slopes is very strange. What does names(x) say? If the
>>  third one is "slopes", why not just say x$slopes? The object does behave
>>  just like a data.frame, after all. I think that you need to do things
>>  strictly step by step, do the slivers first, then associate the correct
>>  attributes with the output polygons.
>>
>>  Roger
>> 
>> > 
>> >  Thanks again
>> > 
>> >  Murray
>> > 
>> > 
>> >  Roger Bivand wrote:
>> > >   On Fri, 15 Aug 2008, Murray Richardson wrote:
>> > > 
>> > > >   Thanks for this Roger.
>> > > > >   One other thing now...if I use unionSpatialPolygons as a 
>> > >  dissolve tool, >  is there a way to then explode distinct polygons 
>> > >  back to individual >  polygons? i.e. once all the slivers are gone and 
>> > >  I have polygons merged >  based on attributes, the result is a 
>> > >  multipart polygon for each ID I >  used for the merge, but I need them 
>> > >  back as separate polys.
>> > > 
>> > >  Each unique ID value should give a separate Polygons object, so look 
>> > >  at
>> > >   the ID vector before going into unionSpatialPolygons() to make sure 
>> > >   it
>> > >   does what you want.
>> > > 
>> > >   Roger
>> > > 
>> > > > >   Thanks
>> > > > >   Murray
>> > > > > > >   Roger Bivand wrote:
>> > > > >    On Wed, 13 Aug 2008, Murray Richardson wrote:
>> > > > > > > >    Hello again r.sig.geo list,
>> > > > > > >    Thanks Roger, for help on my previous question regarding > 
>> > > > iterating >    through a shapefile.
>> > > > > > >    I'm sure once I receive my copy of  "Applied Spatial Data 
>> > > > >   Analysis with >  R" I will find answers to simple questions 
>> > >  like this > >  on my own, but in >  the meantime....
>> > > > > > >    Is it possible to merge sliver polygons that fall below a 
>> > > > > certain >   threshold area with adjacent neighbours (e.g. 
>> > >  perhaps using > >  >  unionSpatialPolygons but without aggregating any 
>> > >  polygons?).  If a > >  >  sliver shares edges with more than one 
>> > >  polygon, it doesn't really > >  matter >  which one it merges with, 
>> > >  but if I had to choose a rule I > >  would have it >  merge with the 
>> > >  largest one.
>> > > > > > >   Not such a simple question ...
>> > > > > > >   Both the Polygon and Polygons objects in the 
>> > >  SpatialPolygons object > >  have
>> > > > >    "area" slots, with different roles. The Polygon objects have a 
>> > > > >    correct
>> > > > >    naive area in the geometry of the coordinates taken as planar. 
>> > >  The
>> > > > >    Polygons objects use the "gross" area of Polygon objects 
>> > >  belonging to
>> > > > >    them, but "only" to provide the plot order (plot from largest 
>> > > to > >   smallest
>> > > > >    to avoid over-painting).
>> > > > > > >    If you "trust" the area slot of the Polygons objects 
>> > >  (beware of hole
>> > > > >   Polygon objects), you can first find your candidate slivers by 
>> > > > >   retrieving
>> > > > >    the areas by:
>> > > > > > >    Polygons_areas <- sapply(slot(SPobj, "polygons"),
>> > > > >      function(x) slot(x, "area"))
>> > > > > > >   and set a cutoff. Then use poly2nb(SPobj, queen=FALSE) in 
>> > >  spdep to > >  find
>> > > > >   the neighbours (rook criterion). Next use the output object to 
>> > > > >   identify
>> > > > >    the largest neighbours of the sliver candidates, and build a 
>> > > > >    "new
>> > > > >    Polygons" ID vector. Finally, use unionSpatialPolygons(). I'm 
>> > > > >   assuming you
>> > > > >    wouldn't have asked if there was useful data in the slivers!
>> > > > > > >    Hope this helps,
>> > > > > > >    Roger
>> > > > > > > > >    Thanks in advance,
>> > > > > > >    Murray Richardson
>> > > > > > 
>>>>>>> _______________________________________________
>> > > > > >    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




More information about the R-sig-Geo mailing list