r.watershed

[a907hart at diamond.sara.nl]

Dan Epstein writes about a problem with r.watershed on a DEC 5100:


> I am running GRASS4.0 on a DEC 5100.  The watershed of concern is very
> large - the elev map is roughly 1700 x 750 cells.  In running
> r.watershed for basin and stream delineation, I dont have enough memory
> to run up front.  As a result, the program resorts to swapping to disk.
>  I managed to free up enough disk space to allow it to do this.
> However, after approx 3 days of running, the program is only 11% done
> with Section 1b- Determining Offmap flow
>
> SECTION 1b (of 5): Determining Offmap Flow. Percent Complete:
> 	11%
>
> At this rate, the program might terminate in 2-3 weeks.
>
> I am running using a min. outside sub-basin area of 150 km2.  The map
> is at roughly 90 meter resolution.
>
> Any suggestions on how to speed this up ?  Or, if I rescale the
> elevation data to 270 meters ( for example), any ideas on how much
> information will be lost ?  At what resolution will the delination
> algorithm not adequately describe the topography ?  Will a 600 X 250
> map run significantly faster ?  Does anyone have experince with scale
> and size in r.watershed ?
>
> Dan Epstein

I had the same problem with r.watershed on an RS/6000 with a slightly
smaller mapset. The problem did not disappear when I reran it on the
mainframe with about 750M internal memory. After some
puzzling I found that there is a fault in memory allocation resulting
from am macro expansion in "ramseg.h", in init_vars.c, if I remember well. I
rewrote that macro as a function (don't forget to #define MAIN in main.c)
and all went well and in negligable time (5-10 minutes).
I enclose this function with the original macro commented out.


#define RAMSEG		int
#define RAMSEGBITS 	4
#define DOUBLEBITS 	8	/* 2 * ramsegbits	*/
#define SEGLENLESS 	15	/* 2 ^ ramsegbits - 1	*/
/*
#define SEG_INDEX(s,r,c) (int) \
   ((((r) >> RAMSEGBITS) * (s) + (((c) >> RAMSEGBITS)) << DOUBLEBITS) \
    + (((r) & SEGLENLESS) << RAMSEGBITS) + ((c) & SEGLENLESS))
*/

#ifdef MAIN
SEG_INDEX(s,r,c)
int s,r,c;
{	return ((((r) >> RAMSEGBITS) * (s) + (((c) >> RAMSEGBITS)) << DOUBLEBITS)
		+ (((r) & SEGLENLESS) << RAMSEGBITS) + ((c) & SEGLENLESS));
}
#endif

Michael Shapiro blamed the IBM prepocessor some time ago on this list.
This could be very well the case, (I don't understand what happens here
at all), but I am not able to test this out further, as this work
on watershed was just a quick hack to get some students working on a particular
problem. Since then we haven't worked with  watershed analysis and the
dataset isn't available any more. Perhaps Dan Epstein's dataset
and a watershed script could be made available at moon for testing on
different architectures?


Jan Hartmann
Faculty of Environmental Sciences
University of Amsterdam