[GRASS5] [bug #870] (grass) r.proj problem: boxes appearing

this bug's URL: http://intevation.de/rt/webrt?serial_num=870
-------------------------------------------------------------------------

Hi Morten
(dear RT bugtracker)

I am sending this to the bug tracker as the information shall
be recorded.

To repeat the bug description: the maps generated by r.proj
- might be shifted by 0.5 to 1 pixel (not sure!)
- the resulting maps contain a sort of "box structure" where
   the information differs from the neighbors. This only applies to
   the border of the boxes, not the box contents. In a Gauss-Boaga
   (transverse mercator) projection the boxes are 500m x 500m
   everywhere like a chess-board.
   Note that this is only visible when generating an aspect map
   in case you re-project an elevation map.

Now more discussion:

On Wed, Nov 28, 2001 at 02:40:32PM +0100, Morten Hulden wrote:

On Wed, 28 Nov 2001, Markus Neteler wrote:

> > Are you running r.proj on a _shaded_ map? I don't think any of the
> > interpolation methods will give good results on anything else than a pure
> > elevation map.
>
> no - of course I used the elevation map. The aspect I only use to verify
> the result as it is often quite useful to identify artefacts.

I should explain what I did:

- r.proj'ed the elevation map
- run r.slope.aspect on the result
- displayed it, finding the artefacts
- sending *this* snapshot to you. The color table of
   the DEM doesn't show much, but d.what.rast seems to
   confirm that boxes are there.

> Sorry: I didn't want to blame you for anything! I was pretty sure that
> I am dealing with a long-term problem introduced by the original
> author. However, you are currently the only person knowing much
> about this module :slight_smile:

Not what i meant. No offence taken. I just wanted to make sure you are
aware that the interpolation methods subroutines in r.proj are separate
from the initial trimming and adjusting in the main routine, and in
bordwalk.c.

Ah - o.k.

And I have really never looked closely into what is happening
inside nearest.c, bipolar.c and cubic,c. Could be some obscure bug there,
but I will have to get my Grass burning again before I can do any further
programming. Right now not possible. Next year when I'm out of work
another story.

> > But I found out (by trial) that the only usable interpolation method is
> > the 'nearest neighbor' (default). The two other methods (bipolar and
> > cubic) tend to change ("interpolate") the coastlines and shape of lakes
> > too much.
>
> Shall I disable them to prevent user confusion?

When I pointed out that 'nearest' is the only acceptable method to one
person earlier, he replied that man pages say that bipolar and cubic
should give better results. Indeed, they take into account more
surrounding cells, so in _theory_ they should be better at interpolating.
But in practice we don't want interpolation of coast lines and of the
shape of lakes - we want familiar looking coasts and lakes. In general,
interpolation around a natural flat surface (sea level, lakes) gives
undesirable results.

Maybe disabling the other methods is too drastic. After all they may have
their uses. But people should be aware, and perhaps man page should be
fixed not to mislead them. I can't argue against theory, what I say is
based on what I leared from testing.

I suggest that I ask someone with a GRASS running and keep you posting.
Do you follow the "grass5" list? Otherwise I can cc to you.

> Perhaps I try again with a synthetic map, a tilted plane and a
> sine curve composition.

I have never seen artifacts like the ones you show, but I'm very
interested to hear you results.
In plain language what r.proj does after all the initial trimming and
adjusting is this

  for each cell in the output location:
  -projects the center of the cell into the input location
  -finds the value of the cell in the input map into which the
  projected point fell (if you use the default 'nearest neighbor' method)
  -sets the value of the cell in the output map to that value

So what can go wrong?

Great difference in resolution between input and output map may result in
some more or less random cell values. But setting the current region to
match the input map's resolution before projecting should fix that. And
you say you did this already, so I just don't know.

Another thing, similar in effect but different reason, is that because of
the narrowing and divergence of meridians in some projections, we may hit
cases where the center of many different cells in the output location
projects into the same cell of the input map. The result is that many
adjacent cells in the output map gets the same value.

So these two cases, resolution difference and projection difference is
what I come up with. Neiter one is really a bug, but rather a result of
the way r.proj works, and restrictions that are inherent to the theory.

Please let me know what your synthetic map tests say.

best regards
Morten

We need more volunteers to test/identify the problem.
What to do: re-project an elevation map between two LOCATIONS in
different projections. Generate the aspect map with r.slop.aspect,
inspect it for any strange systematic error.

Thanks, Morten, for your comments!

Markus

--- Headers Follow ---

From neteler@itc.it Wed Dec 12 16:22:06 2001

Return-Path: <neteler@itc.it>
Delivered-To: grass-bugs@lists.intevation.de
Received: from mail.intevation.de (aktaia [212.95.126.10])
  by lists.intevation.de (Postfix) with ESMTP id 18A7613A02
  for <grass-bugs@lists.intevation.de>; Wed, 12 Dec 2001 16:22:06 +0100 (CET)
Received: from camelot.itc.it (camelot.itc.it [195.223.171.5])
  by mail.intevation.de (Postfix) with ESMTP id A599B1B77C
  for <grass-bugs@intevation.de>; Wed, 12 Dec 2001 16:22:05 +0100 (CET)
Received: from artemide.itc.it (artemide [10.0.10.10])
  by camelot.itc.it (8.11.3/8.11.3) with ESMTP id fBCFLwd23895;
  Wed, 12 Dec 2001 16:21:58 +0100 (MET)
Received: from thuille.itc.it. (thuille [10.40.0.110])
  by artemide.itc.it (8.11.3/8.11.3) with ESMTP id fBCFLs712189;
  Wed, 12 Dec 2001 16:21:55 +0100 (MET)
Received: (from neteler@localhost)
  by thuille.itc.it. (8.11.6/8.11.2) id fBCFLqx04301;
  Wed, 12 Dec 2001 16:21:52 +0100
Date: Wed, 12 Dec 2001 16:21:52 +0100
From: Markus Neteler <neteler@itc.it>
To: Morten Hulden <morten@ngb.se>
Cc: grass-bugs@intevation.de
Subject: r.proj problem: boxes appearing
Message-ID: <20011212162152.A4293@itc.it>
References: <20011128114509.A23640@itc.it> <Pine.LNX.4.33.0111281350200.1197-100000@tor.ngb.se>
Mime-Version: 1.0
Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline
User-Agent: Mutt/1.2.5i
In-Reply-To: <Pine.LNX.4.33.0111281350200.1197-100000@tor.ngb.se>; from morten@ngb.se on Wed, Nov 28, 2001 at 02:40:32PM +0100
X-Spam-Status: No, hits=0 required=5 tests=

-------------------------------------------- Managed by Request Tracker

On Wed, 12 Dec 2001 16:22:07 +0100 (CET), Request Tracker <grass-bugs@intevation.de> wrote:

this bug's URL: http://intevation.de/rt/webrt?serial_num=870
-------------------------------------------------------------------------

Hi Morten
(dear RT bugtracker)

I am sending this to the bug tracker as the information shall
be recorded.

To repeat the bug description: the maps generated by r.proj
- might be shifted by 0.5 to 1 pixel (not sure!)

   Projection introduces error.

- the resulting maps contain a sort of "box structure" where
   the information differs from the neighbors. This only applies to
   the border of the boxes, not the box contents. In a Gauss-Boaga
   (transverse mercator) projection the boxes are 500m x 500m
   everywhere like a chess-board.
   Note that this is only visible when generating an aspect map
   in case you re-project an elevation map.

   Are you certain the original data doesn't contain these artifacts?
   Perhaps the projection exacerbated the effect?

> And I have really never looked closely into what is happening
> inside nearest.c, bipolar.c and cubic,c. Could be some obscure bug there,
> but I will have to get my Grass burning again before I can do any further
> programming. Right now not possible. Next year when I'm out of work
> another story.
>
> > > But I found out (by trial) that the only usable interpolation method is
> > > the 'nearest neighbor' (default). The two other methods (bipolar and
> > > cubic) tend to change ("interpolate") the coastlines and shape of lakes
> > > too much.
> >
> > Shall I disable them to prevent user confusion?

No. Do not disable. Bilinear/Cubic are appropriate for imagery.

--
Eric G. Miller <egm2@jps.net>

On Wed, Dec 12, 2001 at 08:21:29AM -0800, Eric G. Miller wrote:

On Wed, 12 Dec 2001 16:22:07 +0100 (CET), Request Tracker <grass-bugs@intevation.de> wrote:

> this bug's URL: http://intevation.de/rt/webrt?serial_num=870
> -------------------------------------------------------------------------
>
> Hi Morten
> (dear RT bugtracker)
>
> I am sending this to the bug tracker as the information shall
> be recorded.
>
> To repeat the bug description: the maps generated by r.proj
> - might be shifted by 0.5 to 1 pixel (not sure!)

   Projection introduces error.

yes, I know. But the map looks shifted which is too much error.
Again, perhaps someone can try as well.
I have to explain that I am working in two adjacent zones.
I have data for the same geographic location which are projected
to these two zones. My test was to re-project the elevation model
from the second zone to the first zone, then to overlay the
maps to see if any difference is present. For that I diff'ed with
r.mapcalc (showing expected projection errors appearing as distributed
numerical differences). Additionally I generated the aspect
maps which *seem* to show a sort of shift.

To test r.proj, you need data for the same geographic location which are
projected into two different zones or projections.

> - the resulting maps contain a sort of "box structure" where
> the information differs from the neighbors. This only applies to
> the border of the boxes, not the box contents. In a Gauss-Boaga
> (transverse mercator) projection the boxes are 500m x 500m
> everywhere like a chess-board.
> Note that this is only visible when generating an aspect map
> in case you re-project an elevation map.

   Are you certain the original data doesn't contain these artifacts?
   Perhaps the projection exacerbated the effect?

I am very sure that the original data don't contain the artifacts.
If you try with a high-res DEM, you may get them as well.

> > And I have really never looked closely into what is happening
> > inside nearest.c, bipolar.c and cubic,c. Could be some obscure bug there,
> > but I will have to get my Grass burning again before I can do any further
> > programming. Right now not possible. Next year when I'm out of work
> > another story.
> >
> > > > But I found out (by trial) that the only usable interpolation method is
> > > > the 'nearest neighbor' (default). The two other methods (bipolar and
> > > > cubic) tend to change ("interpolate") the coastlines and shape of lakes
> > > > too much.
> > >
> > > Shall I disable them to prevent user confusion?

No. Do not disable. Bilinear/Cubic are appropriate for imagery.

O.k, we'll keep it.

Thanks for your comments,
Markus