[GRASS-user] zonal statistics for multiple areas

GRASS GRASS user
1

Hello to everybody,
I hope my question is appropriate for this mail list. I have created 5 different buffer layers (d=100,200,300,400,500 m) from a point vector. Now i want to compute some zonal statistics using different rasters as input. I know how to do that separately, but i was wondering how could i automate this procedure using a python script. I know python’s basics, but i was not able to do so following guidelines that i found on the internet, e.g GRASS tutorial. Thank you in advance, A

2

Hi Alessandro,

You could try v.rast.bufferstats [1], if you do not have hundrets of thousands of points….

If the buffers don`t overlap, you could use v.rast.stats [2] from the upcoming GRASS 7.6 release, that will allow multiple raster input…

Otherwise have a look at the general introduction to GRASS and Python [3].

Cheers

Stefan

1: https://grass.osgeo.org/grass74/manuals/addons/v.rast.bufferstats.html

2: https://grass.osgeo.org/grass76/manuals/v.rast.stats.html

3: https://grasswiki.osgeo.org/wiki/GRASS_and_Python

···

Hello to everybody,

I hope my question is appropriate for this mail list. I have created 5 different buffer layers (d=100,200,300,400,500 m) from a point vector. Now i want to compute some zonal statistics using different rasters as input. I know how to do that separately, but i was wondering how could i automate this procedure using a python script. I know python’s basics, but i was not able to do so following guidelines that i found on the internet, e.g GRASS tutorial. Thank you in advance, A

3

Some further python specific tips that might help.

To get a list of rasters that match some pattern (say you have many rasters like “band1”, “band2”,…), then run v.rast.stats on each of them in a loop to get zonal stats in the vector “buffers”:

# This must be run from within a grass session

import grass.script as gscript

rast_list = gscript.read_command("g.list", type = "rast", pattern = "band*").strip().split("\n")

# Loop thru the rasters, run v.rast.stats, and specify the column prefix using the raster name.

methods = ("minimum", "maximum", "average", "sum") # add whatever stats you need

for r in rast_list:

gscript.run_command("v.rast.stats",

map_ = "buffers", raster=r, column_prefix=r, method=methods)

Be aware that if you have many rasters, you’ll end up with many,many columns in the “buffers” vector

···

On 10/29/18 9:54 PM, Stefan Blumentrath wrote:

Hi Alessandro,

You could try v.rast.bufferstats [1], if you do not have hundrets of thousands of points….

If the buffers don`t overlap, you could use v.rast.stats [2] from the upcoming GRASS 7.6 release, that will allow multiple raster input…

Otherwise have a look at the general introduction to GRASS and Python [3].

Cheers

Stefan

1: https://grass.osgeo.org/grass74/manuals/addons/v.rast.bufferstats.html

2: https://grass.osgeo.org/grass76/manuals/v.rast.stats.html

3: https://grasswiki.osgeo.org/wiki/GRASS_and_Python

From: grass-user grass-user-bounces@lists.osgeo.org On Behalf Of Alessandro Sebastiani
Sent: mandag 29. oktober 2018 16:06
To: grass-user@lists.osgeo.org
Subject: [GRASS-user] zonal statistics for multiple areas

Hello to everybody,

I hope my question is appropriate for this mail list. I have created 5 different buffer layers (d=100,200,300,400,500 m) from a point vector. Now i want to compute some zonal statistics using different rasters as input. I know how to do that separately, but i was wondering how could i automate this procedure using a python script. I know python’s basics, but i was not able to do so following guidelines that i found on the internet, e.g GRASS tutorial. Thank you in advance, A

_______________________________________________
grass-user mailing list
[grass-user@lists.osgeo.org](mailto:grass-user@lists.osgeo.org)
[https://lists.osgeo.org/mailman/listinfo/grass-user](https://lists.osgeo.org/mailman/listinfo/grass-user)

-- 
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918
4

thank you all for your suggestions.
Dear Nikos, I have some 0-1 rasters that represents presence-absence of different land covers. Each raster’s resolution is 10x10m, thus i want to compute the sum in order to obtain the surface covered by each land cover within my buffers. As regards to the buffer, i want to use the 0-100, 0-200, 0-300, 0-400 and 0-500 for computing stats. I can do it separately with the module v.rast.stats. I just don’t know how to loop it in a python script
Dear Stefan, thank you, i guess that using multiple rasters as input would be a good solution. however i would really like to loop it using a python script since that would be useful in other situation that i’m going to face in the next days
Dear Micha thank you, i’ll try ASAP!

Il giorno lun 29 ott 2018 alle ore 21:07 Micha Silver <tsvibar@gmail.com> ha scritto:

Some further python specific tips that might help.

To get a list of rasters that match some pattern (say you have many rasters like “band1”, “band2”,…), then run v.rast.stats on each of them in a loop to get zonal stats in the vector “buffers”:

# This must be run from within a grass session

import grass.script as gscript

rast_list = gscript.read_command("g.list", type = "rast", pattern = "band*").strip().split("\n")

# Loop thru the rasters, run v.rast.stats, and specify the column prefix using the raster name.

methods = ("minimum", "maximum", "average", "sum") # add whatever stats you need

for r in rast_list:

gscript.run_command("v.rast.stats",

map_ = "buffers", raster=r, column_prefix=r, method=methods)

Be aware that if you have many rasters, you’ll end up with many,many columns in the “buffers” vector

On 10/29/18 9:54 PM, Stefan Blumentrath wrote:

Hi Alessandro,

You could try v.rast.bufferstats [1], if you do not have hundrets of thousands of points….

If the buffers don`t overlap, you could use v.rast.stats [2] from the upcoming GRASS 7.6 release, that will allow multiple raster input…

Otherwise have a look at the general introduction to GRASS and Python [3].

Cheers

Stefan

1: https://grass.osgeo.org/grass74/manuals/addons/v.rast.bufferstats.html

2: https://grass.osgeo.org/grass76/manuals/v.rast.stats.html

3: https://grasswiki.osgeo.org/wiki/GRASS_and_Python

From: grass-user grass-user-bounces@lists.osgeo.org On Behalf Of Alessandro Sebastiani
Sent: mandag 29. oktober 2018 16:06
To: grass-user@lists.osgeo.org
Subject: [GRASS-user] zonal statistics for multiple areas

Hello to everybody,

I hope my question is appropriate for this mail list. I have created 5 different buffer layers (d=100,200,300,400,500 m) from a point vector. Now i want to compute some zonal statistics using different rasters as input. I know how to do that separately, but i was wondering how could i automate this procedure using a python script. I know python’s basics, but i was not able to do so following guidelines that i found on the internet, e.g GRASS tutorial. Thank you in advance, A

_______________________________________________
grass-user mailing list
[grass-user@lists.osgeo.org](mailto:grass-user@lists.osgeo.org)
[https://lists.osgeo.org/mailman/listinfo/grass-user](https://lists.osgeo.org/mailman/listinfo/grass-user)

-- 
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918
5

Resending, as I forgot to Cc the list.

* Alessandro Sebastiani <alessandro.sebastiani@uniroma1.it> [2018-10-29 16:05:54 +0100]:

Hello to everybody,o

Hello Alessandro,

I hope my question is appropriate for this mail list. I have created 5
different buffer layers (d=100,200,300,400,500 m) from a point vector. Now
i want to compute some zonal statistics using different rasters as input.

You are in the right place. And the first guess is that what you need to
do is not difficult using GRASS GIS.

Which rasters exactly? Would you want to use as base maps each of the
'0 to 100 m', '100 to 200 m' and so on "buffers" or '0 to 100 m', '0
to 200 m' and so on up to '0 to 500 m'? And then, the 'other' maps
would be the 'cover' map, each time?

If you share exactly what you need to do, i.e. in form of an algorithm
using pseudo-syntax, it will be easier to understand and potentially
suggest ways to get this done.

For example:

```
for each Buffer in Buffers:
    compute zonal statistics using base=Buffer and cover=RasterX ```

For example, you might be able to get part of this done using raster maps
only and `r.mapcalc` expressions, without the need to construct many loops and
require an extra step in the end to patch maps together.

I know how to do that separately, but i was wondering how could i automate
this procedure using a python script.

If you write down the way to do it separately, then you are just few
steps away in looping over your input rasters. The key point is,
obviously, to express what is to be done in an unambiguous way.

I know python's basics, but i was not
able to do so following guidelines that i found on the internet, e.g GRASS
tutorial.

Please share which guidelines, i.e. which GRASS GIS tutorial you refer
to.

Nikos

6

* Alessandro Sebastiani <alessandro.sebastiani@uniroma1.it> [2018-10-29 21:32:51 +0100]:

thank you all for your suggestions.
Dear Nikos, I have some 0-1 rasters that represents presence-absence of
different land covers. Each raster's resolution is 10x10m, thus i want to
compute the sum in order to obtain the surface covered by each land cover
within my buffers. As regards to the buffer, i want to use the 0-100,
0-200, 0-300, 0-400 and 0-500 for computing stats. I can do it separately
with the module v.rast.stats. I just don't know how to loop it in a python
script

Dear Alessandro,

I don't have time to write back as I would like. Instead, I will send
you off-list some unpublished script. With a bit of reading, there, you
might identify some ideas that will eventually help in what you are
doing.

Kind regards, Nikos

7

Thanks a lot for your help!

Il giorno mar 30 ott 2018 alle ore 11:49 Nikos Alexandris <nik@nikosalexandris.net> ha scritto:

thank you all for your suggestions.
Dear Nikos, I have some 0-1 rasters that represents presence-absence of
different land covers. Each raster’s resolution is 10x10m, thus i want to
compute the sum in order to obtain the surface covered by each land cover
within my buffers. As regards to the buffer, i want to use the 0-100,
0-200, 0-300, 0-400 and 0-500 for computing stats. I can do it separately
with the module v.rast.stats. I just don’t know how to loop it in a python
script

Dear Alessandro,

I don’t have time to write back as I would like. Instead, I will send
you off-list some unpublished script. With a bit of reading, there, you
might identify some ideas that will eventually help in what you are
doing.

Kind regards, Nikos

8

Dear all, i was able to compute zonal statistics for multiple areas using this script. I used

import grass.script as gscript
v_list = gscript.read_command(“g.list” , type = “vect”).strip().split(“\n”)

i have 5 vectors in my list

method = (“sum” , “average”)
for v in v_list:
gscript.run_command(“v.rast.stats”, map_=v, raster = “LC1”, column_prefix = “LC1”, method = methods)

Now i have like 4 rasters that i would like to use as input to compute my statistics. Is there a way to loop this process for rasters and vector at the same time? i tried creating a list for rasters and using the name list as raster input but it seems not to work! cheers, Alessandro

Il giorno mar 30 ott 2018 alle ore 11:52 Alessandro Sebastiani <alessandro.sebastiani@uniroma1.it> ha scritto:

Thanks a lot for your help!

Il giorno mar 30 ott 2018 alle ore 11:49 Nikos Alexandris <nik@nikosalexandris.net> ha scritto:

thank you all for your suggestions.
Dear Nikos, I have some 0-1 rasters that represents presence-absence of
different land covers. Each raster’s resolution is 10x10m, thus i want to
compute the sum in order to obtain the surface covered by each land cover
within my buffers. As regards to the buffer, i want to use the 0-100,
0-200, 0-300, 0-400 and 0-500 for computing stats. I can do it separately
with the module v.rast.stats. I just don’t know how to loop it in a python
script

Dear Alessandro,

I don’t have time to write back as I would like. Instead, I will send
you off-list some unpublished script. With a bit of reading, there, you
might identify some ideas that will eventually help in what you are
doing.

Kind regards, Nikos

9

On 10/31/18 11:59 AM, Alessandro Sebastiani wrote:

Dear all, i was able to compute zonal statistics for multiple areas using this script. I used
> import grass.script as gscript
> v_list = gscript.read_command("g.list" , type = "vect").strip().split("\n")
# i have 5 vectors in my list
> method = ("sum" , "average")
> for v in v_list:
> gscript.run_command("v.rast.stats", map_=v, raster = "LC1", column_prefix = "LC1", method = methods)

Now i have like 4 rasters that i would like to use as input to compute my statistics. Is there a way to loop this process for rasters and vector at the same time? i tried creating a list for rasters and using the name list as raster input but it seems not to work! cheers, Alessandro

What exactly did you try? There should not be any problem to do nested loops. i.e.:

import grass.script as gscript

v_list = gscript.read_command("g.list" , type = "vect").strip().split("\n")
# i have 5 vectors in my list
r_list = gscript.read_command("g.list" , type = "rast", pattern="LC*").strip().split("\n")
method = ("sum" , "average")
for v in v_list:
for r in r_list:
gscript.run_command("v.rast.stats",
map_=v, raster = r, column_prefix = r, method = methods)

Il giorno mar 30 ott 2018 alle ore 11:52 Alessandro Sebastiani <alessandro.sebastiani@uniroma1.it <mailto:alessandro.sebastiani@uniroma1.it>> ha scritto:

    Thanks a lot for your help!

    Il giorno mar 30 ott 2018 alle ore 11:49 Nikos Alexandris
    <nik@nikosalexandris.net <mailto:nik@nikosalexandris.net>> ha scritto:

        * Alessandro Sebastiani <alessandro.sebastiani@uniroma1.it
        <mailto:alessandro.sebastiani@uniroma1.it>> [2018-10-29
        21:32:51 +0100]:

        >thank you all for your suggestions.
        >Dear Nikos, I have some 0-1 rasters that represents
        presence-absence of
        >different land covers. Each raster's resolution is 10x10m,
        thus i want to
        >compute the sum in order to obtain the surface covered by
        each land cover
        >within my buffers. As regards to the buffer, i want to use
        the 0-100,
        >0-200, 0-300, 0-400 and 0-500 for computing stats. I can do
        it separately
        >with the module v.rast.stats. I just don't know how to loop
        it in a python
        >script

        Dear Alessandro,

        I don't have time to write back as I would like. Instead, I
        will send
        you off-list some unpublished script. With a bit of reading,
        there, you
        might identify some ideas that will eventually help in what
        you are
        doing.

        Kind regards, Nikos

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

--
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918

10

I was trying with nested loops and it seems not to be working. Now I tried again (more carefully i have to say) and it worked! Thanks to everybody!

Il giorno mer 31 ott 2018 alle ore 21:55 Micha Silver <tsvibar@gmail.com> ha scritto:

On 10/31/18 11:59 AM, Alessandro Sebastiani wrote:

Dear all, i was able to compute zonal statistics for multiple areas
using this script. I used

import grass.script as gscript
v_list = gscript.read_command(“g.list” , type =
“vect”).strip().split(“\n”)

i have 5 vectors in my list

method = (“sum” , “average”)
for v in v_list:
gscript.run_command(“v.rast.stats”, map_=v, raster = “LC1”,
column_prefix = “LC1”, method = methods)

Now i have like 4 rasters that i would like to use as input to compute
my statistics. Is there a way to loop this process for rasters and
vector at the same time? i tried creating a list for rasters and using
the name list as raster input but it seems not to work! cheers, Alessandro

What exactly did you try? There should not be any problem to do nested
loops. i.e.:

import grass.script as gscript

v_list = gscript.read_command(“g.list” , type = “vect”).strip().split(“\n”)

i have 5 vectors in my list

r_list = gscript.read_command(“g.list” , type = “rast”,
pattern=“LC*”).strip().split(“\n”)
method = (“sum” , “average”)
for v in v_list:
for r in r_list:
gscript.run_command(“v.rast.stats”,
map_=v, raster = r, column_prefix = r, method = methods)

Il giorno mar 30 ott 2018 alle ore 11:52 Alessandro Sebastiani
<alessandro.sebastiani@uniroma1.it
mailto:[alessandro.sebastiani@uniroma1.it](mailto:alessandro.sebastiani@uniroma1.it)> ha scritto:

Thanks a lot for your help!

Il giorno mar 30 ott 2018 alle ore 11:49 Nikos Alexandris
<nik@nikosalexandris.net mailto:[nik@nikosalexandris.net](mailto:nik@nikosalexandris.net)> ha scritto:

thank you all for your suggestions.
Dear Nikos, I have some 0-1 rasters that represents
presence-absence of
different land covers. Each raster’s resolution is 10x10m,
thus i want to
compute the sum in order to obtain the surface covered by
each land cover
within my buffers. As regards to the buffer, i want to use
the 0-100,
0-200, 0-300, 0-400 and 0-500 for computing stats. I can do
it separately
with the module v.rast.stats. I just don’t know how to loop
it in a python
script

Dear Alessandro,

I don’t have time to write back as I would like. Instead, I
will send
you off-list some unpublished script. With a bit of reading,
there, you
might identify some ideas that will eventually help in what
you are
doing.

Kind regards, Nikos

grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

–
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918