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
___
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user