Date: Tue, 29 Jul 2008 20:46:02 +0200
From: Markus Metz <[EMAIL PROTECTED]>
Subject: Re: [GRASS-user] r.watershed speed-up
To: Charles Ehlschlaeger <[EMAIL PROTECTED]>
Cc: [EMAIL PROTECTED]
Message-ID: <[EMAIL PROTECTED]>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Dear Chuck,
r.watershed is a much valued tool in GRASS, for me the best watershed
analsis tool not only in GRASS, therefore I thought about a a way to
keep the results identical too. I am also aware that the closer the
results produced by changes in the algorithm are to the results
produced
by original algorithm, the higher the chances that it will be accepted
by the community.
With regard to your suggestion, I would not adjust DEM values, because
in larger regions the minimum possible increment is already there in
the
data, i.e. there are no gaps in the data distribution that can be
filled
with adjusted values. One theoretical way out would be to read in DEMs
as FCELL or DCELL, but then there is the floating point comparison
problem. (I tried against better knowledge, it doesn't work).
Regarding
the breadth first search, where do you see breadth first <when the DEM
values are different>? You lost me there. I don't see differences in
how
points are searched between the two versions, but maybe I have not
fully
understood the original algorithm. As far as I have understood the
original algorithm, the list of astar_pts following astar_pts.nxt is
kept in ascending order using elevation. If there are already points
with equal elevation, the new point is inserted after all other points
with the same elevation (line 91 in original do_astar.c), so that the
point inserted first (of several points with equal elevation) will be
removed first (line 19 in original do_astar.c). This is still the case
in the new algorithm (insertion: line 136, removal: lines 31 and 192)
most of the time. If the binary heap becomes fairly large and there
are
many points with equal elevation, there might be an exception. Please
let me know if I got something wrong there!
Another possibility to produce the exact same results like in the
original version would be to go recursively down the heap and pick the
point added earliest from all points with elevation equal to the root
point. This is easy to implement, but it would have slowed down the
search algorithm somewhat and I wanted to get something lightening
fast.
I have one main argument why it is not a disaster if the results are
not
100% identical:
The order in which neighbouring cells are added is in both versions,
with respect to the focus cell:
low, up, left, right, upper right corner, lower left corner, lower
right
corner, upper left corner
This order is always kept, irrespective of the already established
flow
direction, thus it is a random order and there is not really a reason
why the algorithm should stick to that order. I think a rare
replacement
of that random order (2% difference of flow direction in Moritz
Lennert's test) with another random order (binary heap shuffling) is
not
a disaster and the result is still valid. I did build in a check to
make
results more similar, but there are still scenarios when this check
doesn't catch.
So my main question to you, the original author of r.watershed, is,
if a
rare violation to the (in my opinion random) order in which neighbours
are added to the list would cause the results to be no longer valid.
The other question is if I should provide now a version that really
produces identical results, or if I first sort out the problem of how
neighbours are (should be) added to and removed from the list. BTW, I
tried to change the order of adding neighbours to the list too, taking
into account the already established flow direction. It produces very
straight lines in flat terrain, which is ok in hydrological terms, but
some randomness looked better. Flat terrain in the DEM must not be
flat
in reality because of problems with DEM resolution and accuracy,
randomness produces there more naturally looking results.
Sorry for the long reply!
Regards,
Markus
