[GRASS-user] importing and cleaning overlapping polygons that are supposed to overlap

Markus Metz markus.metz.giswork at gmail.com
Thu May 16 00:34:47 PDT 2019


Hi Vero,

On Thu, May 16, 2019 at 3:11 AM Veronica Andreo <veroandreo at gmail.com>
wrote:
>
> Dear all,
>
> thanks for the answers...
>
> @Madi, I know, but that's how I got the data from a colleague using
SaTScan to get cluster sizes.
>
> So, these "clusters" are 3, they are represented by circular areas and 2
of them overlap, and it is just fine that they overlap. I just want one
centroid per area to be able to get one value per original area.
>
> I tested your solution, @Micha (thanks much for your time!), but it gives
me 4 values, and I need 3. Moreover, I can no longer recognize which
polygon is which since v.db.select for layer 3 reports cats from 1 to 4,
but d.vect shows something different (I'd assume 2/3 has become 4?). See
screenshot below.
>
> So, is it possible somehow to keep my 3 original polygons? And how can I
get raster values for my original 3 polygons in GRASS? Or is it that this
is not allowed by topology?

in this example, you get 4 areas and 3 different categories. The
overlapping part is represented as a small area with 2 categories. You can
restore the original polygons with e.g. v.extract type=area cat=1,
iterating over categotry values. If you want to get raster values for the
original polygons, you need to v.extract each category, rasterize it to
create a MASK, run r.univar and upload that value to the original vector.
v.rast.stats does not work here because it assumes that each area has only
one category value.

For point queries using centroids, I think Micha's solution is quite
elegant.

HTH,

Markus M
>
> best,
> Vero
>
>
>
> El mié., 15 may. 2019 a las 12:04, Micha Silver (<tsvibar at gmail.com>)
escribió:
>>
>>
>> On 15/05/2019 0:46, Veronica Andreo wrote:
>>
>> Hi all,
>>
>> I was working today with a very simple vector map which corresponds to
clusters (circular polygons) that overlap and it is just fine that they
overlap. So, i received a shapefile with 3 of these clusters. Two of them
overlaped. When I import them into GRASS with v.import I get an extra
centroid and area where 2 of the polygons overlap.
>>
>> Problem arises when I want to query a raster map with those polygons
since originally the attribute table contained only 3 polygons (which is
just fine). However, v.what.rast will only upload values for 2 of those
three polygons because it finds 2 centroids with the same category, AFAIU.
>>
>>
>> Here's how to get the raster values for all polygons (including those
which are overlaps). It involves the somewhat non-intuitive process of
creating another layer. I imported a simple polygon shapefile with three
overlapping areas. Here are the categories:
>>
>>
>> micha at TP480:~$ v.category polys opt=report
>> Layer/table: 1/polys
>> type       count        min        max
>> point          0          0          0
>> line           0          0          0
>> boundary       0          0          0
>> centroid      12          1          3
>> area           0          0          0
>> face           0          0          0
>> kernel         0          0          0
>> all           12          1          3
>> Layer: 2
>> type       count        min        max
>> point          0          0          0
>> line           0          0          0
>> boundary       0          0          0
>> centroid       4          2          3
>> area           0          0          0
>> face           0          0          0
>> kernel         0          0          0
>> all            4          2          3
>>
>>
>> So, as you found, you get two layers, with all original areas split up
at overlaps in layer 1, and just the overlap areas in layer 2.  Now add a
third layer, with new category values (option=add):
>>
>>
>> micha at TP480:~$ v.category polys option=add layer=3 out=polys_3layers
>>
>> Layer/table: 1/polys_3layers
>> type       count        min        max
>> point          0          0          0
>> line           0          0          0
>> boundary       0          0          0
>> centroid      12          1          3
>> area           0          0          0
>> face           0          0          0
>> kernel         0          0          0
>> all           12          1          3
>> Layer: 2
>> type       count        min        max
>> point          0          0          0
>> line           0          0          0
>> boundary       0          0          0
>> centroid       4          2          3
>> area           0          0          0
>> face           0          0          0
>> kernel         0          0          0
>> all            4          2          3
>> Layer: 3
>> type       count        min        max
>> point          0          0          0
>> line           0          0          0
>> boundary       0          0          0
>> centroid       7          1          7
>> area           0          0          0
>> face           0          0          0
>> kernel         0          0          0
>> all            7          1          7
>>
>>
>> This gives me three layers, and the new layer 3 has individual cat
values for each polygon. So I'm ready to add a database table and column to
that layer for the raster value, and then run v.what.rast:
>>
>>
>> micha at TP480:~$ v.db.addtable polys_3layers layer=3 columns="rast_value
DOUBLE"
>>
>> micha at TP480:~$ v.what.rast polys_3layers rast=dem_4m column=rast_value
layer=3 type=centroid
>> Reading features from vector map...
>> Update vector attributes...
>>  100%
>> v.what.rast complete. 7 records updated.
>> micha at TP480:~$ v.db.select polys_3layers layer=3
>> cat|rast_value
>> 1|488.3321
>> 2|492.7044
>> 3|481.2958
>> 4|498.173
>> 5|501.3336
>> 6|493.2202
>> 7|471.7223
>>
>> Does this help?
>>
>>
>>
>> I tried with
>>
>> v.clean input=clusters output=clusters_clean1
tool=break,rmdupl,rmsa,rmdac
>>
>> and
>>
>> v.clean input=clusters type=centroid output=clusters_clean2 tool=rmdupl
>>
>> but nothing seemed to do the kind of cleaning I wanted. I ended up using
brute force and removing the extra centroid manually in the GUI. That
helped with v.what.rast (all cats in the attr table were updated) but part
of the original and correct geometry was clearly gone. I am sure that's not
the right combination nor the right way either.
>>
>> In the wiki [0], I only found this piece of text: "If the input polygons
have logical errors.... You can investigate overlapping areas in the
imported vector with 'd.vect yourmap type=area layer=2' (only overlapping
areas have a category in layer 2 after import). Additionally you may show
the centroids of layer=2 to easier find tiny overlapping areas with 'd.vect
yourmap type=centroid layer=2'"
>>
>> However, it says nothing about how to proceed further as to keep
correctly overlapping polygons, each with its own centroid and remove the
duplicated ones that are generated when importing.
>>
>> Can someone please share the set of steps that should be followed in
these cases? Maybe it's a silly question, but I'm more a raster person so I
am very easily lost with vector issues.
>>
>> Thanks much in advance!
>> Vero
>>
>> [0] https://grasswiki.osgeo.org/wiki/Vector_topology_cleaning
>>
>> _______________________________________________
>> grass-user mailing list
>> grass-user at lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/grass-user
>>
>> --
>> Micha Silver
>> Ben Gurion Univ.
>> Sde Boker, Remote Sensing Lab
>> cell: +972-523-665918
>
> _______________________________________________
> grass-user mailing list
> grass-user at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/grass-user
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20190516/b9f157cb/attachment-0001.html>


More information about the grass-user mailing list