Jose Gomez-Dans wrote:
I have generated rasters from vectors at a very fine resolution
(related to the accuracy of the vector data), and I want to aggregate
these cells to generate rasters at a much coarser resolution. For each
cell in the aggregated raster, the value should be the value that is
majoritary on the high resolution raster, provided it also exceeds an
area threshold. Additionally, it would be nice to have a second raster
with the percentage of area of the chosen category.
There isn't a simple way to do this.
The two approaches which spring to mind are:
1. Use r.statistics with a base map where each cell has a unique
category, e.g.:
g.region res=$low_res
cols=`g.region -p | fgrep cols: | cut -c 6-`
r.mapcalc "base = row() * $cols + col()"
g.region res=$high_res
r.statistics base=base cover=inmap output=mode method=mode
r.mapcalc 'ismode = (inmap == mode)'
r.statistics base=base cover=ismode method=average output=area
r.mapcalc "outmap = if(area > $min_area,mode,null())"
g.remove rast=base,mode,ismode,area
2. Use r.neighbors and downsample the output, e.g.:
g.region res=$high_res
r.neighbors input=inmap output=mode.hi method=mode size=$size
g.region res=$low_res
r.resample input=mode.hi output=mode
g.region res=$high_res
r.mapcalc 'ismode = (inmap == mode)'
r.neighbors inmap=ismode output=area.hi method=average size=$size
g.region res=$low_res
r.resample input=area.hi output=area
g.region res=$high_res
r.mapcalc "outmap = if(area > $min_area,mode,null())"
g.remove rast=mode.hi,mode,ismode,area.hi,area
Note that r.neighbors will only work with odd-sized an neighbourhood.
Enough people ask for this that it would be useful to have something
which combines r.resample and r.neighbors using gridded
neighbourhoods.
In principle, I had thought about using r.mapcalc, but I haven't
sussed out how to select the neighborhood to be considered. There's
the ewres and the nsres operators, but I don't think I understand how
to use them, or how to carry out what I want to do!
r.mapcalc doesn't do neighbourhood operations. You can implement them
manually using multiple offset terms, but this quickly becomes
impractical unless the neighbourhood size is very small.
--
Glynn Clements <glynn@gclements.plus.com>