[GRASS-user] subsampling to a coarser resolution

GRASS GRASS user
1

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.

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?)
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.

For now, I'm just going to do it in my "native" programming language, cause
it really is a trival problem when you're familiar with the software... I'd
really like to know how to in GRASS.

THanks,

JM
--
View this message in context: http://n2.nabble.com/subsampling-to-a-coarser-resolution-tp3563929p3563929.html
Sent from the Grass - Users mailing list archive at Nabble.com.

2

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

3

Hi,

2009/9/2 jamesmcc <jlmccreight@gmail.com>

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.

While you can do this with GRASS, I tend to do these sort of processes with python and GDAL/OGR. If you are interested, I have put some notes here:
< http://sites.google.com/site/spatialpython/>

Jose

4

Thanks Glynn,

I didnt realize that bilinear was what I was doing, though I'm familiar with
it. This should work, though I immediately become concerned by the first
line of the documentation

r.resamp.interp - Resamples raster map layers to a finer grid using
interpolation.

since I'm not going to a finer resolution. I noticed that you are the
author, so I thought I should say something! :slight_smile:

One point of clarification in the equations you wrote, are the xy coords in
the input raster and in the current region are normalized by the resolution
of the input raster?

Lastly, r.what looks nice but I have to kinda scratch my head. While I
appreciate being able to query the raster value at a specific geographic
location, I use GIS so that IT will keep track of the positions. I want to
get the values based on row,column of the raster not north,east... clearly,
if i could figure out how to get north,east from row,column I'd just pipe
this to r.what, but I havent figured that out either.

Since I still cant check, I havent tried to implement any of this yet...

thanks for you help!

Glynn Clements wrote:

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@gclements.plus.com>
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

--
View this message in context: http://n2.nabble.com/subsampling-to-a-coarser-resolution-tp3563929p3567193.html
Sent from the Grass - Users mailing list archive at Nabble.com.

5

On Wed, Sep 2, 2009 at 4:42 PM, jamesmcc<jlmccreight@gmail.com> wrote:

Thanks Glynn,

I didnt realize that bilinear was what I was doing, though I'm familiar with
it. This should work, though I immediately become concerned by the first
line of the documentation

r.resamp.interp - Resamples raster map layers to a finer grid using
interpolation.

since I'm not going to a finer resolution. I noticed that you are the
author, so I thought I should say something! :slight_smile:

see also:
http://grass.osgeo.org/grass64/manuals/html64_user/r.resamp.stats.html
r.resamp.stats - Resamples raster map layers to a coarser grid using aggregation

Markus

6

Hi Markus-

I've used r.resamp.stats alot, but in this case it dosent appear to be
flexible enough to just select the cells I want to do stats on. It's still a
lovely routine, however. It would be nice to be able to supply it a matrix
mask to make it more flexible, then it could do what I want. Then, it could
cook me lunch too. :slight_smile:

I've now tested the bilinear method in r.resamp.interp, exactly as suggested
by Glynn and it works perfectly going to coarser resolution. I imported my
interpolates/aggregates made in IDL and r.stats shows only a zero category
for the difference of the rasters!

Still wondering about r.what and row,col inputs...

thanks!

Markus Neteler wrote:

On Wed, Sep 2, 2009 at 4:42 PM, jamesmcc<jlmccreight@gmail.com> wrote:

Thanks Glynn,

I didnt realize that bilinear was what I was doing, though I'm familiar
with
it. This should work, though I immediately become concerned by the first
line of the documentation

r.resamp.interp - Resamples raster map layers to a finer grid using
interpolation.

since I'm not going to a finer resolution. I noticed that you are the
author, so I thought I should say something! :slight_smile:

see also:
http://grass.osgeo.org/grass64/manuals/html64_user/r.resamp.stats.html
r.resamp.stats - Resamples raster map layers to a coarser grid using
aggregation

Markus
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

--
View this message in context: http://n2.nabble.com/subsampling-to-a-coarser-resolution-tp3563929p3567972.html
Sent from the Grass - Users mailing list archive at Nabble.com.

7

jamesmcc wrote:

I didnt realize that bilinear was what I was doing, though I'm familiar with
it. This should work, though I immediately become concerned by the first
line of the documentation

r.resamp.interp - Resamples raster map layers to a finer grid using
interpolation.

since I'm not going to a finer resolution. I noticed that you are the
author, so I thought I should say something! :slight_smile:

It isn't limited to upsampling, although that's the main use.

One point of clarification in the equations you wrote, are the xy coords in
the input raster and in the current region are normalized by the resolution
of the input raster?

Yes; those equations assume that x and y are in input cells, e.g.
0.5,0.5 is the centre of the top-left cell.

Lastly, r.what looks nice but I have to kinda scratch my head. While I
appreciate being able to query the raster value at a specific geographic
location, I use GIS so that IT will keep track of the positions. I want to
get the values based on row,column of the raster not north,east... clearly,
if i could figure out how to get north,east from row,column I'd just pipe
this to r.what, but I havent figured that out either.

You can obtain the coordinates using the output from r.info, e.g.:

  inmap=...
  row=...
  col=...
  eval `r.info -gs $inmap` # sets north,south,east,west,nsres,ewres
  n=`echo "$north - ($row + 0.5) * $nsres" | bc`
  e=`echo "$west + ($col + 0.5) * $ewres" | bc`
  r.what $inmap east_north=$e,$n

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