[GRASS-user] subsampling to a coarser resolution

[Glynn Clements]

jamesmcc wrote:

> Greetings. I'm new to GRASS and after reading the help over and searching
> this forum I've decided to post. 
> 
> I have a raster at 30m resolution in which is nested another raster at 1.5m
> resolution. I want to average the 4 points in the center of the higher res
> raster (in [row,col] enumerating from 1: [10,10], [10,11], [11,10] and
> [11,11]) and assign this to the nesting cell of the lower res raster. Of
> course I want to do this over all cells in the raster, which is say 10x10. 

You can do this with:

	r.resamp.interp method=bilinear

For each cell in the current region, it takes the cell's centre
coordinates and samples the surface defined by the input map at that
point. For method=bilinear, the sample for a given x,y coordinate is
calculated as:

	x0 = floor(x)
	x1 = x0 + 1
	u = x - x0

	y0 = floor(y)
	y1 = y0 + 1
	v = y - y0

	c00 = input[y0][x0]
	c01 = input[y0][x1]
	c10 = input[y1][x0]
	c11 = input[y1][x1]

	c0 = c00 * (1-u) + c01 * u
	c1 = c10 * (1-u) + c11 * u

	c = c0 * (1-v) + c1 * v

The input map is read at its native resolution.

method=bicubic is similar, but uses a 4x4 window and bicubic
interpolation.

> Seems easy, but 
> 1. I'm not confident in any of the resample routines going from their
> descriptions (probably overlooking something. One problem is that
> r.neighbors and r.mfilter arent clear about going between regions and
> require odd neighborhoods both of which discourage me from trying them. I've
> also done an average resampling to 3m res, but how do i discard all the
> unwanted cells and get to 30m res?)

r.neighbors and r.mfilter generate an output map with the same
resolution as the input. You could subsequently use

	r.resamp.interp method=nearest

to downsample the data, but this is inefficient as you'll be
calculating many values which will subsequently be discarded.

The various r.resamp.* modules perform downsampling, i.e. they read
the input map at its native resolution and generate output at the
current region resolution, and only calculate the cell values which
appear in the output map.

> 2. How in the world would I check this? Can I print out the cells in the
> input raster and check those against the corresponding cells in the output
> raster. 

r.what can be used to obtain individual cell values, while r.out.ascii
will output a raster map as ASCII data. Both of these read their input
at the current region resolution, so you should first use

	g.region rast=...

if you want the raw map data without any resampling.

-- 
Glynn Clements <glynn at gclements.plus.com>