r.terraflow: allow larger maps (#265)

* allow larger maps; patch provided by @metzm
* add testsuite

See: https://trac.osgeo.org/grass/ticket/1421

raster/r.terraflow/3scan.h

@@ -193,7 +193,7 @@ scan3(AMI_STREAM<T> &amis0,
     l_prev = l_crt;
     l_crt = l_next;
     if (i < nr-2) {
-      ae = amis0.new_substream(AMI_READ_STREAM, (i+2)*nc, (i+3)*nc-1, 
+      ae = amis0.new_substream(AMI_READ_STREAM, (off_t)(i+2)*nc, (off_t)(i+3)*nc-1,
 			      &l_next);
       assert(ae == AMI_ERROR_NO_ERROR);
     } else {
@@ -258,7 +258,7 @@ memoryScan(AMI_STREAM<T> &str,
   str.seek(0);
   
   assert(nrows > 1);
-  assert(nrows * ncols == str.stream_len());
+  assert((off_t)nrows * ncols == str.stream_len());
   buf[0] = new T[ncols+2];
   buf[1] = new T[ncols+2];
   buf[2] = new T[ncols+2];
@@ -310,8 +310,8 @@ memoryScan(AMI_STREAM<T1> &str1, AMI_STREAM<T2> &str2,
   str2.seek(0);
 
   assert(nrows > 1);
-  assert(nrows * ncols == str1.stream_len());
-  assert(nrows * ncols == str2.stream_len());
+  assert((off_t)nrows * ncols == str1.stream_len());
+  assert((off_t)nrows * ncols == str2.stream_len());
   buf1[0] = new T1[ncols+2];
   buf1[1] = new T1[ncols+2];
   buf1[2] = new T1[ncols+2];

raster/r.terraflow/ccforest.cpp

@@ -79,8 +79,8 @@ ccforest<T>::~ccforest() {
 /* ---------------------------------------------------------------------- */
 
 template<class T>
-int ccforest<T>::size() {
-  size_t streamLength = edgeStream->stream_len();
+off_t ccforest<T>::size() {
+  off_t streamLength = edgeStream->stream_len();
   return streamLength;
 }
 
@@ -134,11 +134,11 @@ void ccforest<T>::findAllRoots(int depth) {
   EMPQueueAdaptive<cckeyvalue,T> *pq =
 	new EMPQueueAdaptive<cckeyvalue,T>();	/* parent queue */
 
-  size_t streamLength = edgeStream->stream_len();
+  off_t streamLength = edgeStream->stream_len();
   T prevSrc = T(-1);
   T parent = T(-1);
   ccedge prevEdge;
-  for(unsigned int i=0; i<streamLength; i++) {
+  for(off_t i=0; i<streamLength; i++) {
 	ccedge *e;
 	AMI_err ae = edgeStream->read_item(&e);
 	assert(ae == AMI_ERROR_NO_ERROR);

raster/r.terraflow/ccforest.h

@@ -139,7 +139,7 @@ public:
   T findNextRoot(const T& i);	/* find root where i >= prev i */
   void printRootStream();
   void printEdgeStream();
-  int size();
+  off_t size();
 };
 
 

raster/r.terraflow/fill.cpp

@@ -275,7 +275,7 @@ computeFlowDirections(AMI_STREAM<elevation_type>*& elstr,
   assert(boundaryStr->stream_len() == 0);
   delete boundaryStr;
   
-  assert(labeledWater->stream_len() == nrows * ncols);
+  assert(labeledWater->stream_len() == (off_t)nrows * ncols);
   rt_stop(rt);
   if (stats)
     stats->recordTime("generating watersheds", rt);
@@ -315,7 +315,7 @@ computeFlowDirections(AMI_STREAM<elevation_type>*& elstr,
   rt_start(rt);
   filledstr = new AMI_STREAM<elevation_type>();
   commit_fill(labeledWater, raise, labelFactory::getLabelCount(), filledstr);
-  assert(filledstr->stream_len() == nrows * ncols);
+  assert(filledstr->stream_len() == (off_t)nrows * ncols);
   delete [] raise;
   rt_stop(rt);
   if (stats) {
@@ -653,9 +653,9 @@ mergeStreamGridGrid(AMI_STREAM<T1> *grid1,
       ae = outStream->write_item(t4);
       assert(ae == AMI_ERROR_NO_ERROR);
     }
-    /*assert(outStream->stream_len() == (row+1) * cols); */
+    /*assert(outStream->stream_len() == (off_t)(row+1) * cols); */
   }
-  assert(outStream->stream_len() == rows * cols);
+  assert(outStream->stream_len() == (off_t)rows * cols);
   return;
 }
 

raster/r.terraflow/flow.cpp

@@ -210,7 +210,7 @@ fillstr2sweepstr(AMI_STREAM<waterWindowBaseType>* fillStream) {
   if (stats)
     stats->comment("creating sweep stream from fill output stream");
   
-  assert(fillStream->stream_len() == nrows * ncols);
+  assert(fillStream->stream_len() == (off_t)nrows * ncols);
   
   /* create the sweep stream */
   sweepstr = new AMI_STREAM<sweepItem>();

raster/r.terraflow/grass2str.h

@@ -135,9 +135,9 @@ cell2stream(char* cellname, elevation_type T_max_value, long* nodata_count) {
   /* close map files */
   Rast_close (infd);
 
-  G_debug(3, "nrows=%d   ncols=%d    stream_len()=%" PRI_OFF_T, nrows, ncols,
+  G_debug(1, "nrows=%d   ncols=%d    stream_len()=%" PRI_OFF_T, nrows, ncols,
 		str->stream_len());  
-  assert((off_t) nrows * ncols == str->stream_len());
+  assert((off_t)nrows * ncols == str->stream_len());
   rt_stop(rt);
   if (stats)
     stats->recordTime("reading raster map", rt);
@@ -160,7 +160,7 @@ stream2_CELL(AMI_STREAM<T>* str, dimension_type nrows, dimension_type ncols,
   
   rt_start(rt); 
   assert(str);
-  assert(str->stream_len() == nrows*ncols);
+  assert(str->stream_len() == (off_t)nrows*ncols);
   str->seek(0);
   {
     char * foo;
@@ -242,7 +242,7 @@ stream2_CELL(AMI_STREAM<T> *str, dimension_type nrows, dimension_type ncols,
   AMI_err ae; 
   
   assert(str && cellname);
-  /* assert(str->stream_len() == nrows*ncols); */
+  /* assert(str->stream_len() == (off_t)nrows*ncols); */
 
   rt_start(rt); 
   str->seek(0);
@@ -316,7 +316,7 @@ stream2_FCELL(AMI_STREAM<T> *str, dimension_type nrows, dimension_type ncols,
   AMI_err ae; 
   
   assert(str && cellname);
-  /* assert(str->stream_len() == nrows*ncols); */
+  /* assert(str->stream_len() == (off_t)nrows*ncols); */
 
   rt_start(rt); 
   str->seek(0);

raster/r.terraflow/grid.cpp

@@ -31,10 +31,10 @@ grid::grid(dimension_type giMin, dimension_type gjMin,
   iMin(giMin-1), jMin(gjMin-1), label(glabel), size(gsize) {
   width = jMax - jMin + 2;
   height = iMax - iMin + 2;
-  assert(width*height*sizeof(gridElement) < getAvailableMemory());
-  data = new gridElement[width*height];
+  assert((size_t)width*height*sizeof(gridElement) < getAvailableMemory());
+  data = new gridElement[(size_t)width*height];
   assert(data);
-  memset(data, 0, width*height*sizeof(gridElement));
+  memset(data, 0, (size_t)width*height*sizeof(gridElement));
 }
 
 
@@ -60,7 +60,7 @@ grid::load(AMI_STREAM<plateauType> &str) {
 	dimension_type pti, ptj;
 	pti = pt->i - iMin;
 	ptj = pt->j - jMin;
-	gridElement *datap = data + pti * width + ptj;
+	gridElement *datap = data + (size_t)pti * width + ptj;
 	datap->dir = pt->dir;
 	datap->depth = DEPTH_INITIAL;			/* initial depth */
 	datap->valid = 1;
@@ -80,7 +80,7 @@ grid::save(AMI_STREAM<waterType> &str) {
   GRID_DEBUG cout << "saving grid" << endl;
 
   for(dimension_type i=1; i<height-1; i++) {
-    gridElement *rowp = data + i * width;
+    gridElement *rowp = data + (size_t)i * width;
     for(dimension_type j=1; j<width-1; j++) {
       gridElement *datap = rowp + j;
       if(datap->valid) {
@@ -104,8 +104,8 @@ grid::print() {
   for(int j=0; j<height; j++) {
     printf("%3d ", j + iMin);
     for(int i=0; i<width; i++) {
-      if(data[i+width*j].valid) {
-	cout << " " << directionSymbol(data[i+width*j].dir);
+      if(data[i + (size_t)width * j].valid) {
+	cout << " " << directionSymbol(data[i + (size_t)width * j].dir);
       } else {
 	cout << " .";
       }

raster/r.terraflow/main.cpp

@@ -413,20 +413,20 @@ setSinkWatershedColorTable(char* cellname) {
 void
 printMaxSortSize(long nodata_count) {
   char buf[BUFSIZ];
-  long long  fillmaxsize = (long long)nrows*ncols*sizeof(waterWindowType);
-  long long  flowmaxsize = (long long)(nrows*ncols - nodata_count)*sizeof(sweepItem);
-  long long maxneed = (fillmaxsize > flowmaxsize) ? fillmaxsize: flowmaxsize;
+  off_t  fillmaxsize = (off_t)nrows*ncols*sizeof(waterWindowType);
+  off_t  flowmaxsize = ((off_t)nrows*ncols - nodata_count)*sizeof(sweepItem);
+  off_t maxneed = (fillmaxsize > flowmaxsize) ? fillmaxsize: flowmaxsize;
   maxneed =  2*maxneed; /* need 2*N to sort */
 
-  G_debug(1, "total elements=%ld, nodata elements=%ld",
-          (long)nrows * ncols, nodata_count);
+  G_debug(1, "total elements=%lld, nodata elements=%ld",
+          (off_t)nrows * ncols, nodata_count);
   G_debug(1, "largest temporary files: ");
-  G_debug(1, "\t\t FILL: %s [%ld elements, %ldB each]",
+  G_debug(1, "\t\t FILL: %s [%lld elements, %ldB each]",
           formatNumber(buf, fillmaxsize),
-          (long)nrows * ncols, sizeof(waterWindowType));
-  G_debug(1, "\t\t FLOW: %s [%ld elements, %ldB each]",
+          (off_t)nrows * ncols, sizeof(waterWindowType));
+  G_debug(1, "\t\t FLOW: %s [%lld elements, %ldB each]",
           formatNumber(buf, flowmaxsize),
-          (long)nrows * ncols - nodata_count, sizeof(sweepItem));
+          (off_t)nrows * ncols - nodata_count, sizeof(sweepItem));
   G_debug(1, "Will need at least %s space available in %s",
           formatNumber(buf, maxneed),  	  /* need 2*N to sort */
           getenv(STREAM_TMPDIR));
@@ -514,7 +514,7 @@ main(int argc, char *argv[]) {
       record_args(argc, argv);
       {
 	char buf[BUFSIZ];
-	long grid_size = nrows * ncols;
+	off_t grid_size = (off_t) nrows * ncols;
 	*stats << "region size = " <<  formatNumber(buf, grid_size) << " elts "
 	       << "(" << nrows << " rows x " << ncols << " cols)\n";
 
@@ -589,7 +589,7 @@ main(int argc, char *argv[]) {
 
   sprintf(path, "%s/flowStream", streamdir->answer);
   flowStream = new AMI_STREAM<waterWindowBaseType>(path);
-  G_verbose_message(_("flowStream opened: len=%d\n", flowStream->stream_len());
+  G_verbose_message(_("flowStream opened: len=%lld\n", flowStream->stream_len());
   G_verbose_message(_("jumping to flow accumulation computation\n");
 #endif
   

raster/r.terraflow/streamutils.h

@@ -96,7 +96,7 @@ void printGridStream(AMI_STREAM<T> *str,
   fstrm << "rows=" << nrows << endl;
   fstrm << "cols=" << ncols << endl;
 
-  assert(str->stream_len() == nrows * ncols);
+  assert(str->stream_len() == (off_t)nrows * ncols);
   str->seek(0);
   for(dimension_type i=0; i<nrows; i++) {
 	for(dimension_type j=0; j<ncols; j++) {

raster/r.terraflow/sweep.cpp

@@ -146,7 +146,8 @@ initializePQ() {
 #ifdef EM_PQUEUE
   if (stats)
     stats->comment("FLOW_DATASTRUCTURE: ext-memory pqueue");
-  flowpq = new FLOW_DATASTR(nrows * ncols);  
+  /* TODO: long -> off_t in include/iostream/empq.h */
+  flowpq = new FLOW_DATASTR((long)nrows * ncols);
 #endif
 #ifdef EMPQ_ADAPTIVE
   if (opt->verbose && stats) stats->comment("FLOW_DATASTRUCTURE: adaptive pqueue");
@@ -170,7 +171,7 @@ sweep(AMI_STREAM<sweepItem> *sweepstr, const flowaccumulation_type D8CUT,
   sweepItem* crtpoint;
   AMI_err ae;
   flowStructure x;	
-  long nitems;
+  off_t nitems;
   Rtimer rt;
   AMI_STREAM<sweepOutput>* outstr;
 
@@ -200,7 +201,7 @@ sweep(AMI_STREAM<sweepItem> *sweepstr, const flowaccumulation_type D8CUT,
   ae = sweepstr->seek(0);
   assert(ae == AMI_ERROR_NO_ERROR);
   G_important_message(_("Sweeping..."));
-  for (long k = 0; k < nitems; k++) {
+  for (off_t k = 0; k < nitems; k++) {
     
     /* cout << k << endl; cout.flush(); */
     /* read next sweepItem = (prio, elevwin, topoRankwin, dir) */

raster/r.terraflow/testsuite/test_r_terraflow.py

@@ -0,0 +1,117 @@
+"""Test check stability of results from r.terraflow
+@author Stefan Blumentrath, NINA
+"""
+
+import os
+import tempfile
+from grass.gunittest.case import TestCase
+
+class TestTerraflow(TestCase):
+
+    elevation = 'elevation'
+    testdir = os.path.join(tempfile.gettempdir(), 'terraflow_test')
+    teststats = os.path.join(tempfile.gettempdir(), 'terraflow_test_stats.txt')
+  
+    @classmethod
+    def setUpClass(cls):
+        """Use temporary region settings"""
+        cls.use_temp_region()
+        cls.runModule("g.region", flags="p",  raster=cls.elevation)
+
+    @classmethod
+    def tearDownClass(cls):
+        """!Remove the temporary region
+        """
+        cls.del_temp_region()
+
+    def setUp(cls):
+        """Create input data for steady state groundwater flow computation
+        """
+        if not os.path.exists(cls.testdir):
+            os.mkdir(cls.testdir)
+
+    def test_univar_mfd(cls):
+        #compute a steady state groundwater flow
+        cls.assertModule("r.terraflow", overwrite=True, verbose=True,
+        elevation=cls.elevation, filled='terra_flooded',
+        direction='terra_flowdir', swatershed='terra_sink',
+        accumulation='terra_flowaccum', tci='terra_tci',
+        directory=cls.testdir, stats=cls.teststats)
+        
+        # Output of r.univar -g
+        terra_flooded_univar="""n=2025000
+null_cells=0
+cells=2025000
+min=55.5787925720215
+max=156.329864501953
+range=100.751071929932
+mean=110.466900511132
+mean_of_abs=110.466900511132
+stddev=20.2568412316924
+variance=410.339616685993
+coeff_var=18.3374758755462
+sum=223695473.535042"""
+
+        terra_flowdir_univar="""n=2025000
+null_cells=0
+cells=2025000
+min=0
+max=255
+range=255
+mean=114.239481481481
+mean_of_abs=114.239481481481
+stddev=84.9304048144913
+variance=7213.17366195336
+coeff_var=74.344179186649
+sum=231334950"""
+        
+        terra_sink_univar="""n=2025000
+null_cells=0
+cells=2025000
+min=0
+max=7945
+range=7945
+mean=3716.98878864198
+mean_of_abs=3716.98878864198
+stddev=2352.78190064133
+variance=5535582.67198542
+coeff_var=63.2980628790363
+sum=7526902297"""
+
+        terra_flowaccum_univar="""n=2025000
+null_cells=0
+cells=2025000
+min=1
+max=638570.4375
+range=638569.4375
+mean=644.701550164795
+mean_of_abs=644.701550164795
+stddev=10616.1468932394
+variance=112702574.858836
+coeff_var=1646.67618536449
+sum=1305520639.08371"""
+
+        terra_tci_univar="""n=2025000
+null_cells=0
+cells=2025000
+min=1.07463788986206
+max=16.7091903686523
+range=15.6345524787903
+mean=4.11934358476421
+mean_of_abs=4.11934358476421
+stddev=1.97140337926634
+variance=3.88643128378274
+coeff_var=47.8572213922083
+sum=8341670.75914752"""
+
+        #cls.assertRasterFitsUnivar(raster="terra_flooded",  reference=terra_flooded_univar,  precision=3)
+        cls.assertRasterFitsUnivar(raster="terra_flowdir",  reference=terra_flowdir_univar,  precision=3)
+        cls.assertRasterFitsUnivar(raster="terra_sink",  reference=terra_sink_univar,  precision=3)
+        cls.assertRasterFitsUnivar(raster="terra_flowaccum",  reference=terra_flowaccum_univar,  precision=3)
+        cls.assertRasterFitsUnivar(raster="terra_tci",  reference=terra_tci_univar,  precision=3)
+
+if __name__ == '__main__':
+    from grass.gunittest.main import test
+    test()
+
+

raster/r.terraflow/types.h

@@ -26,8 +26,8 @@
 /* input parameters type */
 /* ------------------------------------------------------------ */
 
-typedef short dimension_type; /* represent dimension of the grid */
-static const dimension_type dimension_type_max=SHRT_MAX;
+typedef int dimension_type; /* represent dimension of the grid */
+static const dimension_type dimension_type_max=INT_MAX;
 
 typedef short direction_type;  /* represent the direction of a cell */
 static const direction_type DIRECTION_UNDEF=-1;

raster/r.terraflow/water.cpp

@@ -256,7 +256,7 @@ generateWatersheds(AMI_STREAM<waterWindowType> **waterWindows,
   if (stats)
     stats->comment("generateWatersheds", opt->verbose);
 
-  assert((*waterWindows)->stream_len() == (nrows * ncols));
+  assert((*waterWindows)->stream_len() == ((off_t)nrows * ncols));
 
   WATER_DEBUG cout << "sort waterWindowsStream (by priority): ";
   sort(waterWindows, priorityCmpWaterWindowType());
@@ -273,7 +273,7 @@ generateWatersheds(AMI_STREAM<waterWindowType> **waterWindows,
   if (stats)
     stats->comment("starting generate watersheds main loop", opt->verbose);
   
-  assert((*waterWindows)->stream_len() == (nrows * ncols));
+  assert((*waterWindows)->stream_len() == ((off_t)nrows * ncols));
   /* not really in a grid, so row, col are not valid (but count correct) */
   for(dimension_type row=0; row<nrows; row++) {
     for(dimension_type col=0; col<ncols; col++) {