On 01/06/17 13:07, Johannes Radinger wrote:
Hi,
I have a vector map of 50 km buffers surrounding 4000 points in entire
Europe. The single buffers where calculated using v.buffer with the -t
flag. Thus, many of my buffers overlap (i.e. geometries of buffers are
not merged but split up). Now, I want to obtain for each 50 km buffer
the underlying land use (from a CORINE land use map). For example, I'd
like to know the relative share of "forest" within each buffer and I
want to update this information in the attribute table of the vector
buffer map.
So what would be the GRASS way for such type of analysis involving
overlapping buffers? Do I need to loop over the single buffers, extract
each buffer, and run something like v.rast.stats using maps for the
single land use classes as raster input?
Ideally one should be able to run v.rast.stats directly on the
overlapping buffers, and get the results by cat number, but this does
not work (cf http://trac.osgeo.org/grass/ticket/3300) because of the way
v.rast.stats works.
Currently, the best way I see is to create a correspondance table
between the pieces and the original cat values. You can get that by running
v.category buffers op=add layer=2 out=buffers_2_layers
v.category buffers_2_layers layers=1,2 option=print
and massaging the results a bit to create a new database table out that
(would be a nice enhancement to v.buffer if it created a second layer
and this lookup table as an option). Here's a very raw version of how
this could look like in Python. The correspondance_table would then have
to be written to the database:
correspondance_table = []
for line in g.read_command('v.category', input_='buffers_2layers',
layer='1,2', option='print').splitlines():
layers=line.split('|')... l1 = layers[0].split('/')
l2 = layers[1]
for cat in l1:
correspondance_table.append((cat, l2))
Then just run v.rast.stats once per class (each class raster can just be
a reclass of the original with classnum = 1\n* = NULL") on layer 2 of
the entire vector buffer map. Thus you will get the stats per piece of
buffer.
To calculate total values per buffer you can then use using the
correspondance table, e.g. something like this (warning: untested) to
get the total number of pixels for a given class in your entire 50km
buffers:
db.exectute sql="UPDATE buffer_layer_1_table SET nb_class_X = t.class_X
FROM (SELECT l1.cat, sum(l2.class_X) as class_X FROM
buffer_layer_1_table l1 JOIN correspondance_table c ON (l1.cat =
c.l1_cat) JOIN buffer_layer_2_table l2 ON (l2.cat = c.l1_cat)
GROUP BY l1.cat) t WHERE buffer_layer_1_table.cat = t.cat"
It would be a nice add-on to have a module that calculates shares of
raster classes per polygon, with either the user providing an optional
correspondance_table and layers info if the polygons are overlapping or
the option to just indicate that polygons are overlapping and the module
would do all this internally.
Moritz
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user
