[GRASS-user] Single centroid/point for multiple polygons sharing the same category

[Moritz Lennert]

On 19/06/17 13:02, Johannes Radinger wrote:
> Hi GRASS users,
>
> I am working on a polygon map of France that shows the single areas
> belonging to a specific postal code. I want to extract for each area the
> centroid (using v.extract), transform it to a point map (using v.type)
> and finally add the X and Y coordinates to the attribute table for each
> postal code region.
>
> However, now I have the problem that some "postal code areas" have more
> than one centroid, i.e. two or more closeby but not connected polygons
> belong the same postal region and share a common category and entry in
> the attribute table. So for these cases I get two or more points
> (sharing the same cat) and of course I cannot calculate a single XY pair
> for that postal code region.
>
> Consequently, I need to find a method to get centroids of all the areas
> sharing a common cat or I need to calculate the mean position of points
> afterwards (using v.centerpoint). As far as I understood v.centerpoint
> all points are used to calculate a mean point but not stratified for
> points sharing the same cat!? FYI: I have around 2000 postal code areas
> of which 200 consists of more than a single polygon.
>
> Any suggestions?


Here are two examples using the urbanarea map of the NC dataset.

1) Use average centroid coordinates

v.category urbanarea op=add layer=2 out=myareas
v.db.addtable myareas layer=2 col="lay1cat integer, x double precision, 
y double precision"
v.to.db myareas layer=2 op=query query_layer=1 query_col=cat col=lay1cat
v.to.db myareas layer=2 op=coor col=x,y
v.db.addcolumn myareas col="x double precision, y double precision"
db.execute sql="update myareas set x=(select avg(x) from myareas_2 where 
lay1cat=myareas.cat group by lay1cat), y=(select avg(y) from myareas_2 
where lay1cat=myareas.cat group by lay1cat)"

2) Use centroid coordinates of the largest polygon

v.category urbanarea op=add layer=2 out=myareas
v.db.addtable myareas layer=2 col="lay1cat integer, area double precision"
v.to.db myareas layer=2 op=query query_layer=1 query_col=cat col=lay1cat
v.to.db myareas layer=2 op=area col=area
v.extract myareas layer=2 cat=$(db.select -c sql="select cat from 
myareas_2 group by lay1cat having area = max(area)" | awk '{printf"%s,", 
$1}') out=my_unique_areas type=area

#change layers to get default layer 1 in final output (not really 
necessary, but makes it nicer ot use afterwards)
v.category my_unique_areas op=chlayer layer=2,1 out=final_areas
v.db.connect -d final_areas layer=2
v.db.connect final_areas layer=1 table=final_areas

v.db.addcolumn final_areas col="x double precision, y double precision"
v.to.db final_areas op=coor col=x,y
v.db.join map=myareas layer=1 col=cat otable=final_areas ocol=lay1cat 
subset_col=x,y

Moritz