Updates to add feature for issue #63

Updates to add feature for issue #63

Author: jblindsay

Cargo.lock

@@ -1,7 +1,7 @@
 
 [package]
 name = "whitebox_tools"
-version = "1.1.1"
+version = "1.2.0"
 authors = ["John Lindsay <jlindsay@uoguelph.ca>"]
 description = "A library for analyzing geospatial data."
 keywords = ["geospatial", "GIS", "remote sensing", "geomatics", "image processing", "lidar", "spatial analysis"]

Cargo.toml

@@ -77,6 +77,8 @@ Version 1.1.1 (XX-XX-20XX)
 - The watershed tool now accepts either a set of vector points or a raster for the pour points 
   file. If a raster is specified, all non-zero, non-NoData valued cells will be considered 
   outlet cells and the watershed labels will be assigned based on these values.
+- The D8 and D-infinity flow accumulation tools now take either an input DEM or a flow pointer raster 
+  as inputs.
 
 Version 1.1.0 (09-12-2019)
 - Added the BreachDepressionsLeastCost tool, which performs a modified form of the Lindsay 

readme.txt

@@ -2,7 +2,7 @@
 This tool is part of the WhiteboxTools geospatial analysis library.
 Authors: Dr. John Lindsay
 Created: 26/06/2017
-Last Modified: 21/11/2019
+Last Modified: 21/02/2020
 License: MIT
 */
 
@@ -317,9 +317,15 @@ impl WhiteboxTool for FD8FlowAccumulation {
         let cell_size_y = input.configs.resolution_y;
         let diag_cell_size = (cell_size_x * cell_size_x + cell_size_y * cell_size_y).sqrt();
 
-        // calculate the number of inflowing cells
+        let mut output = Raster::initialize_using_file(&output_file, &input);
+        output.reinitialize_values(1.0);
+        let mut stack = Vec::with_capacity((rows * columns) as usize);
+        let mut num_solved_cells = 0;
+        let mut interior_pit_found = false;
         let mut num_inflowing: Array2D<i8> = Array2D::new(rows, columns, -1, -1)?;
         let num_procs = num_cpus::get() as isize;
+               
+        // calculate the number of inflowing cells
         let (tx, rx) = mpsc::channel();
         for tid in 0..num_procs {
             let input = input.clone();
@@ -328,16 +334,18 @@ impl WhiteboxTool for FD8FlowAccumulation {
                 let d_x = [1, 1, 1, 0, -1, -1, -1, 0];
                 let d_y = [-1, 0, 1, 1, 1, 0, -1, -1];
                 let mut z: f64;
+                let mut zn: f64;
                 let mut count: i8;
                 let mut interior_pit_found = false;
                 for row in (0..rows).filter(|r| r % num_procs == tid) {
                     let mut data: Vec<i8> = vec![-1i8; columns as usize];
                     for col in 0..columns {
-                        z = input[(row, col)];
+                        z = input.get_value(row, col);
                         if z != nodata {
                             count = 0i8;
                             for i in 0..8 {
-                                if input[(row + d_y[i], col + d_x[i])] > z {
+                                zn = input.get_value(row + d_y[i], col + d_x[i]);
+                                if zn > z && zn != nodata {
                                     count += 1;
                                 }
                             }
@@ -352,11 +360,6 @@ impl WhiteboxTool for FD8FlowAccumulation {
             });
         }
 
-        let mut output = Raster::initialize_using_file(&output_file, &input);
-        output.reinitialize_values(1.0);
-        let mut stack = Vec::with_capacity((rows * columns) as usize);
-        let mut num_solved_cells = 0;
-        let mut interior_pit_found = false;
         for r in 0..rows {
             let (row, data, pit) = rx.recv().expect("Error receiving data from thread.");
             num_inflowing.set_row_data(row, data);

src/tools/hydro_analysis/d8_flow_accum.rs

@@ -2,7 +2,7 @@
 This tool is part of the WhiteboxTools geospatial analysis library.
 Authors: Dr. John Lindsay
 Created: 19/02/2020
-Last Modified: 19/02/2020
+Last Modified: 20/02/2020
 License: MIT
 */
 
@@ -251,10 +251,11 @@ impl WhiteboxTool for InsertDams {
         let (mut target_row, mut target_col): (isize, isize);
         let mut profile_intersects_target: bool;
         let mut max_dam_height: f64;
+        let mut dam_z: f64;
         let mut target_cell: usize;
         let (mut dam_row, mut dam_col): (isize, isize);
         let mut dam_dir: usize;
-        let mut best_dam_profile_filled = vec![0f64; dam_profile_length];
+        let mut best_dam_profile_filled: Vec<f64>; // = vec![0f64; dam_profile_length];
         /* Calculate dam heights
         Each cell will be assigned the altitude (ASL) of the highest dam that
         passes through the cell. Potential dams are calculated for each
@@ -264,14 +265,18 @@ impl WhiteboxTool for InsertDams {
             let record = dam_pts.get_record(record_num);
             target_row = input.get_row_from_y(record.points[0].y);
             target_col = input.get_column_from_x(record.points[0].x);
+            dam_z = input.get_value(target_row, target_col);
             dam_row = 0;
             dam_col = 0;
             dam_dir = 0;
             max_dam_height = f64::MIN;
+            // dam_z = f64::MIN;
+            best_dam_profile_filled = vec![0f64; dam_profile_length];
             for row in (target_row-half_dam_length as isize)..=(target_row+half_dam_length as isize) {
                 for col in (target_col-half_dam_length as isize)..=(target_col+half_dam_length as isize) {
                     z = input.get_value(row, col);
                     if z != nodata {
+                        // dam_z = z;
                         for dir in 0..4 {
                             profile_intersects_target = false;
                             target_cell = 0;
@@ -350,7 +355,7 @@ impl WhiteboxTool for InsertDams {
                 }
             }
 
-            if max_dam_height > f64::MIN {
+            if max_dam_height > f64::MIN && max_dam_height > dam_z {
                 // perform the actual damming
                 perp_dir1 = perpendicular1[dam_dir];
                 perp_dir2 = perpendicular2[dam_dir];
@@ -419,7 +424,10 @@ impl WhiteboxTool for InsertDams {
                     }
                 }
             } else {
-                // Error: no dam was found that covered the point
+                // No dam was found that covered the point
+                if verbose {
+                    println!("Warning: No dam could be identified for Point {} due to its position in non-impoundable terrain.", record_num+1);
+                }
             }
 
             if verbose {

src/tools/hydro_analysis/dinf_flow_accum.rs

@@ -382,6 +382,7 @@ def list_tools(self, keywords=[]):
 
 
 
+    
 
 
     ##############
@@ -648,6 +649,20 @@ def raster_to_vector_points(self, i, output, callback=None):
         args.append("--output='{}'".format(output))
         return self.run_tool('raster_to_vector_points', args, callback) # returns 1 if error
 
+    def raster_to_vector_polygons(self, i, output, callback=None):
+        """Converts a raster dataset to a vector of the POLYGON shapetype.
+
+        Keyword arguments:
+
+        i -- Input raster file. 
+        output -- Output vector polygons file. 
+        callback -- Custom function for handling tool text outputs.
+        """
+        args = []
+        args.append("--input='{}'".format(i))
+        args.append("--output='{}'".format(output))
+        return self.run_tool('raster_to_vector_polygons', args, callback) # returns 1 if error
+
     def reinitialize_attribute_table(self, i, callback=None):
         """Reinitializes a vector's attribute table deleting all fields but the feature ID (FID).
 
@@ -3372,24 +3387,28 @@ def burn_streams_at_roads(self, dem, streams, roads, output, width=None, callbac
         if width is not None: args.append("--width='{}'".format(width))
         return self.run_tool('burn_streams_at_roads', args, callback) # returns 1 if error
 
-    def d8_flow_accumulation(self, dem, output, out_type="cells", log=False, clip=False, callback=None):
-        """Calculates a D8 flow accumulation raster from an input DEM.
+    def d8_flow_accumulation(self, i, output, out_type="cells", log=False, clip=False, pntr=False, esri_pntr=False, callback=None):
+        """Calculates a D8 flow accumulation raster from an input DEM or flow pointer.
 
         Keyword arguments:
 
-        dem -- Input raster DEM file. 
+        i -- Input raster DEM or D8 pointer file. 
         output -- Output raster file. 
         out_type -- Output type; one of 'cells' (default), 'catchment area', and 'specific contributing area'. 
         log -- Optional flag to request the output be log-transformed. 
         clip -- Optional flag to request clipping the display max by 1%. 
+        pntr -- Is the input raster a D8 flow pointer rather than a DEM?. 
+        esri_pntr -- Input  D8 pointer uses the ESRI style scheme. 
         callback -- Custom function for handling tool text outputs.
         """
         args = []
-        args.append("--dem='{}'".format(dem))
+        args.append("--input='{}'".format(i))
         args.append("--output='{}'".format(output))
         args.append("--out_type={}".format(out_type))
         if log: args.append("--log")
         if clip: args.append("--clip")
+        if pntr: args.append("--pntr")
+        if esri_pntr: args.append("--esri_pntr")
         return self.run_tool('d8_flow_accumulation', args, callback) # returns 1 if error
 
     def d8_mass_flux(self, dem, loading, efficiency, absorption, output, callback=None):
@@ -3428,26 +3447,28 @@ def d8_pointer(self, dem, output, esri_pntr=False, callback=None):
         if esri_pntr: args.append("--esri_pntr")
         return self.run_tool('d8_pointer', args, callback) # returns 1 if error
 
-    def d_inf_flow_accumulation(self, dem, output, out_type="Specific Contributing Area", threshold=None, log=False, clip=False, callback=None):
+    def d_inf_flow_accumulation(self, i, output, out_type="Specific Contributing Area", threshold=None, log=False, clip=False, pntr=False, callback=None):
         """Calculates a D-infinity flow accumulation raster from an input DEM.
 
         Keyword arguments:
 
-        dem -- Input raster DEM file. 
+        i -- Input raster DEM or D-infinity pointer file. 
         output -- Output raster file. 
         out_type -- Output type; one of 'cells', 'sca' (default), and 'ca'. 
         threshold -- Optional convergence threshold parameter, in grid cells; default is inifinity. 
         log -- Optional flag to request the output be log-transformed. 
         clip -- Optional flag to request clipping the display max by 1%. 
+        pntr -- Is the input raster a D8 flow pointer rather than a DEM?. 
         callback -- Custom function for handling tool text outputs.
         """
         args = []
-        args.append("--dem='{}'".format(dem))
+        args.append("--input='{}'".format(i))
         args.append("--output='{}'".format(output))
         args.append("--out_type={}".format(out_type))
         if threshold is not None: args.append("--threshold='{}'".format(threshold))
         if log: args.append("--log")
         if clip: args.append("--clip")
+        if pntr: args.append("--pntr")
         return self.run_tool('d_inf_flow_accumulation', args, callback) # returns 1 if error
 
     def d_inf_mass_flux(self, dem, loading, efficiency, absorption, output, callback=None):
@@ -3830,6 +3851,24 @@ def impoundment_size_index(self, dem, output, damlength, out_type="depth", callb
         args.append("--damlength='{}'".format(damlength))
         return self.run_tool('impoundment_size_index', args, callback) # returns 1 if error
 
+    def insert_dams(self, dem, dam_pts, output, damlength, callback=None):
+        """Calculates the impoundment size resulting from damming a DEM.
+
+        Keyword arguments:
+
+        dem -- Input raster DEM file. 
+        dam_pts -- Input vector dam points file. 
+        output -- Output file. 
+        damlength -- Maximum length of the dam. 
+        callback -- Custom function for handling tool text outputs.
+        """
+        args = []
+        args.append("--dem='{}'".format(dem))
+        args.append("--dam_pts='{}'".format(dam_pts))
+        args.append("--output='{}'".format(output))
+        args.append("--damlength='{}'".format(damlength))
+        return self.run_tool('insert_dams', args, callback) # returns 1 if error
+
     def isobasins(self, dem, output, size, callback=None):
         """Divides a landscape into nearly equal sized drainage basins (i.e. watersheds).
 
@@ -3894,6 +3933,30 @@ def max_upslope_flowpath_length(self, dem, output, callback=None):
         args.append("--output='{}'".format(output))
         return self.run_tool('max_upslope_flowpath_length', args, callback) # returns 1 if error
 
+    def md_inf_flow_accumulation(self, dem, output, out_type="specific contributing area", exponent=1.1, threshold=None, log=False, clip=False, callback=None):
+        """Calculates an FD8 flow accumulation raster from an input DEM.
+
+        Keyword arguments:
+
+        dem -- Input raster DEM file. 
+        output -- Output raster file. 
+        out_type -- Output type; one of 'cells', 'specific contributing area' (default), and 'catchment area'. 
+        exponent -- Optional exponent parameter; default is 1.1. 
+        threshold -- Optional convergence threshold parameter, in grid cells; default is inifinity. 
+        log -- Optional flag to request the output be log-transformed. 
+        clip -- Optional flag to request clipping the display max by 1%. 
+        callback -- Custom function for handling tool text outputs.
+        """
+        args = []
+        args.append("--dem='{}'".format(dem))
+        args.append("--output='{}'".format(output))
+        args.append("--out_type={}".format(out_type))
+        args.append("--exponent={}".format(exponent))
+        if threshold is not None: args.append("--threshold='{}'".format(threshold))
+        if log: args.append("--log")
+        if clip: args.append("--clip")
+        return self.run_tool('md_inf_flow_accumulation', args, callback) # returns 1 if error
+
     def num_inflowing_neighbours(self, dem, output, callback=None):
         """Computes the number of inflowing neighbours to each cell in an input DEM based on the D8 algorithm.
 
@@ -4092,7 +4155,7 @@ def watershed(self, d8_pntr, pour_pts, output, esri_pntr=False, callback=None):
         Keyword arguments:
 
         d8_pntr -- Input D8 pointer raster file. 
-        pour_pts -- Input vector pour points (outlet) file. 
+        pour_pts -- Input pour points (outlet) file. 
         output -- Output raster file. 
         esri_pntr -- D8 pointer uses the ESRI style scheme. 
         callback -- Custom function for handling tool text outputs.