[R-sig-Geo] merge sliver polygons
Murray Richardson
murray.richardson at utoronto.ca
Wed Nov 12 22:09:35 CET 2008
Dear Roger or others,
You may recall this thread in which you helped me work out a region
merging routine for dealing with sliver polygons, or to simplify the
polygons resulting from an image analysis routine.
Based on your advice, I worked out an iterative routine that merges
adjacent polygons based on the highest degree of similarity. It repeats
until all merges are complete to satisfy a similarity threshold (i.e. no
more adjacent polygons whose attributes of interest are within a certain
similarity tolerance.
I realized however that I am having problems with polygon holes - a
problem you seemed to have foreshadowed towards the beginning of this
thread. The problem occurs when I call unionSpatialPolygons and it
appears that the hole polygons are automatically merged with the
adjacent polygons. If I understand correctly, a polygon that contains
a polygon hole will have a slot dedicated to the polygon hole, with
identical attributes to the non-hole polygon? Then when
unionSpatialPolygons is called they will be merged (at least in my case,
since the attributes are the same?). Perhaps you could confirm whether
this is correct, and if so, can you suggest a way around it?
Essentially, I do not want polygon holes to be included in the merging
process.
Thanks in advance!
Murray
Roger Bivand wrote:
> 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
>>> > > > R-sig-Geo at stat.math.ethz.ch
>>> > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>> > > > > >
>>
>
More information about the R-sig-Geo
mailing list