<div dir="ltr">Hi Vero,<br><div><br>On Thu, May 16, 2019 at 3:11 AM Veronica Andreo <<a href="mailto:veroandreo@gmail.com">veroandreo@gmail.com</a>> wrote:<br>><br>> Dear all,<br>><br>> thanks for the answers...<br>><br>> @Madi, I know, but that's how I got the data from a colleague using SaTScan to get cluster sizes.<br>><br>> 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.<br>><br>> 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.<br>><br>> 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?<br></div><div><br></div><div>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.</div><div><br></div><div>For point queries using centroids, I think Micha's solution is quite elegant.<br></div><div><br></div><div>HTH,</div><div><br></div><div>Markus M<br></div><div>><br>> best,<br>> Vero<br>><br>><br>><br>> El mié., 15 may. 2019 a las 12:04, Micha Silver (<<a href="mailto:tsvibar@gmail.com">tsvibar@gmail.com</a>>) escribió:<br>>><br>>><br>>> On 15/05/2019 0:46, Veronica Andreo wrote:<br>>><br>>> Hi all,<br>>><br>>> 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.<br>>><br>>> 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.<br>>><br>>><br>>> 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:<br>>><br>>><br>>> micha@TP480:~$ v.category polys opt=report<br>>> Layer/table: 1/polys<br>>> type       count        min        max<br>>> point          0          0          0<br>>> line           0          0          0<br>>> boundary       0          0          0<br>>> centroid      12          1          3<br>>> area           0          0          0<br>>> face           0          0          0<br>>> kernel         0          0          0<br>>> all           12          1          3<br>>> Layer: 2<br>>> type       count        min        max<br>>> point          0          0          0<br>>> line           0          0          0<br>>> boundary       0          0          0<br>>> centroid       4          2          3<br>>> area           0          0          0<br>>> face           0          0          0<br>>> kernel         0          0          0<br>>> all            4          2          3<br>>><br>>><br>>> 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):<br>>><br>>><br>>> micha@TP480:~$ v.category polys option=add layer=3 out=polys_3layers<br>>><br>>> Layer/table: 1/polys_3layers<br>>> type       count        min        max<br>>> point          0          0          0<br>>> line           0          0          0<br>>> boundary       0          0          0<br>>> centroid      12          1          3<br>>> area           0          0          0<br>>> face           0          0          0<br>>> kernel         0          0          0<br>>> all           12          1          3<br>>> Layer: 2<br>>> type       count        min        max<br>>> point          0          0          0<br>>> line           0          0          0<br>>> boundary       0          0          0<br>>> centroid       4          2          3<br>>> area           0          0          0<br>>> face           0          0          0<br>>> kernel         0          0          0<br>>> all            4          2          3<br>>> Layer: 3<br>>> type       count        min        max<br>>> point          0          0          0<br>>> line           0          0          0<br>>> boundary       0          0          0<br>>> centroid       7          1          7<br>>> area           0          0          0<br>>> face           0          0          0<br>>> kernel         0          0          0<br>>> all            7          1          7<br>>><br>>><br>>> 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:<br>>><br>>><br>>> micha@TP480:~$ v.db.addtable polys_3layers layer=3 columns="rast_value DOUBLE"<br>>><br>>> micha@TP480:~$ v.what.rast polys_3layers rast=dem_4m column=rast_value layer=3 type=centroid<br>>> Reading features from vector map...<br>>> Update vector attributes...<br>>>  100%<br>>> v.what.rast complete. 7 records updated.<br>>> micha@TP480:~$ v.db.select polys_3layers layer=3<br>>> cat|rast_value<br>>> 1|488.3321<br>>> 2|492.7044<br>>> 3|481.2958<br>>> 4|498.173<br>>> 5|501.3336<br>>> 6|493.2202<br>>> 7|471.7223<br>>><br>>> Does this help?<br>>><br>>><br>>><br>>> I tried with<br>>><br>>> v.clean input=clusters output=clusters_clean1 tool=break,rmdupl,rmsa,rmdac<br>>><br>>> and<br>>><br>>> v.clean input=clusters type=centroid output=clusters_clean2 tool=rmdupl<br>>><br>>> 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.<br>>><br>>> 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'"<br>>><br>>> 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.<br>>><br>>> 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.<br>>><br>>> Thanks much in advance!<br>>> Vero<br>>><br>>> [0] <a href="https://grasswiki.osgeo.org/wiki/Vector_topology_cleaning">https://grasswiki.osgeo.org/wiki/Vector_topology_cleaning</a><br>>><br>>> _______________________________________________<br>>> grass-user mailing list<br>>> <a href="mailto:grass-user@lists.osgeo.org">grass-user@lists.osgeo.org</a><br>>> <a href="https://lists.osgeo.org/mailman/listinfo/grass-user">https://lists.osgeo.org/mailman/listinfo/grass-user</a><br>>><br>>> -- <br>>> Micha Silver<br>>> Ben Gurion Univ.<br>>> Sde Boker, Remote Sensing Lab<br>>> cell: +972-523-665918<br>><br>> _______________________________________________<br>> grass-user mailing list<br>> <a href="mailto:grass-user@lists.osgeo.org">grass-user@lists.osgeo.org</a><br>> <a href="https://lists.osgeo.org/mailman/listinfo/grass-user">https://lists.osgeo.org/mailman/listinfo/grass-user</a></div></div>