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

GRASS GRASS user
1

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?

Any suggestions for a computationally efficient way for such type of analyis?

Best regards,
Johannes

2

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

3

On 01/06/17 16:41, Moritz Lennert wrote:

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

[...]

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.

Instead of looping over the classes, you could also rasterize the layer 2 of buffers with v.to.rast [...] use=cat and then run r.stats on the two layers:

r.stats -c buffers,landuse

Add the results as another temporary table to the db and use it plus the correspondance table to fill your original table.

Moritz

4

Thank you Moritz,

your suggestion using r.stats and the rasterized layer 2 of buffers works really nice. It took me just a while to summarize and join all data and get them back into the db in the right format. However, I think I managed this task now. Thank you for your help!

I also think that either a correpondance table function and/or an enhancement of v.rast.stat to work with overlapping buffers would be a really nice enhancement.

···

On Thu, Jun 1, 2017 at 5:00 PM, Moritz Lennert <mlennert@club.worldonline.be> wrote:

On 01/06/17 16:41, Moritz Lennert wrote:

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

[…]

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.

Instead of looping over the classes, you could also rasterize the layer 2 of buffers with v.to.rast […] use=cat and then run r.stats on the two layers:

r.stats -c buffers,landuse

Add the results as another temporary table to the db and use it plus the correspondance table to fill your original table.

Moritz

5

Hi Johannes,

On 07/06/17 10:20, Johannes Radinger wrote:

Thank you Moritz,

your suggestion using r.stats and the rasterized layer 2 of buffers
works really nice. It took me just a while to summarize and join all
data and get them back into the db in the right format. However, I think
I managed this task now. Thank you for your help!

Would you be willing to share the final version of your approach ? This might be a nice seed for a module that would offer this function.

Moritz

6

Hi Moritz,
Hi All…

here my python code snippets that worked for me. I think there might be many ways to simplify the code…

import grass.script as grass
import sqlite3
import pandas

Def variable names

respondents_coord=“respondents_coord”
respondents_points=“respondents_points”
respondents_buffer=“respondents_buffer”
respondents_buffer_cat_2_layer=“respondents_buffer_cat_2_layer”
respondents_buffer_cat_2_raster=“respondents_buffer_cat_2_raster”
correspondance_table=“correspondance_table”
clc_summary_table=“clc_summary_table”
clc_summary_pivot=“clc_summary_pivot”

Import points from X-Y table

grass.run_command(“v.in.db”,
overwrite=True,
table=respondents_coord,
output=respondents_points)

Calculate buffers around points

grass.run_command(“v.buffer”,
overwrite=True,
flags=“t”,
input=respondents_points,
type=“point”,
output=respondents_buffer,
distance=50000)
grass.run_command(“v.category”,
overwrite=True,
input=respondents_buffer,
option=“add”,
layer=2,
out=respondents_buffer_cat_2_layer)

Get CORINE land use special classes per buffer

correspondance_table_values =
for line in grass.read_command(‘v.category’,
input_=respondents_buffer_cat_2_layer,
layer=‘1,2’,
option=‘print’).splitlines():
layers=line.split(‘|’)
l1 = layers[0].split(‘/’)
l2 = layers[1]
for cat in l1:
correspondance_table_values.append((cat, l2))

grass.run_command(“db.execute”,
sql=“CREATE TABLE correspondance_table (buffer_cat_1 INTEGER, buffer_cat_2 INTEGER)”)

conn = sqlite3.connect(‘/path/to/sqlite/sqlite.db’)
cur = conn.cursor()
cur.executemany(‘INSERT INTO correspondance_table (buffer_cat_1, buffer_cat_2) values (?,?)’,correspondance_table_values)
conn.commit()
conn.close()

grass.run_command(“v.to.rast”,
overwrite=True,
input=respondents_buffer_cat_2_layer,
layer=2,
output=respondents_buffer_cat_2_raster,
use=“cat”)

if clc_summary_pivot in grass.read_command(“db.tables”,flags=“p”).split(“\n”):
grass.run_command(“db.droptable”,
flags=“f”,
table=clc_summary_pivot)

grass.run_command(“db.execute”,
sql=“CREATE TABLE {} (buffer_cat_2 INTEGER, CLC_built_up INTEGER, CLC_arable INTEGER, CLC_permanent_crops INTEGER, CLC_grassland INTEGER, CLC_forest INTEGER, CLC_others INTEGER, CLC_intertidal_coastal INTEGER, CLC_water_bodies INTEGER, CLC_sea INTEGER)”.format(clc_summary_pivot))

clc_values =
for line in grass.read_command(“r.stats”,
flags=“cn”,
input=“{},{}”.format(respondents_buffer_cat_2_raster,“CLC_reclass@Corine_LandCover”)).splitlines():
clc_values.append((line.split(’ ‘)[0], line.split(’ ‘)[1],line.split(’ ')[2]))

df = pandas.DataFrame(clc_values, columns=[‘buffer_cat_2’, ‘CLC_class’, ‘count’])
CLC_class_cols = [“1”,“2”,“3”,“4”,“5”,“6”,“7”,“8”,“9”]
clc_summary_pivot_values = [tuple(x) for x in df.pivot(index=‘buffer_cat_2’, columns=‘CLC_class’, values=‘count’).reindex(columns=CLC_class_cols).fillna(0).astype(int).to_records(index=True)]

conn = sqlite3.connect(‘/path/to/sqlite/sqlite.db’)
cur = conn.cursor()
cur.executemany(‘INSERT INTO {} (buffer_cat_2, CLC_built_up, CLC_arable, CLC_permanent_crops, CLC_grassland, CLC_forest, CLC_others, CLC_intertidal_coastal, CLC_water_bodies, CLC_sea) values (?,?,?,?,?,?,?,?,?,?)’.format(clc_summary_pivot),clc_summary_pivot_values)
conn.commit()
conn.close()

grass.run_command(“db.execute”,
sql=“CREATE TABLE {} AS SELECT buffer_cat_1,SUM(CLC_built_up) AS CLC_built_up, SUM(CLC_arable) AS CLC_arable, SUM(CLC_permanent_crops) AS CLC_permanent_crops, SUM(CLC_grassland) AS CLC_grassland, SUM(CLC_forest) AS CLC_forest, SUM(CLC_others) AS CLC_others, SUM(CLC_intertidal_coastal) AS CLC_intertidal_coastal, SUM(CLC_water_bodies) AS CLC_water_bodies, SUM(CLC_sea) AS CLC_sea FROM {} AS A LEFT JOIN {} AS B ON A.buffer_cat_2=B.buffer_cat_2 GROUP BY buffer_cat_1”.format(clc_summary_table,correspondance_table,clc_summary_pivot))

Join information on land use back to original center points of buffers

grass.run_command(“v.db.join”,
map=respondents_points,
column=“id_INT”,
other_table=clc_summary_table,
other_column=“resp_id”)

cheers, Johannes

···

On Wed, Jun 14, 2017 at 11:10 AM, Moritz Lennert <mlennert@club.worldonline.be> wrote:

Hi Johannes,

On 07/06/17 10:20, Johannes Radinger wrote:

Thank you Moritz,

your suggestion using r.stats and the rasterized layer 2 of buffers
works really nice. It took me just a while to summarize and join all
data and get them back into the db in the right format. However, I think
I managed this task now. Thank you for your help!

Would you be willing to share the final version of your approach ? This might be a nice seed for a module that would offer this function.

Moritz