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

Glynn Clements glynn at gclements.plus.com
Thu Aug 11 11:06:52 EDT 2005

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

> 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>

More information about the grass-user mailing list