[GRASS-user] zonal statistics/metrics - beyond simple statistics

[Moritz Lennert]

On 11/04/18 20:58, Bernardo Santos wrote:
> Hi Helmut and Stefan,
> 
> thanks for sharing that with me, it is good to know that there are 
> similar ongoing approaches to solve related issues.
> 
> Helmut,
> very nice addon. However, it seems to me that the part on calculating 
> statistics is still based on v.rast.stats, am I right?
> 
> Stefan, the question you pose is indeed very related. I took a better 
> look at my code and I noticed that what really takes time is not the 
> process of making a mask, but what comes before: the rasterization of 
> one of the polygons of the input vector.
> What I do is:
> 1) g.region using the vector layer, aligning it to the first of the 
> input rasters
> 2) v.to.rast selecting only one of the polygons at a time from the 
> vector layer
> 3) r.mask using the rasterized polygon

You can actually use r.mask directly with the vector map using the 
vector and cat/where parameters.

> 4) g.region zoom=rasterized_polygon
> 5) calculate any function for all raster maps (in principle using a 
> loop), using this small region and the mask, and attaching the results 
> to the attribute table
> 6) repeat 1-5 for all polygons in the input vector
> 
> However, I noticed now that what takes time is v.rast.stats, since it 
> depends on the whole region selected, which in my case is much bigger 
> than a single polygon (the vector is a set of ~6,000 small polygons).

After step 4, your region should be only the size of the 
rasterized_polygon. If this is not the case there is an issue somewhere. 
If you use r.mask directly with the vector map you should be able to set 
the region to the mask using g.region raster=MASK.

> Do you guys know if there is a way of using a single polygon from a 
> vector to define the region using something like g.region?
> [this would be really great!!]
> something like
> g.region vector=vector_of_interest cat=1

g.region does not have cat/where options, but I agree that this might be 
useful.

More generally, and in support of what Stefan has already said, it might 
be helpful to get a bigger picture of what exactly you are doing and 
what your data looks at exactly, as maybe looping is not your best option.

Also, keeping everything in raster and text file realm as long as 
possible could speed up the process as attribute data handling is often 
a bottle neck. You can see the code of the i.segment.stats addon [1] for 
an example of this, including parallelization of the zonal statistics of 
the different raster maps.

Moritz

[1] 
https://trac.osgeo.org/grass/browser/grass-addons/grass7/imagery/i.segment.stats/i.segment.stats.py