Dears,
some background: I am exploring GRASS GIS' capabilities to run an
algorithm I built, by exclusively using raster maps and `r.*` modules.
Alternatives (i.e. Python based solutions) excluded.
Algorithm
1 B ← {0, .., m-1}
2 T ← {0, .., n-1}
3 W ← 0
4 R ← 0
5 F ← 0
6 MASK ← HQR
7 foreach {b} ⊆ B do
8 X ← 0
9 foreach {t} ⊆ T do
10 Wt ← Et * Ct
11 W ← W⋃{Wt}
12 S ← ∑t∈Wt
13 foreach t ← T do
14 Rt ← Wt / S
15 R ← R⋃{Rt}
16 X ← X⋃{R}
where:
B is a set of (aggregation) boundaries (raster categories)
T is a set of land types (raster categories)
E denotes extent of "elements" in T (number of cells * area of current cell)
C denotes a corresponding weighting factor for elements in T (labels for raster categories)
W is an (empty) set to hold weighted extents of "elements" in T
S denotes the sum of all elements in W
R is an (empty) set to hold ratios of elements in W to S
X is an (empty) set to hold all sets R for each iteration over the elements of B
Running it in GRASS GIS
Lines 10-11 compute a weighted extent for each "element" in T.
To exemplify, given
- a CELL map 'types' with zones of interest, i.e.
r.category types
1 Urban
2 Cropland
3 Woodland and forest
4 Grassland
5 Heathland and shrub
6 Sparsely vegetated land
7 Wetland
- a CELL map depicting areas of interest 'cells_of_interest', i.e.
r.category cells_of_interest
9
- a CELL 'scores' map with *float* values set as labels, i.e.
r.category scores
1 0.1
2 0.2
3 0.3
4 0.4
5 0.5
6 0.6
7 0.7
one way to do this in GRASS is using `r.mapcalc`'s '@' notation:
# number of cells
r.stats.zonal base=types cover=cells_of_interest method=count output=cells -r
# extent
r.mapcalc "extent = @cells * area()"
r.stats.zonal base=types cover=extent method=average output=extent -r --o
# weighted extents
r.mapcalc "weighted = @extent * @scores" --o
r.stats.zonal base=types cover=weighted method=average output=weighted -r --o
The 'weighted' map features now:
r.category weighted
1 15668750.000000
2 54460500.000000
3 78163500.000000
4 15247000.000000
5 -nan
6 15000.000000
7 4294500.000000
Note, also the extension `r.area` can be used somewhere.
Questions
Is it possible to compute the sum of values that are
set as labels to a raster map (here, the 'weighted' map),
using raster modules? [Line 12 of the algorithm]
Essentially, how to compute and re-use the
sum_of = 15000 + 4294500 + 15247000 + 15668750 + 54460500 + 78163500
in an `r.mapcalc` expression?
If this is not directly possible, could the 'sum_of' be, for example, a label
to a raster map's category?
Perhaps, the same value 'sum_of' set as label to all of the categories
involved above?
If the above is doable, then, there should be a way to compute the ratio for
each element in W to S (Line 14). Something like:
r.mapcalc "ratios = @weighted / sum_of"
Thanks, Nikos