From f56b347c2c4f20c8d6fe896abd1c96d32c817718 Mon Sep 17 00:00:00 2001 From: John Lindsay Date: Mon, 9 Dec 2019 18:31:33 -0500 Subject: [PATCH] v1.1.0 v1.1.0 --- .DS_Store | Bin 10244 -> 10244 bytes Cargo.lock | 2 +- Cargo.toml | 2 +- readme.txt | 5 +- src/tools/gis_analysis/tin_gridding.rs | 140 +++++++++++++++++++------ whitebox_tools.py | 137 ++++++++++++++++++------ 6 files changed, 219 insertions(+), 67 deletions(-) diff --git a/.DS_Store b/.DS_Store index 9b5f7c955d59016f4a2273ae397a0a2d055ea31e..211dc2f9009f169315c031a80015536818ee6dad 100644 GIT binary patch delta 62 xcmZn(XbIThB0Bk|U@N<+sg8n)sm0`TqOy~lgj6>3i+$tV%&uU|ibM7@BLIjt6YT&1 delta 78 zcmZn(XbIThBFcDt@-4vzPO<80OC1Fh3zNwgL}eMLPktvV&p2Z"] 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 d21ecad3a..5cea08e13 100644 --- a/readme.txt +++ b/readme.txt @@ -56,7 +56,7 @@ for more details. * Release Notes: * ****************** -Version 1.0.3 (09-12-2019) +Version 1.1.0 (09-12-2019) - Added the BreachDepressionsLeastCost tool, which performs a modified form of the Lindsay and Dhun (2015) impact minimizing breaching algorithm. This modified algorithm is very efficient and can provide an excellent method for creating depressionless DEMs from large @@ -85,6 +85,9 @@ Version 1.0.3 (09-12-2019) - Updated the ConstructVectorTIN and TINGridding tools to include a maximum triangle edge length to help avoid the creation of spurious long and narrow triangles in convex regions along the data boundaries. +- Added the ImageCorrelationNeighbourhoodAnalysis tool for performing correlation analysis + between two input rasters within roving search windows. The tool can be used to perform + Pearson's r, Spearman's Rho, or Kendall's Tau-b correlations. Version 1.0.2 (01-11-2019) - Added the BurnStreamsAtRoads tool. diff --git a/src/tools/gis_analysis/tin_gridding.rs b/src/tools/gis_analysis/tin_gridding.rs index 974712be7..4b95ab565 100644 --- a/src/tools/gis_analysis/tin_gridding.rs +++ b/src/tools/gis_analysis/tin_gridding.rs @@ -19,10 +19,19 @@ use std::io::{Error, ErrorKind}; use std::path; /// Creates a raster grid based on a triangular irregular network (TIN) fitted to vector points -/// and linear interpolation within each triangular-shaped plane. +/// and linear interpolation within each triangular-shaped plane. The TIN creation algorithm is based on +/// [Delaunay triangulation](https://en.wikipedia.org/wiki/Delaunay_triangulation). /// +/// The user must specify the attribute field containing point values (`--field`). Alternatively, if the input Shapefile +/// contains z-values, the interpolation may be based on these values (`--use_z`). Either an output grid resolution +/// (`--cell_size`) must be specified or alternatively an existing base file (`--base`) can be used to determine the +/// output raster's (`--output`) resolution and spatial extent. Natural neighbour interpolation generally produces a +/// satisfactorily smooth surface within the region of data points but can produce spurious breaks in the surface +/// outside of this region. Thus, it is recommended that the output surface be clipped to the convex hull of the input +/// points (`--clip`). +/// /// # See Also -/// `LidarTINGridding`, `ConstructVectorTIN` +/// `LidarTINGridding`, `ConstructVectorTIN`, `NaturalNeighbourInterpolation` pub struct TINGridding { name: String, description: String, @@ -90,7 +99,16 @@ impl TINGridding { description: "Output raster's grid resolution.".to_owned(), parameter_type: ParameterType::Float, default_value: None, - optional: false, + optional: true, + }); + + parameters.push(ToolParameter{ + name: "Base Raster File (optional)".to_owned(), + flags: vec!["--base".to_owned()], + description: "Optionally specified input base raster file. Not used when a cell size is specified.".to_owned(), + parameter_type: ParameterType::ExistingFile(ParameterFileType::Raster), + default_value: None, + optional: true }); parameters.push(ToolParameter { @@ -175,7 +193,8 @@ impl WhiteboxTool for TINGridding { let mut use_z = false; let mut use_field = false; let mut output_file: String = "".to_string(); - let mut grid_res: f64 = 1.0; + let mut grid_res: f64 = 0.0; + let mut base_file = String::new(); let mut max_triangle_edge_length = f64::INFINITY; // read the arguments @@ -232,6 +251,12 @@ impl WhiteboxTool for TINGridding { }; max_triangle_edge_length *= max_triangle_edge_length; // actually squared distance + } else if flag_val == "-base" { + base_file = if keyval { + vec[1].to_string() + } else { + args[i + 1].to_string() + }; } } @@ -300,30 +325,81 @@ impl WhiteboxTool for TINGridding { } } - let west: f64 = input.header.x_min; - let north: f64 = input.header.y_max; - let rows: isize = (((north - input.header.y_min) / grid_res).ceil()) as isize; - let columns: isize = (((input.header.x_max - west) / grid_res).ceil()) as isize; - let south: f64 = north - rows as f64 * grid_res; - let east = west + columns as f64 * grid_res; - let nodata = -32768.0f64; + // let west: f64 = input.header.x_min; + // let north: f64 = input.header.y_max; + // let rows: isize = (((north - input.header.y_min) / grid_res).ceil()) as isize; + // let columns: isize = (((input.header.x_max - west) / grid_res).ceil()) as isize; + // let south: f64 = north - rows as f64 * grid_res; + // let east = west + columns as f64 * grid_res; + // let nodata = -32768.0f64; + + // let mut configs = RasterConfigs { + // ..Default::default() + // }; + // configs.rows = rows as usize; + // configs.columns = columns as usize; + // configs.north = north; + // configs.south = south; + // configs.east = east; + // configs.west = west; + // configs.resolution_x = grid_res; + // configs.resolution_y = grid_res; + // configs.nodata = nodata; + // configs.data_type = DataType::F32; + // configs.photometric_interp = PhotometricInterpretation::Continuous; + + // let mut output = Raster::initialize_using_config(&output_file, &configs); - let mut configs = RasterConfigs { - ..Default::default() + let nodata = -32768.0f64; + let mut output = if !base_file.trim().is_empty() || grid_res == 0f64 { + if !base_file.contains(&sep) && !base_file.contains("/") { + base_file = format!("{}{}", working_directory, base_file); + } + let mut base = Raster::new(&base_file, "r")?; + base.configs.nodata = nodata; + Raster::initialize_using_file(&output_file, &base) + } else { + if grid_res == 0f64 { + return Err(Error::new( + ErrorKind::InvalidInput, + "The specified grid resolution is incorrect. Either a non-zero grid resolution \nor an input existing base file name must be used.", + )); + } + // base the output raster on the grid_res and the + // extent of the input vector. + let west: f64 = input.header.x_min; + let north: f64 = input.header.y_max; + let rows: isize = (((north - input.header.y_min) / grid_res).ceil()) as isize; + let columns: isize = (((input.header.x_max - west) / grid_res).ceil()) as isize; + let south: f64 = north - rows as f64 * grid_res; + let east = west + columns as f64 * grid_res; + + let mut configs = RasterConfigs { + ..Default::default() + }; + configs.rows = rows as usize; + configs.columns = columns as usize; + configs.north = north; + configs.south = south; + configs.east = east; + configs.west = west; + configs.resolution_x = grid_res; + configs.resolution_y = grid_res; + configs.nodata = nodata; + configs.data_type = DataType::F32; + configs.photometric_interp = PhotometricInterpretation::Continuous; + + Raster::initialize_using_config(&output_file, &configs) }; - configs.rows = rows as usize; - configs.columns = columns as usize; - configs.north = north; - configs.south = south; - configs.east = east; - configs.west = west; - configs.resolution_x = grid_res; - configs.resolution_y = grid_res; - configs.nodata = nodata; - configs.data_type = DataType::F32; - configs.photometric_interp = PhotometricInterpretation::Continuous; - - let mut output = Raster::initialize_using_config(&output_file, &configs); + + let west = output.configs.west; + let north = output.configs.north; + output.configs.nodata = nodata; // in case a base image is used with a different nodata value. + output.configs.palette = "spectrum.pal".to_string(); + output.configs.data_type = DataType::F32; + output.configs.photometric_interp = PhotometricInterpretation::Continuous; + let res_x = output.configs.resolution_x; + let res_y = output.configs.resolution_y; let mut points: Vec = vec![]; let mut z_values: Vec = vec![]; @@ -414,15 +490,15 @@ impl WhiteboxTool for TINGridding { left = points[p1].x.min(points[p2].x.min(points[p3].x)); right = points[p1].x.max(points[p2].x.max(points[p3].x)); - bottom_row = ((north - bottom) / grid_res).ceil() as isize; - top_row = ((north - top) / grid_res).floor() as isize; - left_col = ((left - west) / grid_res).floor() as isize; - right_col = ((right - west) / grid_res).ceil() as isize; + bottom_row = ((north - bottom) / res_y).ceil() as isize; + top_row = ((north - top) / res_y).floor() as isize; + left_col = ((left - west) / res_x).floor() as isize; + right_col = ((right - west) / res_x).ceil() as isize; for row in top_row..=bottom_row { for col in left_col..=right_col { - x = west + (col as f64 + 0.5) * grid_res; - y = north - (row as f64 + 0.5) * grid_res; + x = west + (col as f64 + 0.5) * res_x; + y = north - (row as f64 + 0.5) * res_y; if point_in_poly(&Point2D::new(x, y), &tri_points) { // calculate the z values z = -(norm.x * x + norm.y * y + k) / norm.z; diff --git a/whitebox_tools.py b/whitebox_tools.py index 188f1728b..5994d8629 100644 --- a/whitebox_tools.py +++ b/whitebox_tools.py @@ -6,7 +6,7 @@ # This script is part of the WhiteboxTools geospatial library. # Authors: Dr. John Lindsay # Created: 28/11/2017 -# Last Modified: 25/07/2019 +# Last Modified: 09/12/2019 # License: MIT from __future__ import print_function @@ -379,6 +379,7 @@ def list_tools(self, keywords=[]): + ############## @@ -885,7 +886,7 @@ def clump(self, i, output, diag=True, zero_back=False, callback=None): if zero_back: args.append("--zero_back") return self.run_tool('clump', args, callback) # returns 1 if error - def construct_vector_tin(self, i, output, field=None, use_z=False, callback=None): + def construct_vector_tin(self, i, output, field=None, use_z=False, max_triangle_edge_length=None, callback=None): """Creates a vector triangular irregular network (TIN) for a set of vector points. Keyword arguments: @@ -894,6 +895,7 @@ def construct_vector_tin(self, i, output, field=None, use_z=False, callback=None field -- Input field name in attribute table. use_z -- Use the 'z' dimension of the Shapefile's geometry instead of an attribute field?. output -- Output vector polygon file. + max_triangle_edge_length -- Optional maximum triangle edge length; triangles larger than this size will not be gridded. callback -- Custom function for handling tool text outputs. """ args = [] @@ -901,6 +903,7 @@ def construct_vector_tin(self, i, output, field=None, use_z=False, callback=None if field is not None: args.append("--field='{}'".format(field)) if use_z: args.append("--use_z") args.append("--output='{}'".format(output)) + if max_triangle_edge_length is not None: args.append("--max_triangle_edge_length='{}'".format(max_triangle_edge_length)) return self.run_tool('construct_vector_tin', args, callback) # returns 1 if error def create_hexagonal_vector_grid(self, i, output, width, orientation="horizontal", callback=None): @@ -1183,6 +1186,30 @@ def minimum_convex_hull(self, i, output, features=True, callback=None): if features: args.append("--features") return self.run_tool('minimum_convex_hull', args, callback) # returns 1 if error + def natural_neighbour_interpolation(self, i, output, field=None, use_z=False, cell_size=None, base=None, clip=True, callback=None): + """Creates a raster grid based on Sibson's natural neighbour method. + + Keyword arguments: + + i -- Input vector points file. + field -- Input field name in attribute table. + use_z -- Use the 'z' dimension of the Shapefile's geometry instead of an attribute field?. + output -- Output raster file. + cell_size -- Optionally specified cell size of output raster. Not used when base raster is specified. + base -- Optionally specified input base raster file. Not used when a cell size is specified. + clip -- Clip the data to the convex hull of the points?. + callback -- Custom function for handling tool text outputs. + """ + args = [] + args.append("--input='{}'".format(i)) + if field is not None: args.append("--field='{}'".format(field)) + if use_z: args.append("--use_z") + args.append("--output='{}'".format(output)) + if cell_size is not None: args.append("--cell_size='{}'".format(cell_size)) + if base is not None: args.append("--base='{}'".format(base)) + if clip: args.append("--clip") + return self.run_tool('natural_neighbour_interpolation', args, callback) # returns 1 if error + def nearest_neighbour_gridding(self, i, field, output, use_z=False, cell_size=None, base=None, max_dist=None, callback=None): """Creates a raster grid based on a set of vector points and assigns grid values using the nearest neighbour. @@ -1365,7 +1392,7 @@ def smooth_vectors(self, i, output, filter=3, callback=None): args.append("--filter={}".format(filter)) return self.run_tool('smooth_vectors', args, callback) # returns 1 if error - def tin_gridding(self, i, output, resolution, field=None, use_z=False, callback=None): + def tin_gridding(self, i, output, field=None, use_z=False, resolution=None, base=None, max_triangle_edge_length=None, callback=None): """Creates a raster grid based on a triangular irregular network (TIN) fitted to vector points. Keyword arguments: @@ -1375,6 +1402,8 @@ def tin_gridding(self, i, output, resolution, field=None, use_z=False, callback= use_z -- Use the 'z' dimension of the Shapefile's geometry instead of an attribute field?. output -- Output raster file. resolution -- Output raster's grid resolution. + base -- Optionally specified input base raster file. Not used when a cell size is specified. + max_triangle_edge_length -- Optional maximum triangle edge length; triangles larger than this size will not be gridded. callback -- Custom function for handling tool text outputs. """ args = [] @@ -1382,7 +1411,9 @@ def tin_gridding(self, i, output, resolution, field=None, use_z=False, callback= if field is not None: args.append("--field='{}'".format(field)) if use_z: args.append("--use_z") args.append("--output='{}'".format(output)) - args.append("--resolution='{}'".format(resolution)) + if resolution is not None: args.append("--resolution='{}'".format(resolution)) + if base is not None: args.append("--base='{}'".format(base)) + if max_triangle_edge_length is not None: args.append("--max_triangle_edge_length='{}'".format(max_triangle_edge_length)) return self.run_tool('tin_gridding', args, callback) # returns 1 if error def vector_hex_binning(self, i, output, width, orientation="horizontal", callback=None): @@ -3203,6 +3234,30 @@ def breach_depressions(self, dem, output, max_depth=None, max_length=None, flat_ if fill_pits: args.append("--fill_pits") return self.run_tool('breach_depressions', args, callback) # returns 1 if error + def breach_depressions_least_cost(self, dem, output, radius, max_cost=None, min_dist=True, flat_increment=None, fill=True, callback=None): + """Breaches the depressions in a DEM using a least-cost pathway method. + + Keyword arguments: + + dem -- Input raster DEM file. + output -- Output raster file. + radius -- . + max_cost -- Optional maximum breach cost (default is Inf). + min_dist -- Optional flag indicating whether to minimize breach distances. + flat_increment -- Optional elevation increment applied to flat areas. + fill -- Optional flag indicating whether to fill any remaining unbreached depressions. + callback -- Custom function for handling tool text outputs. + """ + args = [] + args.append("--dem='{}'".format(dem)) + args.append("--output='{}'".format(output)) + args.append("--radius='{}'".format(radius)) + if max_cost is not None: args.append("--max_cost='{}'".format(max_cost)) + if min_dist: args.append("--min_dist") + if flat_increment is not None: args.append("--flat_increment='{}'".format(flat_increment)) + if fill: args.append("--fill") + return self.run_tool('breach_depressions_least_cost', args, callback) # returns 1 if error + def breach_single_cell_pits(self, dem, output, callback=None): """Removes single-cell pits from an input DEM by breaching. @@ -3487,7 +3542,7 @@ def fill_burn(self, dem, streams, output, callback=None): args.append("--output='{}'".format(output)) return self.run_tool('fill_burn', args, callback) # returns 1 if error - def fill_depressions(self, dem, output, fix_flats=True, flat_increment=None, callback=None): + def fill_depressions(self, dem, output, fix_flats=True, flat_increment=None, max_depth=None, callback=None): """Fills all of the depressions in a DEM. Depression breaching should be preferred in most cases. Keyword arguments: @@ -3496,6 +3551,7 @@ def fill_depressions(self, dem, output, fix_flats=True, flat_increment=None, cal output -- Output raster file. fix_flats -- Optional flag indicating whether flat areas should have a small gradient applied. flat_increment -- Optional elevation increment applied to flat areas. + max_depth -- Optional maximum depression depth to fill. callback -- Custom function for handling tool text outputs. """ args = [] @@ -3503,8 +3559,27 @@ def fill_depressions(self, dem, output, fix_flats=True, flat_increment=None, cal args.append("--output='{}'".format(output)) if fix_flats: args.append("--fix_flats") if flat_increment is not None: args.append("--flat_increment='{}'".format(flat_increment)) + if max_depth is not None: args.append("--max_depth='{}'".format(max_depth)) return self.run_tool('fill_depressions', args, callback) # returns 1 if error + def fill_depressions_wang_and_lui(self, dem, output, fix_flats=True, flat_increment=None, callback=None): + """Fills all of the depressions in a DEM. Depression breaching should be preferred in most cases. + + Keyword arguments: + + dem -- Input raster DEM file. + output -- Output raster file. + fix_flats -- Optional flag indicating whether flat areas should have a small gradient applied. + flat_increment -- Optional elevation increment applied to flat areas. + callback -- Custom function for handling tool text outputs. + """ + args = [] + args.append("--dem='{}'".format(dem)) + args.append("--output='{}'".format(output)) + if fix_flats: args.append("--fix_flats") + if flat_increment is not None: args.append("--flat_increment='{}'".format(flat_increment)) + return self.run_tool('fill_depressions_wang_and_lui', args, callback) # returns 1 if error + def fill_single_cell_pits(self, dem, output, callback=None): """Raises pit cells to the elevation of their lowest neighbour. @@ -3691,30 +3766,6 @@ def jenson_snap_pour_points(self, pour_pts, streams, output, snap_dist, callback args.append("--snap_dist='{}'".format(snap_dist)) return self.run_tool('jenson_snap_pour_points', args, callback) # returns 1 if error - def least_cost_breach_depressions(self, dem, output, radius, max_cost=None, min_dist=True, flat_increment=None, fill=True, callback=None): - """Breaches the depressions in a DEM using a least-cost pathway method. - - Keyword arguments: - - dem -- Input raster DEM file. - output -- Output raster file. - radius -- . - max_cost -- Optional maximum breach cost (default is Inf). - min_dist -- Optional flag indicating whether to minimize breach distances. - flat_increment -- Optional elevation increment applied to flat areas. - fill -- Optional flag indicating whether to fill any remaining unbreached depressions. - callback -- Custom function for handling tool text outputs. - """ - args = [] - args.append("--dem='{}'".format(dem)) - args.append("--output='{}'".format(output)) - args.append("--radius='{}'".format(radius)) - if max_cost is not None: args.append("--max_cost='{}'".format(max_cost)) - if min_dist: args.append("--min_dist") - if flat_increment is not None: args.append("--flat_increment='{}'".format(flat_increment)) - if fill: args.append("--fill") - return self.run_tool('least_cost_breach_depressions', args, callback) # returns 1 if error - def longest_flowpath(self, dem, basins, output, callback=None): """Delineates the longest flowpaths for a group of subbasins or watersheds. @@ -3795,18 +3846,18 @@ def rho8_pointer(self, dem, output, esri_pntr=False, callback=None): if esri_pntr: args.append("--esri_pntr") return self.run_tool('rho8_pointer', args, callback) # returns 1 if error - def sink(self, dem, output, zero_background=False, callback=None): + def sink(self, i, output, zero_background=False, callback=None): """Identifies the depressions in a DEM, giving each feature a unique identifier. Keyword arguments: - dem -- Input raster DEM file. + i -- Input raster DEM file. output -- Output raster file. zero_background -- Flag indicating whether a background value of zero should be used. callback -- Custom function for handling tool text outputs. """ args = [] - args.append("--dem='{}'".format(dem)) + args.append("--input='{}'".format(i)) args.append("--output='{}'".format(output)) if zero_background: args.append("--zero_background") return self.run_tool('sink', args, callback) # returns 1 if error @@ -6493,6 +6544,28 @@ def image_correlation(self, inputs, output=None, callback=None): if output is not None: args.append("--output='{}'".format(output)) return self.run_tool('image_correlation', args, callback) # returns 1 if error + def image_correlation_neighbourhood_analysis(self, input1, input2, output1, output2, filter=11, stat="pearson", callback=None): + """Performs image correlation on two input images neighbourhood search windows. + + Keyword arguments: + + input1 -- Input raster file. + input2 -- Input raster file. + output1 -- Output correlation (r-value or rho) raster file. + output2 -- Output significance (p-value) raster file. + filter -- Size of the filter kernel. + stat -- Correlation type; one of 'pearson' (default) and 'spearman'. + callback -- Custom function for handling tool text outputs. + """ + args = [] + args.append("--input1='{}'".format(input1)) + args.append("--input2='{}'".format(input2)) + args.append("--output1='{}'".format(output1)) + args.append("--output2='{}'".format(output2)) + args.append("--filter={}".format(filter)) + args.append("--stat={}".format(stat)) + return self.run_tool('image_correlation_neighbourhood_analysis', args, callback) # returns 1 if error + def image_regression(self, input1, input2, output, out_residuals=None, standardize=False, scattergram=False, num_samples=1000, callback=None): """Performs image regression analysis on two input images.