[GRASS5] raster inaccuracy

GRASS GIS GRASS GIS developer
1

Folks,

I use Grass in part to support finite difference ground water models. There
is a substantial inaccuracy in Grass' raster representation that I've known of
for some time. The problem surfaced again recently, caused errors in my
reported results and required me to repeat several hours of work at
substantial costs. I think this problem has been present from early betas of
5.0 through 5.3 (our current installation).

The png file at

http://members.spinn.net/~roger/mismatch.png

illustrates the problem. The image was produced by the Grass PNG driver. It
shows the top left corner of a model grid. The lines are a vector map that
outlines the individual model cells. The colored blocks are a raster map
developed from the vector map using v.to.rast. The raster cells are not
illustrated at the same locations as the cells of the vector map.

The cells in this map are a mile across and the resolution of the region is
set to one mile. The actual size of the mismatch varies according to the
location of the map origin, but even with the origin set exactly at the corner
of the grid (where normally it is set) the raster cells do not correspond to
the cells in the vector map. The size of the mismatch is reduced by
increasing the resolution of the region. The difference is negligible if the
row and column dimension of the region are reduced to a small fraction of the
cell width - 1% worked well in this case, but such detailed resolutions are
often not practical and shouldn't be necessary.

The problem is not limited to the quality of the illustration. My worst
problem arose when the raster file was queried with d.what.rast. The
row-column location and category value produced by d.what.rast are the
location and value of the raster cell as illustrated at a location, not the
actual row-column location and category value at the location queried.

I mapped a sites layer illustrating well locations over the raster map then
queried the raster map to identify the row and column location for several
wells. I entered those wells into the ground water model based on the results
of the query. The wells were sometimes located in the correct model cell, but
in other cases were assigned to one of the 8 surrounding model cells.

I can only ask. *Please* get this problem fixed. Grass is a powerful tool,
but this inaccuracy makes it very difficult for me to use Grass in my work.

Roger Miller
Lee Wilson and Associates

2

I use Grass in part to support finite difference ground water models. There
is a substantial inaccuracy in Grass' raster representation that I've known of
for some time. The problem surfaced again recently, caused errors in my
reported results and required me to repeat several hours of work at
substantial costs. I think this problem has been present from early betas of
5.0 through 5.3 (our current installation).

The png file at

http://members.spinn.net/~roger/mismatch.png

is your region something like this:

e=100025
w=200075
n=300025
s=400075
res=50

or more to the point is your input data like the above and your region like:

e=100000
w=200000
n=300000
s=400000
res=50

?

as a test, do

g.region res=25 # or res="some common denominator"

and see if things start to line up in the XDRIVER.
(the PNG driver is much less tested than XDRIVER, maybe other bugs..)

I've had the same trouble with importing model data; setting the region
resolution to half, or setting the bounds as the first example above (25-75)
fixed it.. but you have to be careful it doesn't inadvertently switch back,
also when you do overlays with maps which use a 0-50 bounds settings the
computer is doing what you thing it is doing...

maybe it?

Hamish

3

Hamish,

The model region is defined by the authors of the model. The model origin is
at a more-or-less arbitrary location and the boundaries are such that it
divides equally into 188X256 rows and columns, each 5280 feet wide.

I'm not sure what the actual boundaries of the region were when I produced the
illustration. That region was set by zooming in on the corner then panning
to the northwest to get a clear view of the corner. The actual size of the
mismatch varies with the location of the origin and with the resolution of
the region. I can make the raster representation accurate by going to very
fine resolutions. Otherwise, if there is any situation in which the raster
map is accurate then it would probably -- as you suggest -- depend on a
getting the right combination of origin, resolution and cell size.

The view in the X monitor is essentially identical to the illustration
provided by the png monitor.

This is not a problem with importing data. The vector representation of the
model grid was constructed as a Grass ascii vector file. Using d.where the
grid looks correctly located and dimensioned. The raster file in this case
is not a model file. The raster file was created directly off the vector
representation of the grid using v.to.rast. The raster file and the vector
file represent the same data, but the raster representation is inaccurate,
with the size of the inaccuracy varying with the region resolution and
boundaries.

Roger Miller

On Tuesday 09 November 2004 16:45, you wrote:

> I use Grass in part to support finite difference ground water models.
> There is a substantial inaccuracy in Grass' raster representation that
> I've known of for some time. The problem surfaced again recently,
> caused errors in my reported results and required me to repeat several
> hours of work at substantial costs. I think this problem has been
> present from early betas of 5.0 through 5.3 (our current installation).
>
> The png file at
>
> http://members.spinn.net/~roger/mismatch.png

is your region something like this:

e=100025
w=200075
n=300025
s=400075
res=50

or more to the point is your input data like the above and your region
like:

e=100000
w=200000
n=300000
s=400000
res=50

?

as a test, do

g.region res=25 # or res="some common denominator"

and see if things start to line up in the XDRIVER.
(the PNG driver is much less tested than XDRIVER, maybe other bugs..)

I've had the same trouble with importing model data; setting the region
resolution to half, or setting the bounds as the first example above
(25-75) fixed it.. but you have to be careful it doesn't inadvertently
switch back, also when you do overlays with maps which use a 0-50 bounds
settings the computer is doing what you thing it is doing...

maybe it?

Hamish

4

Roger Miller wrote:

I use Grass in part to support finite difference ground water models. There
is a substantial inaccuracy in Grass' raster representation that I've known of
for some time. The problem surfaced again recently, caused errors in my
reported results and required me to repeat several hours of work at
substantial costs. I think this problem has been present from early betas of
5.0 through 5.3 (our current installation).

The png file at

http://members.spinn.net/~roger/mismatch.png

illustrates the problem. The image was produced by the Grass PNG driver. It
shows the top left corner of a model grid. The lines are a vector map that
outlines the individual model cells. The colored blocks are a raster map
developed from the vector map using v.to.rast. The raster cells are not
illustrated at the same locations as the cells of the vector map.

The cells in this map are a mile across and the resolution of the region is
set to one mile. The actual size of the mismatch varies according to the
location of the map origin, but even with the origin set exactly at the corner
of the grid (where normally it is set) the raster cells do not correspond to
the cells in the vector map. The size of the mismatch is reduced by
increasing the resolution of the region. The difference is negligible if the
row and column dimension of the region are reduced to a small fraction of the
cell width - 1% worked well in this case, but such detailed resolutions are
often not practical and shouldn't be necessary.

The problem is not limited to the quality of the illustration. My worst
problem arose when the raster file was queried with d.what.rast. The
row-column location and category value produced by d.what.rast are the
location and value of the raster cell as illustrated at a location, not the
actual row-column location and category value at the location queried.

I mapped a sites layer illustrating well locations over the raster map then
queried the raster map to identify the row and column location for several
wells. I entered those wells into the ground water model based on the results
of the query. The wells were sometimes located in the correct model cell, but
in other cases were assigned to one of the 8 surrounding model cells.

I can only ask. *Please* get this problem fixed. Grass is a powerful tool,
but this inaccuracy makes it very difficult for me to use Grass in my work.

First, I suspect that any discrepancy between vectors and rasters is
more likely to be in the vector code.

Also, bear in mind that GRASS' raster coordinate system is
"grid-registered" rather than "node-registered". E.g. the north-east
edge of the region is the location of the north-east corner of the
north-east cell, not of the centre of the north-east cell. If you are
given the region coordinates using a node-registered convention, you
have to expand the region by half a cell in each direction for use
with GRASS.

[BTW, GDAL has been known to get this wrong on occastion, i.e. failing
to add half a cell or adding it when it shouldn't; if you are
importing any data with GDAL, ensure that you have a recent version.]

For similar reasons, display modules which draw lines tend to be
problematic, as the display architecture uses integer coordinates, so
you can only place vertices on the corners of cells, not on their
centres.

In any case, in order to look into this any further, I would need
specific cases, i.e. a specific sequence of commands which can be
performed on data which I either have (e.g. spearfish, global) or can
generate.

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

5

On Wed, 10 Nov 2004 14:49:24 +0000, glynn wrote

First, I suspect that any discrepancy between vectors and rasters is
more likely to be in the vector code.

The output from d.where, d.what.vect and d.what.rast consistently show that
where there is a difference the vector map is in about the right place and the
raster map is in the wrong place.

Also, bear in mind that GRASS' raster coordinate system is
"grid-registered" rather than "node-registered". E.g. the north-east
edge of the region is the location of the north-east corner of the
north-east cell, not of the centre of the north-east cell. If you are
given the region coordinates using a node-registered convention, you
have to expand the region by half a cell in each direction for use
with GRASS.

This difference doesn't effect my current problem, but I'm aware that it does
need to be kept in mind.

For similar reasons, display modules which draw lines tend to be
problematic, as the display architecture uses integer coordinates, so
you can only place vertices on the corners of cells, not on their
centres.

The cells in this case are thousands of units across, so the difference
between integer and floating point coordinates should be negligible.

In any case, in order to look into this any further, I would need
specific cases, i.e. a specific sequence of commands which can be
performed on data which I either have (e.g. spearfish, global) or can
generate.

I put together several examples that can be had from

http://members.spinn.net/~roger/example.tgz

The tar file produces a GRASS XY location named "example." The examples are
in the PERMANENT map set.

Here is the step-by-step procedure:

copy the tar file to our Grass data file hierarchy and unarchive the file

start Grass, using location "example" and mapset "PERMANENT"

"g.list vect" should list 6 vector files. Each file contains a simple grid
with square cells each 5280x5280.

Start with 2x2example. The procedure for all of the examples is identical.

"v.info 2x2example" to get the limits of the grid.

Use "g.region" to set the region to match the grid. In the case of 2x2example
the limits are n=10560 s=0 w=0 e=10560 nsres=5280 ewres=5280. In all examples
the resolution is 5280. The region so defined should have two rows and two
columns.

Ideally it seems like you should be able to set the region with "g.region
vect=2x2example." That may work, but in most examples it does not. If the
region is not set exactly to the dimensions of the grid then the next step
will produce unwanted results.

"v.to.rast i=2x2example o=2x2example"
"d.mon start=x0"
"d.rast -o 2x2example"
"d.vect 2x2example"

This should produce an accurate map, with the lines of the vector grid exactly
outlining the colored cells in the raster map.

"d.pan" and click on the image to recenter (e.g., click your pointer in the
center of the northwest cell). Repeat this step several times with different
center locations. With some centers the map is accurate. With other centers
the map is drawn with the raster file offset by half a cell width in both the
x and y directions. When there is an offset you can use d.where, d.what.rast
or d.what.vect to confirm that the vector map is correctly located (southwest
corner at x=0, y=0) and the raster map is offset.

As an aside, d.what.vect in my 5.3 installation does not work correctly. It
identifies the two cells on the top row correctly and reports consistent
measurements. It misidentifies the two cells on the lower row and reports -0
for all dimensions.

The additional examples illustrate different effects.

"3x3example" uses a 3x3 grid. Grids with odd-numbered dimensions don't have
the same problem with inaccuracy.

"10x10example and 11x11example" are similar to 2x2example and 3x3example and
are just there to show that the same effect occurs on slightly large grids.

"14x11example" has 14 columns and 11 rows. In this case the grid and the
raster map always agree in the y direction (odd number of rows) but often
disagree in the x direction (even number of columns)

"shiftexample" has 14 columns and 11 rows, but the origin is shifted by 880 in
both x and y. The shift does not alter the results.

If the region is not set correctly when the raster file is created (for
instance, if "g.region vect=..." is used to set the region) then the raster
map will have additional problems and will not correspond accurately to the
vector grid. That seems to be a different problem.

Roger Miller

6

Roger Miller wrote:

I put together several examples that can be had from

http://members.spinn.net/~roger/example.tgz

The tar file produces a GRASS XY location named "example." The examples are
in the PERMANENT map set.

Here is the step-by-step procedure:

copy the tar file to our Grass data file hierarchy and unarchive the file

start Grass, using location "example" and mapset "PERMANENT"

"g.list vect" should list 6 vector files. Each file contains a simple grid
with square cells each 5280x5280.

Start with 2x2example. The procedure for all of the examples is identical.

"v.info 2x2example" to get the limits of the grid.

Use "g.region" to set the region to match the grid. In the case of 2x2example
the limits are n=10560 s=0 w=0 e=10560 nsres=5280 ewres=5280. In all examples
the resolution is 5280. The region so defined should have two rows and two
columns.

Ideally it seems like you should be able to set the region with "g.region
vect=2x2example." That may work, but in most examples it does not. If the
region is not set exactly to the dimensions of the grid then the next step
will produce unwanted results.

"v.to.rast i=2x2example o=2x2example"
"d.mon start=x0"
"d.rast -o 2x2example"
"d.vect 2x2example"

This should produce an accurate map, with the lines of the vector grid exactly
outlining the colored cells in the raster map.

"d.pan" and click on the image to recenter (e.g., click your pointer in the
center of the northwest cell). Repeat this step several times with different
center locations. With some centers the map is accurate. With other centers
the map is drawn with the raster file offset by half a cell width in both the
x and y directions. When there is an offset you can use d.where, d.what.rast
or d.what.vect to confirm that the vector map is correctly located (southwest
corner at x=0, y=0) and the raster map is offset.

OK; I understand you now. At first glance, this appears to be a
problem with the display library. I'll look into it, and get back to
you.

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

7

Glynn Clements wrote:

> I put together several examples that can be had from
>
> http://members.spinn.net/~roger/example.tgz
>
> The tar file produces a GRASS XY location named "example." The examples are
> in the PERMANENT map set.
>
> Here is the step-by-step procedure:
>
> copy the tar file to our Grass data file hierarchy and unarchive the file
>
> start Grass, using location "example" and mapset "PERMANENT"
>
> "g.list vect" should list 6 vector files. Each file contains a simple grid
> with square cells each 5280x5280.
>
> Start with 2x2example. The procedure for all of the examples is identical.
>
> "v.info 2x2example" to get the limits of the grid.
>
> Use "g.region" to set the region to match the grid. In the case of 2x2example
> the limits are n=10560 s=0 w=0 e=10560 nsres=5280 ewres=5280. In all examples
> the resolution is 5280. The region so defined should have two rows and two
> columns.
>
> Ideally it seems like you should be able to set the region with "g.region
> vect=2x2example." That may work, but in most examples it does not. If the
> region is not set exactly to the dimensions of the grid then the next step
> will produce unwanted results.
>
> "v.to.rast i=2x2example o=2x2example"
> "d.mon start=x0"
> "d.rast -o 2x2example"
> "d.vect 2x2example"
>
> This should produce an accurate map, with the lines of the vector grid exactly
> outlining the colored cells in the raster map.
>
> "d.pan" and click on the image to recenter (e.g., click your pointer in the
> center of the northwest cell). Repeat this step several times with different
> center locations. With some centers the map is accurate. With other centers
> the map is drawn with the raster file offset by half a cell width in both the
> x and y directions. When there is an offset you can use d.where, d.what.rast
> or d.what.vect to confirm that the vector map is correctly located (southwest
> corner at x=0, y=0) and the raster map is offset.

OK; I understand you now. At first glance, this appears to be a
problem with the display library. I'll look into it, and get back to
you.

I can see why this might be confusing; it had me confused for a while.

But this is the correct behaviour (from d.rast; d.pan is another
issue) It may not necessarily the behaviour which you want, but d.rast
is doing what it is supposed to.

The key point is that, in common with almost every other program which
operates upon rasters, d.rast isn't reading (or displaying) the
original raster, but the raster after it has been resampled according
to the current region settings.

If the current region has the same resolution as the raster, but isn't
aligned to the same grid, the result is that the raster which is
displayed on the monitor will be shifted by up to +/- half a cell.

This isn't an issue for vectors, as they aren't aligned to a grid.

If there's a flaw here, it's that d.pan appears to be aligning the
region to a half-cell grid. I have no idea why it does this, but it
probably shouldn't be performing any alignment. IOW, after panning,
the vector and the raster should be misaligned by an arbitrary amount
(in the range +/- half a cell), rather than by either half a cell or
nothing.

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

8

If there's a flaw here, it's that d.pan appears to be aligning the
region to a half-cell grid. I have no idea why it does this, but it
probably shouldn't be performing any alignment. IOW, after panning,
the vector and the raster should be misaligned by an arbitrary amount
(in the range +/- half a cell), rather than by either half a cell or
nothing.

'd.pan' is superseded by 'd.zoom -p' in 5.7, so I wouldn't bother with
it too much. But does 'd.zoom -p' (and d.zoom I guess) do the right thing?

Hamish

9

Glynn, thanks for looking into this for me.

On Thu, 11 Nov 2004 18:15:28 +0000, Glynn Clements wrote

The key point is that, in common with almost every other program
which operates upon rasters, d.rast isn't reading (or displaying)
the original raster, but the raster after it has been resampled according
to the current region settings.

If the current region has the same resolution as the raster, but
isn't aligned to the same grid, the result is that the raster which
is displayed on the monitor will be shifted by up to +/- half a cell.

I see what you mean. It leaves me with a problem. In the general case grass
raster displays and grass raster queries are not accurate; they may be grossly
inaccurate. That is because the original data when displayed are resampled to
a grid that may have little in common with the original data array. In my
examples the data were displayed in a grid that overlapped as little as 25% of
the actual data cell. The d.what.rast query returns the data *as displayed*
not the original data, so it propogates the problem.

It appears to me that there are two (and only two) cases in which the raster
display and the raster queries will always be accurate:

1) The region resolution is set to the pixel size for the current display.
For instance, in my 2x2 example with the default 640x480 monitor size that
means setting the row and column sizes to 22 units.

2) The region boundaries are set to correspond exactly to boundaries between
rows and columns in the original data and the region resolution is set so
there is an integral number of region rows and columns for each row and column
in the original data.

I don't understand what gets displayed (hence queried) in the case where the
data cells are smaller than the pixel size. That seems to be a largely
different problem, but another problem with accuracy.

The first solution is easy. The simple solution is generally useful for the
purpose of display and query but it would not work for most analyses of raster
data -- through r.mapcalc, for instance.

The second case is more complicated, but I think it would work for all cases
where the pixel size is much smaller than the resolution of the data.

It is at least a burden on users to require them to make sure those conditions
are met, and probably really too much to expect from most users. Certainly I
don't like the prospect of manually checking and adjusting the region every
time I need to query a raster map.

Is there some way to impose those limits on the region definitions in an
automatic fashion?

This isn't an issue for vectors, as they aren't aligned to a grid.

If there's a flaw here, it's that d.pan appears to be aligning the
region to a half-cell grid. I have no idea why it does this, but it
probably shouldn't be performing any alignment. IOW, after panning,
the vector and the raster should be misaligned by an arbitrary amount
(in the range +/- half a cell), rather than by either half a cell or
nothing.

d.pan's behavior might be a problem, but it is a different problem. It did
have the advantage of maximizing the error in my examples. The same (but less
extreme) problem arises with the region is manually adjusted by changing the
region boundaries and/or resolution.

Roger Miller

10

Roger Miller wrote:

> The key point is that, in common with almost every other program
> which operates upon rasters, d.rast isn't reading (or displaying)
> the original raster, but the raster after it has been resampled according
> to the current region settings.
>
> If the current region has the same resolution as the raster, but
> isn't aligned to the same grid, the result is that the raster which
> is displayed on the monitor will be shifted by up to +/- half a cell.

I see what you mean. It leaves me with a problem. In the general case grass
raster displays and grass raster queries are not accurate; they may be grossly
inaccurate. That is because the original data when displayed are resampled to
a grid that may have little in common with the original data array. In my
examples the data were displayed in a grid that overlapped as little as 25% of
the actual data cell. The d.what.rast query returns the data *as displayed*
not the original data, so it propogates the problem.

However, this isn't just an issue for d.rast. It applies to (more or
less) r.<everything>.

One way to look at it is: "What You See Is What GRASS Gets", i.e. what
is displayed on the monitor is what a raster (r.*) program will get if
it reads the raster while the current region settings are in effect
(if it doesn't, that *is* a bug).

It appears to me that there are two (and only two) cases in which the raster
display and the raster queries will always be accurate:

1) The region resolution is set to the pixel size for the current display.
For instance, in my 2x2 example with the default 640x480 monitor size that
means setting the row and column sizes to 22 units.

2) The region boundaries are set to correspond exactly to boundaries between
rows and columns in the original data and the region resolution is set so
there is an integral number of region rows and columns for each row and column
in the original data.

Correct.

I don't understand what gets displayed (hence queried) in the case where the
data cells are smaller than the pixel size. That seems to be a largely
different problem, but another problem with accuracy.

Actually, it's exactly the same problem. The sampling algorthim is
identical regardless of whether region's resolution is higher or lower
than that of the original raster.

Possibly the easiest way to look at it is this: a raster map defines a
function over the plane (the set R^2), i.e. it has infinite resolution
and infinite extent (although any non-null values are contained in a
finite subregion).

The value at any point on the plane is the value from the raster cell
which contains the point, i.e. the cell whose centre is closest to
that point. If the point lies outside of the raster, the value is
null.

However, a program can't read an infinite amount of data, so it
obtains a finite number of samples, forming a rectangular grid. The
program can choose the dimensions and spacing of the grid, although
most allow the user to choose this via the current region (i.e. the
WIND file).

If you partition the rectangle defined by the N/S/E/W values into
rows*cols "tiles", the sample points are at the centre of each tile,
i.e. at the coordinates:

  x = west + (east - west ) * (col + 1/2) / cols
  y = north - (north - south) * (row + 1/2) / rows

for col in [0..cols-1] and row in [0..rows-1] inclusive.

If the region's resolution (spacing) is finer (smaller distance) than
that of the raster, several adjacent samples may fall into the same
cell in the original raster. OTOH, if the region's resolution is
coarser, some cells in the original raster will be omitted from the
sample. If the resolutions are identical, each cell should be sampled
exactly once.

The first solution is easy. The simple solution is generally useful for the
purpose of display and query but it would not work for most analyses of raster
data -- through r.mapcalc, for instance.

The second case is more complicated, but I think it would work for all cases
where the pixel size is much smaller than the resolution of the data.

It is at least a burden on users to require them to make sure those conditions
are met, and probably really too much to expect from most users. Certainly I
don't like the prospect of manually checking and adjusting the region every
time I need to query a raster map.

Is there some way to impose those limits on the region definitions in an
automatic fashion?

"g.region rast=..." will trivially force the second case (the
"integral number" being 1). Although it might be useful to have a
scale= option for g.region to allow rescaling by an integral amount
("g.region res=..." suffices when the resolution is an integer, but
not in general).

Actually, it might also be useful to have rows= and cols= options, so
that you can set those directly rather than setting res=.

The first case is somewhat more complex, as the monitor's aspect ratio
won't generally match that of the raster, so you first have to
determine whether to "pad" the region horizontally or vertically, then
determine the amount of padding.

A third option would be to provide an option to d.rast to
automatically use the original raster's resolution, so that it is only
resampled once (in general, it gets resampled from the original grid
to the current region, and then again to the monitor's pixel grid).

However, this is risky, because what you're seeing isn't what r.* will
get. It's probably better to force the user to explicitly set the
region so that the display matches the "view" that r.* will receive.

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

11

On Thursday 11 November 2004 20:05, Glynn Clements wrote:

"g.region rast=..." will trivially force the second case (the
"integral number" being 1).

That should be true, anyway. When I was putting together those examples I
came across some cases (possibly using a vector file to set the region) where
the bounds of the region did not match the bounds of the file.

Although it might be useful to have a
scale= option for g.region to allow rescaling by an integral amount
("g.region res=..." suffices when the resolution is an integer, but
not in general).

There are a number of possible solutions. It seems like the biggest problem
would arise if the display contained more than one raster file with different
boundaries and resolutions.

Actually, it might also be useful to have rows= and cols= options, so
that you can set those directly rather than setting res=.

I think we have ewres= and nsres= for that purpose.

The first case is somewhat more complex, as the monitor's aspect ratio
won't generally match that of the raster, so you first have to
determine whether to "pad" the region horizontally or vertically, then
determine the amount of padding.

I don't think that would be all that difficult. Your d.info code can be used
to produce the monitor dimensions (nx, ny) and the map dimensions are X, Y.
The correct pixel size (region row or column width) should be the larger of
X/nx or Y/ny.

A third option would be to provide an option to d.rast to
automatically use the original raster's resolution, so that it is only
resampled once (in general, it gets resampled from the original grid
to the current region, and then again to the monitor's pixel grid).

However, this is risky, because what you're seeing isn't what r.* will
get. It's probably better to force the user to explicitly set the
region so that the display matches the "view" that r.* will receive.

As conditions are now the problem can only be solved with user education,
which is probably the hardest solution of all. It is probably better in the
long run to build some kind of solution in the code.

For my part I do a lot of my work these days through a custom front end that I
wrote mostly to simplify some of the display and query procedures. I think
that for now at least I will try setting it up to keep the region row and
column dimensions equal to the pixel size. The region will have to be
manually reset before using r.mapcalc, r.stats, and so on, but that is
normally true even without making any further changes.

Roger Miller