Re: [GRASS-user] Obtain land use within 50 km buffers around vector points
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=5)
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 <
mlenn...@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
>
>
___
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Obtain land use within 50 km buffers around vector points
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 ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Obtain land use within 50 km buffers around vector points
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 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 > > ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Obtain land use within 50 km buffers around vector points
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 ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Obtain land use within 50 km buffers around vector points
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
[GRASS-user] Obtain land use within 50 km buffers around vector points
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 ___ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user
