diff --git a/Cargo.lock b/Cargo.lock index 927ca34eb..660113392 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -569,7 +569,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" [[package]] name = "whitebox_tools" -version = "1.1.1" +version = "1.2.0" dependencies = [ "byteorder 1.3.1 (registry+https://github.com/rust-lang/crates.io-index)", "chrono 0.4.6 (registry+https://github.com/rust-lang/crates.io-index)", diff --git a/Cargo.toml b/Cargo.toml index cbfe097ce..7105d2eda 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,7 +1,7 @@ [package] name = "whitebox_tools" -version = "1.1.1" +version = "1.2.0" authors = ["John Lindsay "] description = "A library for analyzing geospatial data." keywords = ["geospatial", "GIS", "remote sensing", "geomatics", "image processing", "lidar", "spatial analysis"] diff --git a/readme.txt b/readme.txt index f3c0fc264..90830d73c 100644 --- a/readme.txt +++ b/readme.txt @@ -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 diff --git a/src/tools/hydro_analysis/d8_flow_accum.rs b/src/tools/hydro_analysis/d8_flow_accum.rs index fdffca7b3..1acbbddc8 100644 --- a/src/tools/hydro_analysis/d8_flow_accum.rs +++ b/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: 26/016/2017 -Last Modified: 18/10/2019 +Last Modified: 21/02/2020 License: MIT */ @@ -22,11 +22,14 @@ use std::thread; /// D8 (O'Callaghan and Mark, 1984) algorithm. This algorithm is an example of single-flow-direction /// (SFD) method because the flow entering each grid cell is routed to only one downslope neighbour, /// i.e. flow divergence is not permitted. The user must specify the name of the input digital -/// elevation model (DEM). The DEM must have been hydrologically corrected to remove all spurious -/// depressions and flat areas. DEM pre-processing is usually achieved using the `BreachDepressions` or -/// `FillDepressions` tools. +/// elevation model (DEM) or D8 flow pointer (`DInfPointer`) raster (`--input`). If an input DEM is used, it must have +/// been hydrologically corrected to remove all spurious depressions and flat areas. DEM pre-processing +/// is usually achieved using the `BreachDepressionsLeastCost` or `FillDepressions` tools. If a D8 pointer +/// raster is input, the user must also specify the optional `--pntr` flag. If the D8 pointer follows +/// the Esri pointer scheme, rather than the default WhiteboxTools scheme, the user must also specify the +/// optional `--esri_pntr` flag. /// -/// In addition to the input DEM, the user must specify the output type. The output flow-accumulation +/// In addition to the input DEM/pointer, the user must specify the output type. The output flow-accumulation /// can be 1) `cells` (i.e. the number of inflowing grid cells), `catchment area` (i.e. the upslope area), /// or `specific contributing area` (i.e. the catchment area divided by the flow width. The default value /// is `cells`. The user must also specify whether the output flow-accumulation grid should be @@ -39,11 +42,11 @@ use std::thread; /// however, log-transformed flow-accumulation grids must not be used to estimate other secondary terrain /// indices, such as the wetness index, or relative stream power index. /// -/// Grid cells possessing the **NoData** value in the input flow-pointer grid are assigned the **NoData** +/// Grid cells possessing the **NoData** value in the input DEM/pointer raster are assigned the **NoData** /// value in the output flow-accumulation image. /// /// # See Also: -/// `DInfFlowAccumulation`, `BreachDepressions`, `FillDepressions` +/// `DInfPointer`, `DInfFlowAccumulation`, `BreachDepressionsLeastCost`, `FillDepressions` pub struct D8FlowAccumulation { name: String, description: String, @@ -57,13 +60,13 @@ impl D8FlowAccumulation { // public constructor let name = "D8FlowAccumulation".to_string(); let toolbox = "Hydrological Analysis".to_string(); - let description = "Calculates a D8 flow accumulation raster from an input DEM.".to_string(); + let description = "Calculates a D8 flow accumulation raster from an input DEM or flow pointer.".to_string(); let mut parameters = vec![]; parameters.push(ToolParameter { - name: "Input DEM File".to_owned(), - flags: vec!["-i".to_owned(), "--dem".to_owned()], - description: "Input raster DEM file.".to_owned(), + name: "Input DEM or D8 Pointer File".to_owned(), + flags: vec!["-i".to_owned(), "--input".to_owned()], + description: "Input raster DEM or D8 pointer file.".to_owned(), parameter_type: ParameterType::ExistingFile(ParameterFileType::Raster), default_value: None, optional: false, @@ -105,6 +108,24 @@ impl D8FlowAccumulation { optional: true, }); + parameters.push(ToolParameter { + name: "Is the input raster a D8 flow pointer?".to_owned(), + flags: vec!["--pntr".to_owned()], + description: "Is the input raster a D8 flow pointer rather than a DEM?".to_owned(), + parameter_type: ParameterType::Boolean, + default_value: None, + optional: true, + }); + + parameters.push(ToolParameter { + name: "If a pointer is input, does it use the ESRI pointer scheme?".to_owned(), + flags: vec!["--esri_pntr".to_owned()], + description: "Input D8 pointer uses the ESRI style scheme.".to_owned(), + parameter_type: ParameterType::Boolean, + default_value: Some("false".to_owned()), + optional: true, + }); + let sep: String = path::MAIN_SEPARATOR.to_string(); let p = format!("{}", env::current_dir().unwrap().display()); let e = format!("{}", env::current_exe().unwrap().display()); @@ -116,8 +137,8 @@ impl D8FlowAccumulation { if e.contains(".exe") { short_exe += ".exe"; } - let usage = format!(">>.*{0} -r={1} -v --wd=\"*path*to*data*\" --dem=DEM.tif -o=output.tif --out_type='cells' ->>.*{0} -r={1} -v --wd=\"*path*to*data*\" --dem=DEM.tif -o=output.tif --out_type='specific catchment area' --log --clip", short_exe, name).replace("*", &sep); + let usage = format!(">>.*{0} -r={1} -v --wd=\"*path*to*data*\" --input=DEM.tif -o=output.tif --out_type='cells' +>>.*{0} -r={1} -v --wd=\"*path*to*data*\" --input=DEM.tif -o=output.tif --out_type='specific catchment area' --log --clip", short_exe, name).replace("*", &sep); D8FlowAccumulation { name: name, @@ -168,6 +189,8 @@ impl WhiteboxTool for D8FlowAccumulation { let mut out_type = String::from("sca"); let mut log_transform = false; let mut clip_max = false; + let mut pntr_input = false; + let mut esri_style = false; if args.len() == 0 { return Err(Error::new( @@ -184,23 +207,20 @@ impl WhiteboxTool for D8FlowAccumulation { if vec.len() > 1 { keyval = true; } - if vec[0].to_lowercase() == "-i" - || vec[0].to_lowercase() == "--input" - || vec[0].to_lowercase() == "--dem" - { + let flag_val = vec[0].to_lowercase().replace("--", "-"); + if flag_val == "-i" || flag_val == "-input" || flag_val == "-dem" { if keyval { input_file = vec[1].to_string(); } else { input_file = args[i + 1].to_string(); } - } else if vec[0].to_lowercase() == "-o" || vec[0].to_lowercase() == "--output" { + } else if flag_val == "-o" || flag_val == "-output" { if keyval { output_file = vec[1].to_string(); } else { output_file = args[i + 1].to_string(); } - } else if vec[0].to_lowercase() == "-out_type" || vec[0].to_lowercase() == "--out_type" - { + } else if flag_val == "-out_type" { if keyval { out_type = vec[1].to_lowercase(); } else { @@ -213,14 +233,23 @@ impl WhiteboxTool for D8FlowAccumulation { } else { out_type = String::from("ca"); } - } else if vec[0].to_lowercase() == "-log" || vec[0].to_lowercase() == "--log" { + } else if flag_val == "-log" { if vec.len() == 1 || !vec[1].to_string().to_lowercase().contains("false") { log_transform = true; } - } else if vec[0].to_lowercase() == "-clip" || vec[0].to_lowercase() == "--clip" { + } else if flag_val == "-clip" { if vec.len() == 1 || !vec[1].to_string().to_lowercase().contains("false") { clip_max = true; } + } else if flag_val == "-pntr" { + if vec.len() == 1 || !vec[1].to_string().to_lowercase().contains("false") { + pntr_input = true; + } + } else if flag_val == "-esri_pntr" || flag_val == "-esri_style" { + if vec.len() == 1 || !vec[1].to_string().to_lowercase().contains("false") { + esri_style = true; + pntr_input = true; + } } } @@ -248,7 +277,6 @@ impl WhiteboxTool for D8FlowAccumulation { let input = Arc::new(Raster::new(&input_file, "r")?); - // calculate the flow direction let start = Instant::now(); let rows = input.configs.rows as isize; let columns = input.configs.columns as isize; @@ -257,108 +285,198 @@ impl WhiteboxTool for D8FlowAccumulation { let cell_size_x = input.configs.resolution_x; 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(); - - let mut flow_dir: Array2D = Array2D::new(rows, columns, -1, -1)?; + // -2 indicates NoData, -1 indicates no downslope neighbour, 0-7 indicate flow to one neighbour. + let mut flow_dir: Array2D = Array2D::new(rows, columns, -2, -2)?; + let mut interior_pit_found = false; let num_procs = num_cpus::get() as isize; - let (tx, rx) = mpsc::channel(); - for tid in 0..num_procs { - let input = input.clone(); - let tx = tx.clone(); - thread::spawn(move || { - let nodata = input.configs.nodata; - let dx = [1, 1, 1, 0, -1, -1, -1, 0]; - let dy = [-1, 0, 1, 1, 1, 0, -1, -1]; - let grid_lengths = [ - diag_cell_size, - cell_size_x, - diag_cell_size, - cell_size_y, - diag_cell_size, - cell_size_x, - diag_cell_size, - cell_size_y, - ]; - let (mut z, mut z_n): (f64, f64); - let (mut max_slope, mut slope): (f64, f64); - let mut dir: i8; - let mut neighbouring_nodata: bool; - let mut interior_pit_found = false; - for row in (0..rows).filter(|r| r % num_procs == tid) { - let mut data: Vec = vec![-1i8; columns as usize]; - for col in 0..columns { - z = input[(row, col)]; - if z != nodata { - dir = 0i8; - max_slope = f64::MIN; - neighbouring_nodata = false; - for i in 0..8 { - z_n = input[(row + dy[i], col + dx[i])]; - if z_n != nodata { - slope = (z - z_n) / grid_lengths[i]; - if slope > max_slope && slope > 0f64 { - max_slope = slope; - dir = i as i8; + + if !pntr_input { + // calculate the flow direction from the input DEM + let (tx, rx) = mpsc::channel(); + for tid in 0..num_procs { + let input = input.clone(); + let tx = tx.clone(); + thread::spawn(move || { + let nodata = input.configs.nodata; + let dx = [1, 1, 1, 0, -1, -1, -1, 0]; + let dy = [-1, 0, 1, 1, 1, 0, -1, -1]; + let grid_lengths = [ + diag_cell_size, + cell_size_x, + diag_cell_size, + cell_size_y, + diag_cell_size, + cell_size_x, + diag_cell_size, + cell_size_y, + ]; + let (mut z, mut z_n): (f64, f64); + let (mut max_slope, mut slope): (f64, f64); + let mut dir: i8; + let mut neighbouring_nodata: bool; + let mut interior_pit_found = false; + for row in (0..rows).filter(|r| r % num_procs == tid) { + let mut data: Vec = vec![-2i8; columns as usize]; + for col in 0..columns { + z = input.get_value(row, col); + if z != nodata { + dir = 0i8; + max_slope = f64::MIN; + neighbouring_nodata = false; + for i in 0..8 { + z_n = input[(row + dy[i], col + dx[i])]; + if z_n != nodata { + slope = (z - z_n) / grid_lengths[i]; + if slope > max_slope && slope > 0f64 { + max_slope = slope; + dir = i as i8; + } + } else { + neighbouring_nodata = true; } + } + if max_slope >= 0f64 { + data[col as usize] = dir; } else { - neighbouring_nodata = true; + data[col as usize] = -1i8; + if !neighbouring_nodata { + interior_pit_found = true; + } } } - if max_slope >= 0f64 { - data[col as usize] = dir; - } else { - data[col as usize] = -1i8; - if !neighbouring_nodata { - interior_pit_found = true; + } + tx.send((row, data, interior_pit_found)).unwrap(); + } + }); + } + + for r in 0..rows { + let (row, data, pit) = rx.recv().expect("Error receiving data from thread."); + flow_dir.set_row_data(row, data); + if pit { + interior_pit_found = true; + } + if verbose { + progress = (100.0_f64 * r as f64 / (rows - 1) as f64) as usize; + if progress != old_progress { + println!("Flow directions: {}%", progress); + old_progress = progress; + } + } + } + } else { // The input raster is a D8 flow pointer + // map the pointer values into 0-7 style pointer vlaues + let (tx, rx) = mpsc::channel(); + for tid in 0..num_procs { + let input = input.clone(); + let tx = tx.clone(); + thread::spawn(move || { + let nodata = input.configs.nodata; + let mut z: f64; + let mut interior_pit_found = false; + let dx = [1, 1, 1, 0, -1, -1, -1, 0]; + let dy = [-1, 0, 1, 1, 1, 0, -1, -1]; + let mut neighbouring_nodata: bool; + // Create a mapping from the pointer values to cells offsets. + // This may seem wasteful, using only 8 of 129 values in the array, + // but the mapping method is far faster than calculating z.ln() / ln(2.0). + // It's also a good way of allowing for different point styles. + let mut pntr_matches: [i8; 129] = [-2i8; 129]; + if !esri_style { + // This maps Whitebox-style D8 pointer values + // onto the cell offsets in d_x and d_y. + pntr_matches[1] = 0i8; + pntr_matches[2] = 1i8; + pntr_matches[4] = 2i8; + pntr_matches[8] = 3i8; + pntr_matches[16] = 4i8; + pntr_matches[32] = 5i8; + pntr_matches[64] = 6i8; + pntr_matches[128] = 7i8; + } else { + // This maps Esri-style D8 pointer values + // onto the cell offsets in d_x and d_y. + pntr_matches[1] = 1i8; + pntr_matches[2] = 2i8; + pntr_matches[4] = 3i8; + pntr_matches[8] = 4i8; + pntr_matches[16] = 5i8; + pntr_matches[32] = 6i8; + pntr_matches[64] = 7i8; + pntr_matches[128] = 0i8; + } + for row in (0..rows).filter(|r| r % num_procs == tid) { + let mut data: Vec = vec![-1i8; columns as usize]; + for col in 0..columns { + z = input.get_value(row, col); + if z != nodata { + if z > 0f64 { + data[col as usize] = pntr_matches[z as usize]; + } else { + data[col as usize] = -1i8; + // is this no-flow cell interior? + neighbouring_nodata = false; + for i in 0..8 { + if input.get_value(row + dy[i], col + dx[i]) == nodata { + neighbouring_nodata = true; + } + } + if !neighbouring_nodata { + interior_pit_found = true; + } } } - } else { - data[col as usize] = -1i8; } + tx.send((row, data, interior_pit_found)).unwrap(); } - tx.send((row, data, interior_pit_found)).unwrap(); - } - }); - } - - let mut interior_pit_found = false; - for r in 0..rows { - let (row, data, pit) = rx.recv().expect("Error receiving data from thread."); - flow_dir.set_row_data(row, data); //(data.0, data.1); - if pit { - interior_pit_found = true; + }); } - if verbose { - progress = (100.0_f64 * r as f64 / (rows - 1) as f64) as usize; - if progress != old_progress { - println!("Flow directions: {}%", progress); - old_progress = progress; + + for r in 0..rows { + let (row, data, pit) = rx.recv().expect("Error receiving data from thread."); + flow_dir.set_row_data(row, data); + if pit { + interior_pit_found = true; + } + if verbose { + progress = (100.0_f64 * r as f64 / (rows - 1) as f64) as usize; + if progress != old_progress { + println!("Flow directions: {}%", progress); + old_progress = progress; + } } } } + let mut output = Raster::initialize_using_file(&output_file, &input); + output.configs.photometric_interp = PhotometricInterpretation::Continuous; // if the input is a pointer, this may not be the case by default. + output.reinitialize_values(1.0); + drop(input); + // calculate the number of inflowing cells let flow_dir = Arc::new(flow_dir); let mut num_inflowing: Array2D = Array2D::new(rows, columns, -1, -1)?; let (tx, rx) = mpsc::channel(); for tid in 0..num_procs { - let input = input.clone(); + // let input = input.clone(); let flow_dir = flow_dir.clone(); let tx = tx.clone(); thread::spawn(move || { let dx = [1, 1, 1, 0, -1, -1, -1, 0]; let dy = [-1, 0, 1, 1, 1, 0, -1, -1]; let inflowing_vals: [i8; 8] = [4, 5, 6, 7, 0, 1, 2, 3]; - let mut z: f64; + // let mut z: f64; let mut count: i8; for row in (0..rows).filter(|r| r % num_procs == tid) { let mut data: Vec = vec![-1i8; columns as usize]; for col in 0..columns { - z = input[(row, col)]; - if z != nodata { + // z = input.get_value(row, col); + // if z != nodata { + if flow_dir.get_value(row, col) != -2i8 { count = 0i8; for i in 0..8 { - if flow_dir[(row + dy[i], col + dx[i])] == inflowing_vals[i] { + if flow_dir.get_value(row + dy[i], col + dx[i]) == inflowing_vals[i] { count += 1; } } @@ -372,8 +490,6 @@ impl WhiteboxTool for D8FlowAccumulation { }); } - 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; for r in 0..rows { @@ -400,7 +516,6 @@ impl WhiteboxTool for D8FlowAccumulation { let dy = [-1, 0, 1, 1, 1, 0, -1, -1]; let (mut row, mut col): (isize, isize); let (mut row_n, mut col_n): (isize, isize); - // let mut cell: (isize, isize); let mut dir: i8; let mut fa: f64; while !stack.is_empty() { @@ -467,7 +582,8 @@ impl WhiteboxTool for D8FlowAccumulation { if log_transform { for row in 0..rows { for col in 0..columns { - if input[(row, col)] == nodata { + // if input[(row, col)] == nodata { + if flow_dir.get_value(row, col) == -2 { output[(row, col)] = nodata; } else { let dir = flow_dir[(row, col)]; @@ -492,15 +608,15 @@ impl WhiteboxTool for D8FlowAccumulation { } else { for row in 0..rows { for col in 0..columns { - if input[(row, col)] == nodata { + // if input[(row, col)] == nodata { + if flow_dir.get_value(row, col) == -2 { output[(row, col)] = nodata; } else { - let dir = flow_dir[(row, col)]; + let dir = flow_dir.get_value(row, col); if dir >= 0 { - output[(row, col)] = - output[(row, col)] * cell_area / flow_widths[dir as usize]; + output.set_value(row, col, output.get_value(row, col) * cell_area / flow_widths[dir as usize]); } else { - output[(row, col)] = output[(row, col)] * cell_area / flow_widths[3]; + output.set_value(row, col, output.get_value(row, col) * cell_area / flow_widths[3]); } } } diff --git a/src/tools/hydro_analysis/dinf_flow_accum.rs b/src/tools/hydro_analysis/dinf_flow_accum.rs index b87700b59..f62d1ae79 100644 --- a/src/tools/hydro_analysis/dinf_flow_accum.rs +++ b/src/tools/hydro_analysis/dinf_flow_accum.rs @@ -2,7 +2,7 @@ This tool is part of the WhiteboxTools geospatial analysis library. Authors: Dr. John Lindsay Created: 24/06/2017 -Last Modified: 13/02/2020 +Last Modified: 21/02/2020 License: MIT */ @@ -22,11 +22,12 @@ use std::thread; /// This tool is used to generate a flow accumulation grid (i.e. contributing area) using the D-infinity algorithm /// (Tarboton, 1997). This algorithm is an examples of a multiple-flow-direction (MFD) method because the flow entering /// each grid cell is routed to one or two downslope neighbour, i.e. flow divergence is permitted. The user must -/// specify the name of the input digital elevation model (`--dem`). The DEM should have been hydrologically corrected +/// specify the name of the input digital elevation model or D-infinity pointer raster (`--input`). If an input DEM is +/// specified, the DEM should have been hydrologically corrected /// to remove all spurious depressions and flat areas. DEM pre-processing is usually achieved using the -/// `BreachDepressions` or `FillDepressions` tool. +/// `BreachDepressionsLeastCost` or `FillDepressions` tool. /// -/// In addition to the input DEM grid name, the user must specify the output type (`--out_type`). The output +/// In addition to the input DEM/pointer raster name, the user must specify the output type (`--out_type`). The output /// flow-accumulation /// can be 1) specific catchment area (SCA), which is the upslope contributing area divided by the contour length (taken /// as the grid resolution), 2) total catchment area in square-metres, or 3) the number of upslope grid cells. The user @@ -39,7 +40,7 @@ use std::thread; /// flow-accumulation grids must not be used to estimate other secondary terrain indices, such as the wetness index, or /// relative stream power index. /// -/// Grid cells possessing the NoData value in the input DEM raster are assigned the NoData value in the output +/// Grid cells possessing the NoData value in the input DEM/pointer raster are assigned the NoData value in the output /// flow-accumulation image. The output raster is of the float data type and continuous data scale. /// /// # Reference @@ -47,7 +48,7 @@ use std::thread; /// elevation models. Water resources research, 33(2), 309-319. /// /// # See Also -/// `DInfPointer`, `MDInfFlowAccumulation`, `BreachDepressions`, `FillDepressions`,` +/// `DInfPointer`, `MDInfFlowAccumulation`, `BreachDepressionsLeastCost`, `FillDepressions`,` pub struct DInfFlowAccumulation { name: String, description: String, @@ -66,9 +67,9 @@ impl DInfFlowAccumulation { let mut parameters = vec![]; parameters.push(ToolParameter { - name: "Input DEM File".to_owned(), - flags: vec!["-i".to_owned(), "--dem".to_owned()], - description: "Input raster DEM file.".to_owned(), + name: "Input DEM or D-inf Pointer File".to_owned(), + flags: vec!["-i".to_owned(), "--input".to_owned()], + description: "Input raster DEM or D-infinity pointer file.".to_owned(), parameter_type: ParameterType::ExistingFile(ParameterFileType::Raster), default_value: None, optional: false, @@ -125,6 +126,15 @@ impl DInfFlowAccumulation { optional: true, }); + parameters.push(ToolParameter { + name: "Is the input raster a D-inf flow pointer?".to_owned(), + flags: vec!["--pntr".to_owned()], + description: "Is the input raster a D-infinity flow pointer rather than a DEM?".to_owned(), + parameter_type: ParameterType::Boolean, + default_value: None, + optional: true, + }); + let sep: String = path::MAIN_SEPARATOR.to_string(); let p = format!("{}", env::current_dir().unwrap().display()); let e = format!("{}", env::current_exe().unwrap().display()); @@ -136,8 +146,8 @@ impl DInfFlowAccumulation { if e.contains(".exe") { short_exe += ".exe"; } - let usage = format!(">>.*{0} -r={1} -v --wd=\"*path*to*data*\" --dem=DEM.tif -o=output.tif --out_type=sca ->>.*{0} -r={1} -v --wd=\"*path*to*data*\" --dem=DEM.tif -o=output.tif --out_type=sca --threshold=10000 --log --clip", short_exe, name).replace("*", &sep); + let usage = format!(">>.*{0} -r={1} -v --wd=\"*path*to*data*\" --input=DEM.tif -o=output.tif --out_type=sca +>>.*{0} -r={1} -v --wd=\"*path*to*data*\" --input=DEM.tif -o=output.tif --out_type=sca --threshold=10000 --log --clip", short_exe, name).replace("*", &sep); DInfFlowAccumulation { name: name, @@ -189,6 +199,7 @@ impl WhiteboxTool for DInfFlowAccumulation { let mut convergence_threshold = f64::INFINITY; let mut log_transform = false; let mut clip_max = false; + let mut pntr_input = false; if args.len() == 0 { return Err(Error::new( @@ -251,6 +262,10 @@ impl WhiteboxTool for DInfFlowAccumulation { if vec.len() == 1 || !vec[1].to_string().to_lowercase().contains("false") { clip_max = true; } + } else if flag_val == "-pntr" { + if vec.len() == 1 || !vec[1].to_string().to_lowercase().contains("false") { + pntr_input = true; + } } } @@ -289,134 +304,138 @@ impl WhiteboxTool for DInfFlowAccumulation { // calculate the flow directions let mut flow_dir: Array2D = Array2D::new(rows, columns, nodata, nodata)?; - + let mut interior_pit_found = false; let num_procs = num_cpus::get() as isize; - let (tx, rx) = mpsc::channel(); - for tid in 0..num_procs { - let input = input.clone(); - let tx = tx.clone(); - thread::spawn(move || { - let nodata = input.configs.nodata; - let grid_res = (cell_size_x + cell_size_y) / 2.0; - let mut dir: f64; - let mut max_slope: f64; - let mut e0: f64; - let mut af: f64; - let mut ac: f64; - let (mut e1, mut r, mut s1, mut s2, mut s, mut e2): (f64, f64, f64, f64, f64, f64); - - let ac_vals = [0f64, 1f64, 1f64, 2f64, 2f64, 3f64, 3f64, 4f64]; - let af_vals = [1f64, -1f64, 1f64, -1f64, 1f64, -1f64, 1f64, -1f64]; - - let e1_col = [1, 0, 0, -1, -1, 0, 0, 1]; - let e1_row = [0, -1, -1, 0, 0, 1, 1, 0]; - - let e2_col = [1, 1, -1, -1, -1, -1, 1, 1]; - let e2_row = [-1, -1, -1, -1, 1, 1, 1, 1]; - - let atanof1 = 1.0f64.atan(); - - let mut neighbouring_nodata: bool; - let mut interior_pit_found = false; - const HALF_PI: f64 = PI / 2f64; - for row in (0..rows).filter(|r| r % num_procs == tid) { - let mut data: Vec = vec![nodata; columns as usize]; - for col in 0..columns { - e0 = input[(row, col)]; - if e0 != nodata { - dir = 360.0; - max_slope = f64::MIN; - neighbouring_nodata = false; - for i in 0..8 { - ac = ac_vals[i]; - af = af_vals[i]; - e1 = input[(row + e1_row[i], col + e1_col[i])]; - e2 = input[(row + e2_row[i], col + e2_col[i])]; - if e1 != nodata && e2 != nodata { - if e0 > e1 && e0 > e2 { - s1 = (e0 - e1) / grid_res; - s2 = (e1 - e2) / grid_res; - r = if s1 != 0f64 { - (s2 / s1).atan() - } else { - PI / 2.0 - }; - s = (s1 * s1 + s2 * s2).sqrt(); - if s1 < 0.0 && s2 < 0.0 { - s *= -1.0; - } - if s1 < 0.0 && s2 == 0.0 { - s *= -1.0; - } - if s1 == 0.0 && s2 < 0.0 { - s *= -1.0; - } - if r < 0.0 || r > atanof1 { - if r < 0.0 { + + if !pntr_input { + let (tx, rx) = mpsc::channel(); + for tid in 0..num_procs { + let input = input.clone(); + let tx = tx.clone(); + thread::spawn(move || { + let nodata = input.configs.nodata; + let grid_res = (cell_size_x + cell_size_y) / 2.0; + let mut dir: f64; + let mut max_slope: f64; + let mut e0: f64; + let mut af: f64; + let mut ac: f64; + let (mut e1, mut r, mut s1, mut s2, mut s, mut e2): (f64, f64, f64, f64, f64, f64); + + let ac_vals = [0f64, 1f64, 1f64, 2f64, 2f64, 3f64, 3f64, 4f64]; + let af_vals = [1f64, -1f64, 1f64, -1f64, 1f64, -1f64, 1f64, -1f64]; + + let e1_col = [1, 0, 0, -1, -1, 0, 0, 1]; + let e1_row = [0, -1, -1, 0, 0, 1, 1, 0]; + + let e2_col = [1, 1, -1, -1, -1, -1, 1, 1]; + let e2_row = [-1, -1, -1, -1, 1, 1, 1, 1]; + + let atanof1 = 1.0f64.atan(); + + let mut neighbouring_nodata: bool; + let mut interior_pit_found = false; + const HALF_PI: f64 = PI / 2f64; + for row in (0..rows).filter(|r| r % num_procs == tid) { + let mut data: Vec = vec![nodata; columns as usize]; + for col in 0..columns { + e0 = input[(row, col)]; + if e0 != nodata { + dir = 360.0; + max_slope = f64::MIN; + neighbouring_nodata = false; + for i in 0..8 { + ac = ac_vals[i]; + af = af_vals[i]; + e1 = input[(row + e1_row[i], col + e1_col[i])]; + e2 = input[(row + e2_row[i], col + e2_col[i])]; + if e1 != nodata && e2 != nodata { + if e0 > e1 && e0 > e2 { + s1 = (e0 - e1) / grid_res; + s2 = (e1 - e2) / grid_res; + r = if s1 != 0f64 { + (s2 / s1).atan() + } else { + PI / 2.0 + }; + s = (s1 * s1 + s2 * s2).sqrt(); + if s1 < 0.0 && s2 < 0.0 { + s *= -1.0; + } + if s1 < 0.0 && s2 == 0.0 { + s *= -1.0; + } + if s1 == 0.0 && s2 < 0.0 { + s *= -1.0; + } + if r < 0.0 || r > atanof1 { + if r < 0.0 { + r = 0.0; + s = s1; + } else { + r = atanof1; + s = (e0 - e2) / diag_cell_size; + } + } + if s >= max_slope && s != 0.00001 { + max_slope = s; + dir = af * r + ac * HALF_PI; + } + } else if e0 > e1 || e0 > e2 { + if e0 > e1 { r = 0.0; - s = s1; + s = (e0 - e1) / grid_res; } else { r = atanof1; s = (e0 - e2) / diag_cell_size; } + if s >= max_slope && s != 0.00001 { + max_slope = s; + dir = af * r + ac * HALF_PI; + } } - if s >= max_slope && s != 0.00001 { - max_slope = s; - dir = af * r + ac * HALF_PI; - } - } else if e0 > e1 || e0 > e2 { - if e0 > e1 { - r = 0.0; - s = (e0 - e1) / grid_res; - } else { - r = atanof1; - s = (e0 - e2) / diag_cell_size; - } - if s >= max_slope && s != 0.00001 { - max_slope = s; - dir = af * r + ac * HALF_PI; - } + } else { + neighbouring_nodata = true; } - } else { - neighbouring_nodata = true; } - } - if max_slope > 0f64 { - dir = 360.0 - dir.to_degrees() + 90.0; - if dir > 360.0 { - dir = dir - 360.0; + if max_slope > 0f64 { + dir = 360.0 - dir.to_degrees() + 90.0; + if dir > 360.0 { + dir = dir - 360.0; + } + data[col as usize] = dir; + } else { + data[col as usize] = -1f64; + if !neighbouring_nodata { + interior_pit_found = true; + } } - data[col as usize] = dir; } else { data[col as usize] = -1f64; - if !neighbouring_nodata { - interior_pit_found = true; - } } - } else { - data[col as usize] = -1f64; } + tx.send((row, data, interior_pit_found)).unwrap(); } - tx.send((row, data, interior_pit_found)).unwrap(); - } - }); - } - - let mut interior_pit_found = false; - for r in 0..rows { - let (row, data, pit) = rx.recv().expect("Error receiving data from thread."); - flow_dir.set_row_data(row, data); - if pit { - interior_pit_found = true; + }); } - if verbose { - progress = (100.0_f64 * r as f64 / (rows - 1) as f64) as usize; - if progress != old_progress { - println!("Flow directions: {}%", progress); - old_progress = progress; + + for r in 0..rows { + let (row, data, pit) = rx.recv().expect("Error receiving data from thread."); + flow_dir.set_row_data(row, data); + if pit { + interior_pit_found = true; + } + if verbose { + progress = (100.0_f64 * r as f64 / (rows - 1) as f64) as usize; + if progress != old_progress { + println!("Flow directions: {}%", progress); + old_progress = progress; + } } } + } else { // pointer file input + flow_dir = input.get_data_as_array2d(); } // calculate the number of inflowing cells diff --git a/src/tools/hydro_analysis/fd8_flow_accum.rs b/src/tools/hydro_analysis/fd8_flow_accum.rs index 78d664a3b..6818036e0 100644 --- a/src/tools/hydro_analysis/fd8_flow_accum.rs +++ b/src/tools/hydro_analysis/fd8_flow_accum.rs @@ -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 = 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 = 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); diff --git a/src/tools/hydro_analysis/insert_dams.rs b/src/tools/hydro_analysis/insert_dams.rs index fe84f080d..97b76cafd 100644 --- a/src/tools/hydro_analysis/insert_dams.rs +++ b/src/tools/hydro_analysis/insert_dams.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; // = 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 { diff --git a/whitebox_tools.py b/whitebox_tools.py index 5cf9633d9..8af670576 100644 --- a/whitebox_tools.py +++ b/whitebox_tools.py @@ -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.