Skip to content

Commit

Permalink
Added the InsertDams tool
Browse files Browse the repository at this point in the history 
Added the InsertDams tool
  • Loading branch information
@jblindsay
jblindsay committed Feb 19, 2020
1 parent c3a8712 commit b0c3471
Show file tree
Hide file tree
Showing 8 changed files with 533 additions and 148 deletions.
3 changes: 3 additions & 0 deletions readme.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,9 @@ Version 1.1.1 (XX-XX-20XX)
- Added the RasterPerimeter tool to measure the perimeter of raster polygons.
- Added the MDInfFlowAccumulation tool to perform the MD-infinity flow accumulation of Seibert
and McGlynn (2007).
- Added the InsertDams tool, which can be used to insert impoundment features at a set of points
of interest into a DEM. This tool can be used in combination with the ImpoundmentSizeIndex tool
to create artificial reservoirs/depressions.
- Modified the LidarRbfInterpolation tool to improve efficiency.
- Fixed an issue with how floating point attributes were written in Shapefile attribute tables.
- Updated the LidarSegmentation tool, which now used RANSAC to fit planar models to points.
Expand Down
4 changes: 2 additions & 2 deletions src/tools/data_tools/raster_to_vector_points.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,12 +14,12 @@ use std::f64;
use std::io::{Error, ErrorKind};
use std::path;

/// Converts a raster dataset to a vector of the POINT shapetype. The user must specify
/// Converts a raster data set to a vector of the POINT shapetype. The user must specify
/// the name of a raster file (`--input`) and the name of the output vector (`--output`). Points will correspond
/// with grid cell centre points. All grid cells containing non-zero, non-NoData values
/// will be considered a point. The vector's attribute table will contain a field called
/// 'VALUE' that will contain the cell value for each point feature.
///
///
/// # See Also
/// `RasterToVectorPolygons`, `RasterToVectorLines`
pub struct RasterToVectorPoints {
Expand Down
105 changes: 54 additions & 51 deletions src/tools/data_tools/raster_to_vector_polygons.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,12 +18,12 @@ use std::io::{Error, ErrorKind};
use std::path;
const EPSILON: f64 = std::f64::EPSILON;

/// Converts a raster dataset to a vector of the POLYGON geometry type. The user must specify
/// Converts a raster data set to a vector of the POLYGON geometry type. The user must specify
/// the name of a raster file (`--input`) and the name of the output (`--output`) vector. All grid cells containing
/// non-zero, non-NoData values will be considered a point. The vector's attribute table
/// will contain a field called 'VALUE' that will contain the cell value for each point
/// feature.
///
/// non-zero, non-NoData values will be considered part of a polygon feature. The vector's attribute table
/// will contain a field called 'VALUE' that will contain the cell value for each polygon
/// feature, in addition to the standard feature ID (FID) attribute.
///
/// # See Also
/// `RasterToVectorPoints`, `RasterToVectorLines`
pub struct RasterToVectorPolygons {
Expand Down Expand Up @@ -131,42 +131,6 @@ impl WhiteboxTool for RasterToVectorPolygons {
working_directory: &'a str,
verbose: bool,
) -> Result<(), Error> {
/* Diagram 1:
* Cell Numbering
* _____________
* | | |
* | 0 | 1 |
* |_____|_____|
* | | |
* | 2 | 3 |
* |_____|_____|
*
*/

/* Diagram 2:
* Edge Numbering (shared edges between cells)
* _____________
* | | |
* | 3 |
* |__2__|__0__|
* | | |
* | 1 |
* |_____|_____|
*
*/

/* Diagram 3:
* Cell Edge Numbering
*
* ___0___
* | |
* | |
* 3 1
* | |
* |___2___|
*
*/

let mut input_file = String::new();
let mut output_file = String::new();

Expand Down Expand Up @@ -236,7 +200,7 @@ impl WhiteboxTool for RasterToVectorPolygons {
let west = input.configs.west;
let north = input.configs.north;

let get_x_from_column = |col| -> f64 { west + half_res_x + col as f64 * res_x};
let get_x_from_column = |col| -> f64 { west + half_res_x + col as f64 * res_x };
let get_y_from_row = |row| -> f64 { north - half_res_y - row as f64 * res_y };

let mut output = Shapefile::new(&output_file, ShapeType::Polygon)?;
Expand Down Expand Up @@ -303,6 +267,30 @@ impl WhiteboxTool for RasterToVectorPolygons {
drop(input);
drop(visited);

/* Diagram 1:
* Edge Numbering (shared edges between cells)
* _____________
* | | |
* | 3 |
* |__2__|__0__|
* | | |
* | 1 |
* |_____|_____|
*
*/

/* Diagram 2:
* Cell Edge Numbering
*
* ___0___
* | |
* | |
* 3 1
* | |
* |___2___|
*
*/

let mut edges: Array2D<u8> = Array2D::new(rows, columns, 0u8, 0u8)?;
let mut num_edges: Array2D<u8> = Array2D::new(rows, columns, 0u8, 0u8)?;
let mut cell_edges: u8;
Expand All @@ -317,7 +305,8 @@ impl WhiteboxTool for RasterToVectorPolygons {
zn = clumps.get_value(row + dy[n], col + dx[n]);
if z != zn {
cell_edges |= 1u8 << n;
if n < 4 { // Edges are only counted on the non-diagonal cells
if n < 4 {
// Edges are only counted on the non-diagonal cells
num_edges.increment(row, col, 1u8);
}
}
Expand All @@ -335,7 +324,8 @@ impl WhiteboxTool for RasterToVectorPolygons {
}
}

let mut geometries = vec![ShapefileGeometry::new(ShapeType::Polygon); clump_val as usize-1];
let mut geometries =
vec![ShapefileGeometry::new(ShapeType::Polygon); clump_val as usize - 1];

let mut cell: (isize, isize);
let (mut x, mut y): (f64, f64);
Expand All @@ -350,7 +340,7 @@ impl WhiteboxTool for RasterToVectorPolygons {
let edge_offsets_pt2_x = [half_res_x, half_res_x, -half_res_x, -half_res_x];
let edge_offsets_pt2_y = [half_res_y, -half_res_y, -half_res_y, half_res_y];
let (mut p1, mut p2, mut p3): (Point2D, Point2D, Point2D);
let prec = (5f64*EPSILON).tan();
let prec = (5f64 * EPSILON).tan();
for row in 0..rows {
for col in 0..columns {
z = clumps.get_value(row, col);
Expand Down Expand Up @@ -384,7 +374,10 @@ impl WhiteboxTool for RasterToVectorPolygons {
// is there an edge with the diagonal?
if (edges_val >> (edge + 4)) & 1u8 == 0u8 {
// There's no edge with the diagonal. Move to the diagonal cell.
cell = (cell.0 + dy[edge as usize + 4], cell.1 + dx[edge as usize + 4]);
cell = (
cell.0 + dy[edge as usize + 4],
cell.1 + dx[edge as usize + 4],
);
edge = edge_update[edge as usize];
} else {
// there is an edge with the diagonal
Expand Down Expand Up @@ -414,11 +407,18 @@ impl WhiteboxTool for RasterToVectorPolygons {

if points.len() > 1 {
// Remove unnecessary points
for a in (1..points.len()-1).rev() {
p1 = points[a-1];
for a in (1..points.len() - 1).rev() {
p1 = points[a - 1];
p2 = points[a];
p3 = points[a+1];
if ((p2.y-p1.y)*(p3.x-p2.x)-(p3.y-p2.y)*(p2.x-p1.x)).abs() <= (((p2.x-p1.x)*(p3.x-p2.x)+(p2.y-p1.y)*(p3.y-p2.y))).abs() * prec {
p3 = points[a + 1];
if ((p2.y - p1.y) * (p3.x - p2.x)
- (p3.y - p2.y) * (p2.x - p1.x))
.abs()
<= ((p2.x - p1.x) * (p3.x - p2.x)
+ (p2.y - p1.y) * (p3.y - p2.y))
.abs()
* prec
{
points.remove(a);
}
}
Expand Down Expand Up @@ -449,7 +449,10 @@ impl WhiteboxTool for RasterToVectorPolygons {
for fid in 0..geometries.len() {
output.add_record(geometries[fid].clone());
output.attributes.add_record(
vec![FieldData::Int(fid as i32 + 1), FieldData::Real(clump_to_value[fid + 1])],
vec![
FieldData::Int(fid as i32 + 1),
FieldData::Real(clump_to_value[fid + 1]),
],
false,
);

Expand Down
2 changes: 1 addition & 1 deletion src/tools/hydro_analysis/impoundment_index.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ use std::path::Path;
/// Online resource: [Whitebox Blog](https://whiteboxgeospatial.wordpress.com/2015/04/29/modelling-the-spatial-pattern-of-potential-impoundment-size-from-dems/)
///
/// # See Also
/// `StochasticDepressionAnalysis`
/// `InsertDams`, `StochasticDepressionAnalysis`
pub struct ImpoundmentSizeIndex {
name: String,
description: String,
Expand Down
Loading

0 comments on commit b0c3471

Please sign in to comment.