[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