[GRASS-user] Obtain land use within 50 km buffers around vector points

[Moritz Lennert]

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