[Qgis-user] Zonal Statistics multiple classe

Nikos Alexandris nik at nikosalexandris.net
Mon Feb 12 11:52:00 PST 2018


(Cc in the list, someone might pick this up and see why you can't go
through)

Some examples in the end using r.stats.zonal. I think they are useful
for your working case. I am sorry that I don't include specific examples for
v.rast.stats. But they don't differ much.

You could convert your vector polygons in a raster map. Think about the
resolution of the computational region before doing so.


* Nicolas Cadieux <nicolas.cadieux at archeotec.ca> [2018-02-12 10:04:33 -0500]:

>Hi Nikos,
>I am having problems replying to the user list.  The Email comes back 
>to me as unsent (I need to contact the list admin). I am therefore 
>replying using your email address.  Can you confirm that you got this 
>email.
>Sorry if you received multiple copies.
>Nicolas
>
>
>Thanks for your help,
>
>I started playing with the v.rast.stats but with difficulty.  I Will 
>work on it more today.  I figured it may have to do with the -r flag.
>
>How can I convert the int16 to a category format in grass if I was not 
>using the -r flag?  I am trying to understand the data type in Grass 
>for the category format?  I was hoping int16 would cover it.

I am sure you know all of the following, but just to be on the same ground:

GRASS understands all integers (automatically). If your raster map is of
type integer, you don't need to do anything about it. Just take care to
set the computational region using `g.region`.  This is something I
(still) frequently forget. The region is important for when using
`r.report`! See also: https://grasswiki.osgeo.org/wiki/GRASS_raster_semantics.

With the `r.category` module one can associate category values (in which
case "integer raster values") and category labels. See the manual for
details and examples.


>“If I understand correctly, you want to get statistics per polygon, not
>>per class. I.e., min, max, sum and the rest, for each class inside
>>each vector polygon. “
>
>Basically I am looking to getting the total cover and percent cover 
>for each class (the classes are crop data), for each polygons 
>(drainage area).  So basically, how many square km of corn, forest, 
>wheat... do I have in each drainage area polygons? Does that make more 
>sense for you?

The modules mentioned are a perfect fit for this task.

All you need is to import in GRASS your raster map with (integer) land
cover categories and your vector map.  The land cover map will be the
"cover". The vector map, in your case, will be the base.

>I was hoping the LecoS plugin would work but I think It need to work 
>with one crop type at a time.  Working on that too!  Patience is a 
>virtue... I guess!
>
>Thanks for your help
>Nicolas

Thanks to the GRASS community.

Compute and report (and even output as an ASCII file) zonal statistics per land
cover category:

```
r.stats.zonal method=average base=landcover cover=temperature output=average_temperature_per_landcover -r
r.report average_temperature_per_landcover -nh units=me,c,p --q output=zonal_statistics_report
```

Sample report where the 1, 2, 10, .., 50 are land cover classes:

+-----------------------------------------------------------------------------+
|             Category Information               |       square|   cell|   %  |
| #|description                                  |       meters|  count| cover|
|-----------------------------------------------------------------------------|
| 1|28.447180. . . . . . . . . . . . . . . . . . |   25,967,465|  28853|  0.80|
| 2|38.517164. . . . . . . . . . . . . . . . . . |    3,245,371|   3606|  0.10|
|10|30.085036. . . . . . . . . . . . . . . . . . |  289,801,877| 322005|  8.94|
|15|31.038448. . . . . . . . . . . . . . . . . . |   16,441,051|  18268|  0.51|
|20|26.774910. . . . . . . . . . . . . . . . . . |2,630,504,294|2922809| 81.11|
|25|29.068074. . . . . . . . . . . . . . . . . . |   28,675,540|  31862|  0.88|
|30|38.678536. . . . . . . . . . . . . . . . . . |   38,565,551|  42851|  1.19|
|35|37.814765. . . . . . . . . . . . . . . . . . |   17,395,043|  19328|  0.54|
|40|33.512648. . . . . . . . . . . . . . . . . . |  132,926,997| 147698|  4.10|
|41|-nan . . . . . . . . . . . . . . . . . . . . |      112,499|    125|  0.00|
|45|35.253821. . . . . . . . . . . . . . . . . . |    4,467,560|   4964|  0.14|
|50|39.191611. . . . . . . . . . . . . . . . . . |   54,872,503|  60970|  1.69|
|-----------------------------------------------------------------------------|
|TOTAL                                           |3,242,975,751|3603339|100.00|
+-----------------------------------------------------------------------------+


Compute without the -r flag and report:

```
r.stats.zonal method=average base=landcover cover=temperature output=average_temperature_per_landcover
r.report average_temperature_per_landcover -nh units=me,c,p --q output=zonal_statistics_report
```

Sample output where the categories are temperature strata:

+-----------------------------------------------------------------------------+
|               Category Information                 |    square|  cell|   %  |
|                  #|description                     |    meters| count| cover|
|-----------------------------------------------------------------------------|
|17.611052-17.619049|from  to . . . . . . . . . . . .| 1,166,400|  1296|  1.25|
| 18.226892-18.23489|from  to . . . . . . . . . . . .|   695,700|   773|  0.75|
|18.874724-18.882722|from  to . . . . . . . . . . . .| 2,393,100|  2659|  2.57|
|18.898717-18.906715|from  to . . . . . . . . . . . .|   551,700|   613|  0.59|
|18.978697-18.986695|from  to . . . . . . . . . . . .|   451,800|   502|  0.49|
| 19.034682-19.04268|from  to . . . . . . . . . . . .| 8,695,800|  9662|  9.34|
| 19.04268-19.050678|from  to . . . . . . . . . . . .|32,199,300| 35777| 34.60|
|19.138655-19.146653|from  to . . . . . . . . . . . .| 5,391,000|  5990|  5.79|
|19.210636-19.218634|from  to . . . . . . . . . . . .| 8,170,200|  9078|  8.78|
|19.570543-19.578541|from  to . . . . . . . . . . . .|17,059,500| 18955| 18.33|
|19.634526-19.642524|from  to . . . . . . . . . . . .| 1,031,400|  1146|  1.11|
|19.642524-19.650522|from  to . . . . . . . . . . . .|15,250,500| 16945| 16.39|
|-----------------------------------------------------------------------------|
|TOTAL                                               |93,056,400|103396|100.00|
+-----------------------------------------------------------------------------+


Building up on this, loop over a series of descriptive statistics:

for METHOD in min max range average variance ;do r.stats.zonal method=$METHOD base=landcover cover=temperature output=average_temperature_per_landcover -r ;done


A somewhat elaborated example to loop over cities and corresponding
landcover map tiles (filenames including the .tif extension):

--%<---
for CITY in $(v.db.select -c cities column=name,tile |sort -u) ;do

    # positional parameters
    set -- $(echo $CITY |tr "|" " ")

    # Take cate about Region & MASK
    g.region vector=cities -pa --q &&
    r.mask --o vector=cities where="name='$(echo ${1})'"
    g.region zoom=MASK &&

    # loop over Aggergated maps
    for AGGREGATE in $(t.rast.list aggregates -u column=name) ;do

        # loop over Methods
        for METHOD in min max range average variance ;do

            echo "$1 $2 $AGGREGATE $METHOD:"
            r.stats.zonal method=$METHOD base=${2%.tif} cover=$AGGREGATE output=${1}_${2%.tif}_${AGGREGATE}_${METHOD}_categories -r --o

        done

    done
    r.mask -r --v

done
--->%--

Nikos

>
>
>On 2018-02-12 4:38 AM, Nikos Alexandris wrote:
>>* Nicolas Cadieux <nicolas.cadieux at archeotec.ca> [2018-02-11 
>>23:37:13 -0500]:
>>
>>>Hi,
>>>I have a raster (int16) with multiple class (ex 1 = grass, 2 = 
>>>trees...) and a vector polygon file.  I want to have the zonal 
>>>statistic (cum, min, max, sum...) for every class.  What is the 
>>>best way to get this done in QGis 2.18.16 on Windows 10?
>>>
>>>Thanks,
>>>Nicolas
>>
>>Nicolas,
>>
>>not a direct answer, yet if you have access to GRASS' tools, then 
>>you could look
>>into v.rast.stats.html [0] (or even r.stats.zonal [1] if you have a
>>raster version of your "polygons").
>>
>>If I understand correctly, you want to get statistics per polygon, not
>>per class. I.e., min, max, sum and the rest, for each class inside
>>each vector polygon.
>>
>>v.rast.stats will derive the statistics in question from a raster map
>>for the parts covered by the vector map. Note, internally, the 
>>vector map is
>>converted to a raster.
>>
>>If you have two rasters, then your base map is the "multiple classes 
>>map".
>>Your cover map will be the "rasterised vector polygons", which you 
>>consider as
>>the zones to compute the statistics in question within. Very useful,
>>here, is the `-r` flag: "Create reclass map with statistics as category
>>labels".
>>
>>Nikos
>>
>>[0] https://grass.osgeo.org/grass74/manuals/v.rast.stats.html
>>[1] https://grass.osgeo.org/grass74/manuals/r.stats.zonal.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 228 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/qgis-user/attachments/20180212/3cd6f8af/attachment.sig>


More information about the Qgis-user mailing list