[GRASSLIST:7859] Re: Generating rasters at different resolutions

[Glynn Clements]

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 at gclements.plus.com>