On Fri, 5 Sep 2008, Murray Richardson wrote:

Thanks Roger,

I have singletons with areas. I think the slivers will be quite obvious based on their areas relative to the rest of the polygon areas. I could even compute some perimeter to area index I supposed since they tend to be long a skinny.

So, my main uncertainty is still just how to work with the neighbour list to merge the sliver polygons with the LARGEST non-sliver polygon they share an edge with. i.e. to generate the new ids to use in unionSpatialPolygons. I can't seem to find a good example of how to use the nb object in such a way.

OK, so the Polygons are all singletons. This isn't tried:

is_a_sliver <- AREAS < very_small
library(spdep)
nb <- poly2nb(your_SPDF)
# check length(is_a_sliver) against length(nb)
crd <- card(nb)
IDs <- sapply(slot(your_SPDF, "polygons"), function(x) slot(x, "ID"))
for (i in seq(along=nb)) {
  if (is_a_sliver[i] && crd[i] > 0) {
    nbi <- nb[[i]]
    max_area <- nbi[which.max(AREAS[nbi])]
    IDs[i] <- IDs[max_area]
  }
}

then use IDs in unionSpatialPolygons()

I repeat, untried, but something like this should work.

Roger

PS. If you have the perimeter, you could try it as a sanity check if area alone give spurious results.


cheers

Murray


Roger Bivand wrote:
 On Fri, 5 Sep 2008, Murray Richardson wrote:

>  Hi Roger
> > I'm just getting around to trying this out. I must say, I'm not clear on > how to work with the neighbour list object to accomplish this. > Specifically, how would I identify the largest neighbours of the sliver > candidates? > > Furthermore, do I first subset the polygons according to the area > threshold, and then construct the neighbour list for only these > candidate slivers? I'm not clear on how to identify neighbours for a > subset of polygons.

 My idea (IIRC!) was to look at a distribution of the areas of the Polygon
 objects in the dataset. A good deal will depend on whether the slivers are
 represented as separate Polygons objects (note the _s_), or whether they
 are Polygon objects belonging to Polygons objects. So nimbleness will be
 needed!

 library(maptools)
 xx <- readShapeSpatial(system.file("shapes/sids.shp",
 package="maptools")[1],
       IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))
 pls <- slot(xx, "polygons")
 AREAS <- sapply(pls, function(x) sapply(slot(x, "Polygons"), function(y)
       slot(y, "area")))
 summary(sapply(AREAS, length))
 summary(unlist(AREAS))

 is not a very good example, because the "areas" are planar in non-planar
 coordinates, but I hope you see the point. If you have a clear break in
 the distribution, you can find out which they are, and where they are. If
 all the Polygons objects (the ones in the "polygons" slot of your
 SpatialPolygons object, pls) are singletons, things get much easier,
 because the merging will be between Polygons objects, provided directly in
 maptools.

 This should be a start - once you have singleton Polygons objects with
 areas as attributes, we can continue.

 Hope this helps,

 Roger

> > Sorry to bother you with this again. It will be very useful if it works > as I would also like to use it for region merging based on other > attributes. > > Still waiting on that book...!

 PS. Edzer, Virgilio and I did see the only copy in Europe in Dortmund at
 useR! three weeks ago. It turns out that it was printed in the US, and is
 at or awaited at the publisher's warehouse for dispatch. Things seem to
 take time in the "real" world!

> > 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
> > >   [email protected]
> > >   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: [EMAIL PROTECTED]

_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Reply via email to