From b569c7575f3d9225df6f3dee07d1c142128252a2 Mon Sep 17 00:00:00 2001 From: John Lindsay Date: Fri, 4 Sep 2020 09:43:09 -0400 Subject: [PATCH] Fixes for #107 and #109 --- .DS_Store | Bin 10244 -> 10244 bytes Cargo.lock | 2 +- Cargo.toml | 2 +- readme.txt | 7 +- src/tools/image_analysis/mosaic.rs | 38 ++-- src/tools/lidar_analysis/lidar_dsm.rs | 52 +++++ .../lidar_analysis/lidar_idw_interpolation.rs | 34 +++ src/tools/lidar_analysis/lidar_nn_gridding.rs | 50 ++++- .../lidar_analysis/lidar_segmentation.rs | 2 + .../lidar_analysis/lidar_tin_gridding.rs | 36 ++++ src/tools/terrain_analysis/horizon_angle.rs | 8 + src/tools/terrain_analysis/map_otos.rs | 11 +- .../terrain_analysis/max_diff_from_mean.rs | 2 +- .../terrain_analysis/max_elev_deviation.rs | 2 +- .../remove_off_terrain_objects.rs | 2 +- .../terrain_analysis/time_in_daylight.rs | 38 +++- whitebox_tools.py | 193 +++++++++++++----- 17 files changed, 398 insertions(+), 81 deletions(-) diff --git a/.DS_Store b/.DS_Store index e562abf11514d5c89308619a61ba7de5c2ece00f..a293feee3a97c5f54d40f8775bf179e8424ae1e4 100644 GIT binary patch delta 42 ycmZn(XbISmCOkP+w3Xe+R7b(Y)OhkoVcE%iVrrXBL@se}W>@&bvYAPgkr@CfEey*5 delta 75 zcmZn(XbISmCd_zZa=LJXhIn"] 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 b1c253e29..9c68a986b 100644 --- a/readme.txt +++ b/readme.txt @@ -56,7 +56,7 @@ for more details. * Release Notes: * ****************** -Version 1.4.0 (03-09-2020) +Version 1.4.0 (04-09-2020) - Added the TimeInDaylight model tool for modelling the proportion of daytime that a location is not in shadow. - Added the MapOffTerrainObjects tool. - Added the FilterRasterFeaturesByArea tool. @@ -68,6 +68,11 @@ Version 1.4.0 (03-09-2020) - The Resample tool has been modified so that it does not require a 'destination' raster. Instead, it will create a new output raster either based on a user-specified target cell resolution or an optional base raster, much like the vector-to-raster conversion tools. +- Tools that input a z_factor conversion no longer override user input with geographic coordinates + (see issue #113). +- The StreamLinkIdentifier tool now outputs a 32-bit integer format, increasing the maximum allowable + number of streams (see issue #110). +- Fixed a bug with cubic-convolution and bilinear resampling in the Mosaic tool (see issue #109). Version 1.3.1 (23-07-2020) - Added the HypsometricallyTintedHillshade tool to create hypsometric tinted hillshades. diff --git a/src/tools/image_analysis/mosaic.rs b/src/tools/image_analysis/mosaic.rs index d83c6fef8..62b2bc8fd 100644 --- a/src/tools/image_analysis/mosaic.rs +++ b/src/tools/image_analysis/mosaic.rs @@ -2,7 +2,7 @@ This tool is part of the WhiteboxTools geospatial analysis library. Authors: Dr. John Lindsay Created: 02/01/2018 -Last Modified: 02/01/2018 +Last Modified: 03/09/2020 License: MIT */ @@ -535,22 +535,20 @@ impl WhiteboxTool for Mosaic { let ret = tree .locate_all_at_point(&[x[col as usize], y[row as usize]]) .collect::>(); + + let mut flag = true; for a in 0..ret.len() { i = ret[a].data; - // row_src = inputs[i].get_row_from_y(y[row as usize]); - // col_src = inputs[i].get_column_from_x(x[col as usize]); + if !flag { + break; + } + row_src = (inputs[i].configs.north - y[row as usize]) / inputs[i].configs.resolution_y; col_src = (x[col as usize] - inputs[i].configs.west) / inputs[i].configs.resolution_x; - origin_row = row_src.floor() as isize; origin_col = col_src.floor() as isize; - - // z = inputs[i].get_value(origin_row, origin_col); - - // if origin_row < 0 || origin_col < 0 { break; } - // if origin_row > inputs[i].configs.rows as isize || origin_col > inputs[i].configs.columns as isize { break; } sum_dist = 0f64; for n in 0..num_neighbours { row_n = origin_row + shift_y[n]; @@ -566,8 +564,7 @@ impl WhiteboxTool for Mosaic { neighbour[n][1] = 0f64; } else { data[col as usize] = neighbour[n][0]; - // flag = false; - break; + flag = false; } } @@ -577,9 +574,7 @@ impl WhiteboxTool for Mosaic { z += (neighbour[n][0] * neighbour[n][1]) / sum_dist; } data[col as usize] = z; - - // flag = false; - break; + flag = false; } } } @@ -633,19 +628,18 @@ impl WhiteboxTool for Mosaic { let ret = tree .locate_all_at_point(&[x[col as usize], y[row as usize]]) .collect::>(); + let mut flag = true; for a in 0..ret.len() { i = ret[a].data; - // if !flag { - // break; - // } + if !flag { + break; + } row_src = (inputs[i].configs.north - y[row as usize]) / inputs[i].configs.resolution_y; col_src = (x[col as usize] - inputs[i].configs.west) / inputs[i].configs.resolution_x; origin_row = row_src.floor() as isize; origin_col = col_src.floor() as isize; - // if origin_row < 0 || origin_col < 0 { break; } - // if origin_row > inputs[i].configs.rows as isize || origin_col > inputs[i].configs.columns as isize { break; } sum_dist = 0f64; for n in 0..num_neighbours { row_n = origin_row + shift_y[n]; @@ -661,8 +655,7 @@ impl WhiteboxTool for Mosaic { neighbour[n][1] = 0f64; } else { data[col as usize] = neighbour[n][0]; - break; - // flag = false; + flag = false; } } @@ -672,8 +665,7 @@ impl WhiteboxTool for Mosaic { z += (neighbour[n][0] * neighbour[n][1]) / sum_dist; } data[col as usize] = z; - break; - // flag = false; + flag = false; } } } diff --git a/src/tools/lidar_analysis/lidar_dsm.rs b/src/tools/lidar_analysis/lidar_dsm.rs index 7e6ab3a92..dab0d0f9e 100644 --- a/src/tools/lidar_analysis/lidar_dsm.rs +++ b/src/tools/lidar_analysis/lidar_dsm.rs @@ -21,6 +21,58 @@ use std::sync::Arc; use std::{env, f64, fs, path, thread}; // use rayon::prelude::*; +/// This tool creates a digital surface model (DSM) from a LiDAR point cloud. A DSM reflects the elevation of the tops +/// of all off-terrain objects (i.e. non-ground features) contained within the data set. For example, a DSM will model +/// the canopy top as well as building roofs. This is in stark contrast to a bare-earth digital elevation model (DEM), +/// which models the ground surface without off-terrain objects present. Bare-earth DEMs can be derived from LiDAR data +/// by interpolating last-return points using one of the other LiDAR interpolators (e.g. `LidarTINGridding`). The algorithm +/// used for interpolation in this tool is based on gridding a triangulation (TIN) fit to top-level points in the +/// input LiDAR point cloud. All points in the input LiDAR data set that are below other neighbouring points, within +/// a specified search radius (`--radius`), and that have a large inter-point slope, are filtered out. Thus, this tool +/// will remove the ground surface beneath as well as any intermediate points within a forest canopy, leaving only the +/// canopy top surface to be interpolated. Similarly, building wall points and any ground points beneath roof overhangs +/// will also be remove prior to interpolation. Note that because the ground points beneath overhead wires and utility +/// lines are filtered out by this operation, these features tend to be appear as 'walls' in the output DSM. If these +/// points are classified in the input LiDAR file, you may wish to filter them out before using this tool (`FilterLidarClasses`). +/// +/// The following images show the differences between creating a DSM using the `LidarDigitalSurfaceModel` +/// and by interpolating first-return points only using the `LidarTINGridding` tool respectively. Note, the images show +/// `TimeInDaylight`, which is a more effective way of hillshading DSMs than the traditional `Hillshade` method. Compare +/// how the DSM created `LidarDigitalSurfaceModel` tool (above) has far less variability in areas of tree-cover, more +/// effectively capturing the canopy top. As well, notice how building rooftops are more extensive and straighter in +/// the `LidarDigitalSurfaceModel` DSM image. This is because this method eliminates ground returns beneath roof overhangs +/// before the triangulation operation. +/// +/// ![](../../doc_img/DSM.png) +/// +/// ![](../../doc_img/FirstReturnTIN.png) +/// +/// The user must specify the grid resolution of the output raster (`--resolution`), and optionally, the name of the +/// input LiDAR file (`--input`) and output raster (`--output`). Note that if an input LiDAR file (`--input`) is not +/// specified by the user, the tool will search for all valid LiDAR (*.las, *.zlidar) contained within the current +/// working directory. This feature can be very useful when you need to interpolate a DSM for a large number of LiDAR +/// files. Not only does this batch processing mode enable the tool to run in a more optimized parallel manner, but it +/// will also allow the tool to include a small buffer of points extending into adjacent tiles when interpolating an +/// individual file. This can significantly reduce edge-effects when the output tiles are later mosaicked together. +/// When run in this batch mode, the output file (`--output`) also need not be specified; the tool will instead create +/// an output file with the same name as each input LiDAR file, but with the .tif extension. This can provide a very +/// efficient means for processing extremely large LiDAR data sets. +/// +/// Users may also exclude points from the interpolation if they fall below or above the minimum (`--minz`) or +/// maximum (`--maxz`) thresholds respectively. This can be a useful means of excluding anomalously high or low +/// points. Note that points that are classified as low points (LAS class 7) or high noise (LAS class 18) are +/// automatically excluded from the interpolation operation. +/// +/// Triangulation will generally completely fill the convex hull containing the input point data. This can sometimes +/// result in very long and narrow triangles at the edges of the data or connecting vertices on either side of void +/// areas. In LiDAR data, these void areas are often associated with larger waterbodies, and triangulation can result +/// in very unnatural interpolated patterns within these areas. To avoid this problem, the user may specify a the +/// maximum allowable triangle edge length (`max_triangle_edge_length`) and all grid cells within triangular facets +/// with edges larger than this threshold are simply assigned the NoData values in the output DSM. These NoData areas +/// can later be better dealt with using the `FillMissingData` tool after interpolation. +/// +/// # See Also +/// `LidarTINGridding`, `FilterLidarClasses`, `FillMissingData`, `TimeInDaylight` pub struct LidarDigitalSurfaceModel { name: String, description: String, diff --git a/src/tools/lidar_analysis/lidar_idw_interpolation.rs b/src/tools/lidar_analysis/lidar_idw_interpolation.rs index f860a03d9..617e56d8c 100644 --- a/src/tools/lidar_analysis/lidar_idw_interpolation.rs +++ b/src/tools/lidar_analysis/lidar_idw_interpolation.rs @@ -25,6 +25,40 @@ use std::sync::mpsc; use std::sync::{Arc, Mutex}; use std::thread; +/// This tool interpolates LiDAR files using [inverse-distance weighting](https://en.wikipedia.org/wiki/Inverse_distance_weighting) +/// (IDW) scheme. The user must specify the value of the IDW weight parameter (`--weight`). The output grid can be +/// based on any of the stored LiDAR point parameters (`--parameter`), including elevation +/// (in which case the output grid is a digital elevation model, DEM), intensity, class, return number, number of +/// returns, scan angle, RGB (colour) values, and user data values. Similarly, the user may specify which point +/// return values (`--returns`) to include in the interpolation, including all points, last returns (including single +/// return points), and first returns (including single return points). +/// +/// The user must specify the grid resolution of the output raster (`--resolution`), and optionally, the name of the +/// input LiDAR file (`--input`) and output raster (`--output`). Note that if an input LiDAR file (`--input`) is not +/// specified by the user, the tool will search for all valid LiDAR (*.las, *.zlidar) contained within the current +/// working directory. This feature can be very useful when you need to interpolate a DEM for a large number of LiDAR +/// files. Not only does this batch processing mode enable the tool to run in a more optimized parallel manner, but it +/// will also allow the tool to include a small buffer of points extending into adjacent tiles when interpolating an +/// individual file. This can significantly reduce edge-effects when the output tiles are later mosaicked together. +/// When run in this batch mode, the output file (`--output`) also need not be specified; the tool will instead create +/// an output file with the same name as each input LiDAR file, but with the .tif extension. This can provide a very +/// efficient means for processing extremely large LiDAR data sets. +/// +/// Users may excluded points from the interpolation based on point classification values, which follow the LAS +/// classification scheme. Excluded classes are specified using the `--exclude_cls` parameter. For example, +/// to exclude all vegetation and building classified points from the interpolation, use --exclude_cls='3,4,5,6'. +/// Users may also exclude points from the interpolation if they fall below or above the minimum (`--minz`) or +/// maximum (`--maxz`) thresholds respectively. This can be a useful means of excluding anomalously high or low +/// points. Note that points that are classified as low points (LAS class 7) or high noise (LAS class 18) are +/// automatically excluded from the interpolation operation. +/// +/// The tool will search for the nearest input LiDAR point to each grid cell centre, up to a maximum search distance +/// (`--radius`). If a grid cell does not have a LiDAR point within this search distance, it will be assigned the +/// NoData value in the output raster. In LiDAR data, these void areas are often associated with larger waterbodies. +/// These NoData areas can later be better dealt with using the `FillMissingData` tool after interpolation. +/// +/// # See Also +/// `LidarTINGridding`, `LidarNearestNeighbourGridding`, `LidarTINGridding` pub struct LidarIdwInterpolation { name: String, description: String, diff --git a/src/tools/lidar_analysis/lidar_nn_gridding.rs b/src/tools/lidar_analysis/lidar_nn_gridding.rs index 4a89cf4dc..d9da99250 100644 --- a/src/tools/lidar_analysis/lidar_nn_gridding.rs +++ b/src/tools/lidar_analysis/lidar_nn_gridding.rs @@ -8,7 +8,7 @@ License: MIT NOTES: 1. This tool is designed to work either by specifying a single input and output file or a working directory containing multiple input LAS files. -2. Need to add the ability to exclude points based on max scan angle divation. +2. Need to add the ability to exclude points based on max scan angle devation. */ use crate::lidar::*; @@ -25,6 +25,52 @@ use std::sync::mpsc; use std::sync::{Arc, Mutex}; use std::thread; +/// This tool grids LiDAR files using nearest-neighbour (NN) scheme, that is, each grid cell in the output image will +/// be assigned the parameter value of the point nearest the grid cell centre. This method should not be confused +/// for the similarly named [natural-neighbour interpolation](https://en.wikipedia.org/wiki/Natural_neighbor_interpolation) +/// (a.k.a Sibson's method). Nearest neighbour gridding is generally regarded as a poor way of interpolating surfaces +/// from low-density point sets and results in the creation of a [Voronoi +/// diagram](https://en.wikipedia.org/wiki/Voronoi_diagram). However, this method has several advantages when +/// applied to LiDAR data. NN gridding is one of the fastest methods for generating raster surfaces from large +/// LiDAR data sets. NN gridding is one of the few interpolation methods, along with triangulation, that will preserve +/// vertical breaks-in-slope, such as occur at the edges of building. This characteristic can be important when using +/// some post-processing methods, such as the `RemoveOffTerrainObjects` tool. Furthermore, because most LiDAR data +/// sets have remarkably high point densities compared with other types of geographic data, this approach does often +/// produce a satisfactory result; this is particularly true when the point density is high enough that there are multiple +/// points in the majority of grid cells. +/// +/// The output grid can be based on any of the stored LiDAR point parameters (`--parameter`), including elevation +/// (in which case the output grid is a digital elevation model, DEM), intensity, class, return number, number of +/// returns, scan angle, RGB (colour) values, and user data values. Similarly, the user may specify which point +/// return values (`--returns`) to include in the interpolation, including all points, last returns (including single +/// return points), and first returns (including single return points). +/// +/// The user must specify the grid resolution of the output raster (`--resolution`), and optionally, the name of the +/// input LiDAR file (`--input`) and output raster (`--output`). Note that if an input LiDAR file (`--input`) is not +/// specified by the user, the tool will search for all valid LiDAR (*.las, *.zlidar) contained within the current +/// working directory. This feature can be very useful when you need to interpolate a DEM for a large number of LiDAR +/// files. Not only does this batch processing mode enable the tool to run in a more optimized parallel manner, but it +/// will also allow the tool to include a small buffer of points extending into adjacent tiles when interpolating an +/// individual file. This can significantly reduce edge-effects when the output tiles are later mosaicked together. +/// When run in this batch mode, the output file (`--output`) also need not be specified; the tool will instead create +/// an output file with the same name as each input LiDAR file, but with the .tif extension. This can provide a very +/// efficient means for processing extremely large LiDAR data sets. +/// +/// Users may excluded points from the interpolation based on point classification values, which follow the LAS +/// classification scheme. Excluded classes are specified using the `--exclude_cls` parameter. For example, +/// to exclude all vegetation and building classified points from the interpolation, use --exclude_cls='3,4,5,6'. +/// Users may also exclude points from the interpolation if they fall below or above the minimum (`--minz`) or +/// maximum (`--maxz`) thresholds respectively. This can be a useful means of excluding anomalously high or low +/// points. Note that points that are classified as low points (LAS class 7) or high noise (LAS class 18) are +/// automatically excluded from the interpolation operation. +/// +/// The tool will search for the nearest input LiDAR point to each grid cell centre, up to a maximum search distance +/// (`--radius`). If a grid cell does not have a LiDAR point within this search distance, it will be assigned the +/// NoData value in the output raster. In LiDAR data, these void areas are often associated with larger waterbodies. +/// These NoData areas can later be better dealt with using the `FillMissingData` tool after interpolation. +/// +/// # See Also +/// `LidarTINGridding`, `LidarIdwInterpolation`, `LidarTINGridding`, `RemoveOffTerrainObjects`, `FillMissingData` pub struct LidarNearestNeighbourGridding { name: String, description: String, @@ -38,7 +84,7 @@ impl LidarNearestNeighbourGridding { // public constructor let name = "LidarNearestNeighbourGridding".to_string(); let toolbox = "LiDAR Tools".to_string(); - let description = "Grids LAS files using nearest-neighbour scheme. When the input/output parameters are not specified, the tool grids all LAS files contained within the working directory.".to_string(); + let description = "Grids LiDAR files using nearest-neighbour scheme. When the input/output parameters are not specified, the tool grids all LAS files contained within the working directory.".to_string(); let mut parameters = vec![]; parameters.push(ToolParameter { diff --git a/src/tools/lidar_analysis/lidar_segmentation.rs b/src/tools/lidar_analysis/lidar_segmentation.rs index 150cdaf47..e0e7f5986 100644 --- a/src/tools/lidar_analysis/lidar_segmentation.rs +++ b/src/tools/lidar_analysis/lidar_segmentation.rs @@ -65,6 +65,8 @@ use std::thread; /// Planes that have slopes greater than this threshold are rejected by the algorithm. This has the side-effect /// of removing building walls however. /// +/// ![](../../doc_img/LidarSegmentation.png) +/// /// # References /// Fischler MA and Bolles RC. 1981. Random sample consensus: a paradigm for model fitting with applications /// to image analysis and automated cartography. Commun. ACM, 24(6):381–395. diff --git a/src/tools/lidar_analysis/lidar_tin_gridding.rs b/src/tools/lidar_analysis/lidar_tin_gridding.rs index 8a04b92b7..452719936 100644 --- a/src/tools/lidar_analysis/lidar_tin_gridding.rs +++ b/src/tools/lidar_analysis/lidar_tin_gridding.rs @@ -21,6 +21,42 @@ use std::sync::mpsc; use std::sync::Arc; use std::{env, f64, fs, path, thread}; +/// This tool creates a raster grid based on a Delaunay triangular irregular network (TIN) fitted to LiDAR points. +/// The output grid can be based on any of the stored LiDAR point parameters (`--parameter`), including elevation +/// (in which case the output grid is a digital elevation model, DEM), intensity, class, return number, number of +/// returns, scan angle, RGB (colour) values, and user data values. Similarly, the user may specify which point +/// return values (`--returns`) to include in the interpolation, including all points, last returns (including single +/// return points), and first returns (including single return points). +/// +/// The user must specify the grid resolution of the output raster (`--resolution`), and optionally, the name of the +/// input LiDAR file (`--input`) and output raster (`--output`). Note that if an input LiDAR file (`--input`) is not +/// specified by the user, the tool will search for all valid LiDAR (*.las, *.zlidar) contained within the current +/// working directory. This feature can be very useful when you need to interpolate a DEM for a large number of LiDAR +/// files. Not only does this batch processing mode enable the tool to run in a more optimized parallel manner, but it +/// will also allow the tool to include a small buffer of points extending into adjacent tiles when interpolating an +/// individual file. This can significantly reduce edge-effects when the output tiles are later mosaicked together. +/// When run in this batch mode, the output file (`--output`) also need not be specified; the tool will instead create +/// an output file with the same name as each input LiDAR file, but with the .tif extension. This can provide a very +/// efficient means for processing extremely large LiDAR data sets. +/// +/// Users may excluded points from the interpolation based on point classification values, which follow the LAS +/// classification scheme. Excluded classes are specified using the `--exclude_cls` parameter. For example, +/// to exclude all vegetation and building classified points from the interpolation, use --exclude_cls='3,4,5,6'. +/// Users may also exclude points from the interpolation if they fall below or above the minimum (`--minz`) or +/// maximum (`--maxz`) thresholds respectively. This can be a useful means of excluding anomalously high or low +/// points. Note that points that are classified as low points (LAS class 7) or high noise (LAS class 18) are +/// automatically excluded from the interpolation operation. +/// +/// Triangulation will generally completely fill the convex hull containing the input point data. This can sometimes +/// result in very long and narrow triangles at the edges of the data or connecting vertices on either side of void +/// areas. In LiDAR data, these void areas are often associated with larger waterbodies, and triangulation can result +/// in very unnatural interpolated patterns within these areas. To avoid this problem, the user may specify a the +/// maximum allowable triangle edge length (`max_triangle_edge_length`) and all grid cells within triangular facets +/// with edges larger than this threshold are simply assigned the NoData values in the output DSM. These NoData areas +/// can later be better dealt with using the `FillMissingData` tool after interpolation. +/// +/// # See Also +/// `LidarIdwInterpolation`, `LidarNearestNeighbourGridding`, `LidarTINGridding`, `FilterLidarClasses`, `FillMissingData` pub struct LidarTINGridding { name: String, description: String, diff --git a/src/tools/terrain_analysis/horizon_angle.rs b/src/tools/terrain_analysis/horizon_angle.rs index 1850f0358..00fa136b4 100644 --- a/src/tools/terrain_analysis/horizon_angle.rs +++ b/src/tools/terrain_analysis/horizon_angle.rs @@ -41,6 +41,14 @@ use std::thread; /// /// Ray-tracing is a highly computationally intensive task and therefore this tool may take considerable time to operate for /// larger sized DEMs. Maximum upwind slope is best displayed using a Grey scale palette that is inverted. +/// +/// Horizon angle is best visualized using a white-to-black palette and rescaled from approximately -10 to 70 (see below for +/// an example of horizon angle calculated at a 150-degree azimuth). +/// +/// ![](../../doc_img/HorizonAngle.png) +/// +/// # See Also +/// `TimeInDaylight` pub struct HorizonAngle { name: String, description: String, diff --git a/src/tools/terrain_analysis/map_otos.rs b/src/tools/terrain_analysis/map_otos.rs index 6ab389c8e..ee7c1715e 100644 --- a/src/tools/terrain_analysis/map_otos.rs +++ b/src/tools/terrain_analysis/map_otos.rs @@ -14,10 +14,17 @@ use std::io::{Error, ErrorKind}; use std::path; use std::{f32, f64}; -/// This tool can be used to map off-terrain objects in a digital elevation model (DEM) based on cell-to-cell differences -/// in elevations. +/// This tool can be used to map off-terrain objects in a digital surface model (DSM) based on cell-to-cell differences +/// in elevations and local slopes. The algorithm works by using a region-growing operation to connect neighbouring grid +/// cells outwards from seed cells. Two neighbouring cells are considered connected if the slope between the two cells +/// is less than the user-specified maximum slope value (`--max_slope`). Mapped segments that are less than the minimum +/// feature size (`--min_size`), in grid cells, are assigned a common background value. Note that this method of mapping +/// off-terrain objects, and thereby separating ground cells from non-ground objects in DSMs, works best with fine-resolution +/// DSMs that have been interpolated using a non-smoothing method, such as triangulation (TINing) or nearest-neighbour +/// interpolation. /// /// # See Also +/// `RemoveOffTerrainObjects` pub struct MapOffTerrainObjects { name: String, description: String, diff --git a/src/tools/terrain_analysis/max_diff_from_mean.rs b/src/tools/terrain_analysis/max_diff_from_mean.rs index eb044ba06..c7136abed 100644 --- a/src/tools/terrain_analysis/max_diff_from_mean.rs +++ b/src/tools/terrain_analysis/max_diff_from_mean.rs @@ -307,7 +307,7 @@ impl WhiteboxTool for MaxDifferenceFromMean { let num_loops = (max_scale - min_scale) / step; let mut loop_num = 0; - for midpoint in (min_scale..max_scale).filter(|s| (s - min_scale) % step == 0) { + for midpoint in (min_scale..=max_scale).filter(|s| (s - min_scale) % step == 0) { // .step_by(step) { once step_by is stabilized loop_num += 1; let (tx, rx) = mpsc::channel(); diff --git a/src/tools/terrain_analysis/max_elev_deviation.rs b/src/tools/terrain_analysis/max_elev_deviation.rs index a3ae6c24d..4aa38c968 100644 --- a/src/tools/terrain_analysis/max_elev_deviation.rs +++ b/src/tools/terrain_analysis/max_elev_deviation.rs @@ -357,7 +357,7 @@ impl WhiteboxTool for MaxElevationDeviation { let num_loops = (max_scale - min_scale) / step; let mut loop_num = 0; - for midpoint in (min_scale..max_scale).filter(|s| (s - min_scale) % step == 0) { + for midpoint in (min_scale..=max_scale).filter(|s| (s - min_scale) % step == 0) { // .step_by(step) { once step_by is stabilized loop_num += 1; let (tx, rx) = mpsc::channel(); diff --git a/src/tools/terrain_analysis/remove_off_terrain_objects.rs b/src/tools/terrain_analysis/remove_off_terrain_objects.rs index 19caf781b..3b6c9aec7 100644 --- a/src/tools/terrain_analysis/remove_off_terrain_objects.rs +++ b/src/tools/terrain_analysis/remove_off_terrain_objects.rs @@ -39,7 +39,7 @@ use std::thread; /// models. Available online, DOI: [10.13140/RG.2.2.21226.62401](https://www.researchgate.net/publication/323003064_A_new_method_for_the_removal_of_off-terrain_objects_from_LiDAR-derived_raster_surface_models) /// /// # See Also -/// `TophatTransform`, `LidarGroundPointFilter` +/// `MapOffTerrainObjects`, `TophatTransform`, `LidarGroundPointFilter` pub struct RemoveOffTerrainObjects { name: String, description: String, diff --git a/src/tools/terrain_analysis/time_in_daylight.rs b/src/tools/terrain_analysis/time_in_daylight.rs index f8dc991ee..94d1aec82 100644 --- a/src/tools/terrain_analysis/time_in_daylight.rs +++ b/src/tools/terrain_analysis/time_in_daylight.rs @@ -22,7 +22,41 @@ use rayon::prelude::*; use chrono::prelude::*; use chrono::{Date, FixedOffset, NaiveTime, TimeZone}; -/// This tool calculates the proportion of time a location is within daylight (i.e. outside of an area of shadow cast by a local object). +/// This tool calculates the proportion of time a location is within daylight. That is, it calculates the +/// proportion of time, during a user-defined time frame, that a grid cell in an input digital elevation +/// model (`--dem`) is outside of an area of shadow cast by a local object. The input DEM should truly be +/// a digital surface model (DSM) that contains significant off-terrain objects. Such a model, for example, +/// could be created using the first-return points of a LiDAR data set, or using the `LidarDigitalSurfaceModel` +/// tool. +/// +/// The tool operates by calculating a solar almanac, which estimates the sun's position for the location, in +/// latitude and longitude coordinate (`--lat`, `--long`), of the input DSM. The algorithm then calculates +/// horizon angle (see `HorizonAngle`) rasters from the DSM based on the user-specified azimuth fraction (`--az_fraction`). +/// For example, if an azimuth fraction of 15-degrees is specified, horizon angle rasters could be calculated for +/// the solar azimuths 0, 15, 30, 45... In reality, horizon angle rasters are only calculated for azimuths for which +/// the sun is above the horizon for some time during the tested time period. A horizon angle raster evaluates +/// the vertical angle between each grid cell in a DSM and a distant obstacle (e.g. a mountain ridge, building, tree, etc.) that +/// blocks the view along a specified direction. In calculating horizon angle, the user must specify the maximum search +/// distance (`--max_dist`) beyond which the query for higher, more distant objects will cease. This parameter strongly +/// impacts the performance of the tool, with larger values resulting in significantly longer run-times. Users are advised +/// to set the `--max_dist` based on the maximum shadow length expected in an area. For example, in a relatively flat +/// urban landscape, the tallest building will likely determine the longest shadow lengths. All grid cells for which the +/// calculated solar positions throughout the time frame are higher than the cell's horizon angle are deemed to be +/// illuminated during the time the sun is in the corresponding azimuth fraction. +/// +/// By default, the tool calculates time-in-daylight for a time-frame spanning an entire year. That is, the solar almanac +/// is calculated for each hour, at 10-second intervals, and for each day of the year. Users may alternatively restrict the +/// time of year over which time-in-daylight is calculated by specifying a starting day (1-365; `--start_day`) and ending day +/// (1-365; `--end_day`). Similarly, by specifying start time (`--start_time`) and end time (`--end_time`) parameters, +/// the user is able to measure time-in-daylight for specific ranges of the day (e.g. for the morning or afternoon hours). +/// These time parameters must be specified in 24-hour time (HH:MM:SS), e.g. 15:30:00. `sunrise` and `sunset` are also +/// acceptable inputs for the start time and end time respectively. The timing of sunrise and sunset on each day in the +/// tested time-frame will be determined using the solar almanac. +/// +/// ![](../../doc_img/TimeInDaylight.png) +/// +/// # See Also +/// `LidarDigitalSurfaceModel`, `HorizonAngle` pub struct TimeInDaylight { name: String, description: String, @@ -37,7 +71,7 @@ impl TimeInDaylight { let name = "TimeInDaylight".to_string(); let toolbox = "Geomorphometric Analysis".to_string(); let description = - "Calculates the proportion of time a location is within an area of shadow." + "Calculates the proportion of time a location is not within an area of shadow." .to_string(); let mut parameters = vec![]; diff --git a/whitebox_tools.py b/whitebox_tools.py index 9623c8ac7..ad9a16bf3 100644 --- a/whitebox_tools.py +++ b/whitebox_tools.py @@ -421,6 +421,7 @@ def list_tools(self, keywords=[]): + ############## @@ -1103,6 +1104,24 @@ def extract_raster_values_at_points(self, inputs, points, out_text=False, callba if out_text: args.append("--out_text") return self.run_tool('extract_raster_values_at_points', args, callback) # returns 1 if error + def filter_raster_features_by_area(self, i, output, threshold, background="zero", callback=None): + """Removes small-area features from a raster. + + Keyword arguments: + + i -- Input raster file. + output -- Output raster file. + threshold -- Remove features with fewer grid cells than this threshold value. + background -- Background value. + callback -- Custom function for handling tool text outputs. + """ + args = [] + args.append("--input='{}'".format(i)) + args.append("--output='{}'".format(output)) + args.append("--threshold='{}'".format(threshold)) + args.append("--background={}".format(background)) + return self.run_tool('filter_raster_features_by_area', args, callback) # returns 1 if error + def find_lowest_or_highest_points(self, i, output, out_type="lowest", callback=None): """Locates the lowest and/or highest valued cells in a raster. @@ -2297,7 +2316,7 @@ def shape_complexity_index_raster(self, i, output, callback=None): # Geomorphometric Analysis # ############################ - def aspect(self, dem, output, zfactor=1.0, callback=None): + def aspect(self, dem, output, zfactor=None, callback=None): """Calculates an aspect raster from an input DEM. Keyword arguments: @@ -2310,7 +2329,7 @@ def aspect(self, dem, output, zfactor=1.0, callback=None): args = [] args.append("--dem='{}'".format(dem)) args.append("--output='{}'".format(output)) - args.append("--zfactor={}".format(zfactor)) + if zfactor is not None: args.append("--zfactor='{}'".format(zfactor)) return self.run_tool('aspect', args, callback) # returns 1 if error def average_normal_vector_angular_deviation(self, dem, output, filter=11, callback=None): @@ -2465,7 +2484,7 @@ def downslope_index(self, dem, output, drop=2.0, out_type="tangent", callback=No args.append("--out_type={}".format(out_type)) return self.run_tool('downslope_index', args, callback) # returns 1 if error - def edge_density(self, dem, output, filter=11, norm_diff=5.0, zfactor=1.0, callback=None): + def edge_density(self, dem, output, filter=11, norm_diff=5.0, zfactor=None, callback=None): """Calculates the density of edges, or breaks-in-slope within DEMs. Keyword arguments: @@ -2482,7 +2501,7 @@ def edge_density(self, dem, output, filter=11, norm_diff=5.0, zfactor=1.0, callb args.append("--output='{}'".format(output)) args.append("--filter={}".format(filter)) args.append("--norm_diff={}".format(norm_diff)) - args.append("--zfactor={}".format(zfactor)) + if zfactor is not None: args.append("--zfactor='{}'".format(zfactor)) return self.run_tool('edge_density', args, callback) # returns 1 if error def elev_above_pit(self, dem, output, callback=None): @@ -2549,7 +2568,7 @@ def elev_relative_to_watershed_min_max(self, dem, watersheds, output, callback=N args.append("--output='{}'".format(output)) return self.run_tool('elev_relative_to_watershed_min_max', args, callback) # returns 1 if error - def feature_preserving_smoothing(self, dem, output, filter=11, norm_diff=15.0, num_iter=3, max_diff=0.5, zfactor=1.0, callback=None): + def feature_preserving_smoothing(self, dem, output, filter=11, norm_diff=15.0, num_iter=3, max_diff=0.5, zfactor=None, callback=None): """Reduces short-scale variation in an input DEM using a modified Sun et al. (2007) algorithm. Keyword arguments: @@ -2570,7 +2589,7 @@ def feature_preserving_smoothing(self, dem, output, filter=11, norm_diff=15.0, n args.append("--norm_diff={}".format(norm_diff)) args.append("--num_iter={}".format(num_iter)) args.append("--max_diff={}".format(max_diff)) - args.append("--zfactor={}".format(zfactor)) + if zfactor is not None: args.append("--zfactor='{}'".format(zfactor)) return self.run_tool('feature_preserving_smoothing', args, callback) # returns 1 if error def fetch_analysis(self, dem, output, azimuth=0.0, hgt_inc=0.05, callback=None): @@ -2627,7 +2646,7 @@ def find_ridges(self, dem, output, line_thin=True, callback=None): if line_thin: args.append("--line_thin") return self.run_tool('find_ridges', args, callback) # returns 1 if error - def hillshade(self, dem, output, azimuth=315.0, altitude=30.0, zfactor=1.0, callback=None): + def hillshade(self, dem, output, azimuth=315.0, altitude=30.0, zfactor=None, callback=None): """Calculates a hillshade raster from an input DEM. Keyword arguments: @@ -2644,25 +2663,25 @@ def hillshade(self, dem, output, azimuth=315.0, altitude=30.0, zfactor=1.0, call args.append("--output='{}'".format(output)) args.append("--azimuth={}".format(azimuth)) args.append("--altitude={}".format(altitude)) - args.append("--zfactor={}".format(zfactor)) + if zfactor is not None: args.append("--zfactor='{}'".format(zfactor)) return self.run_tool('hillshade', args, callback) # returns 1 if error - def horizon_angle(self, dem, output, azimuth=0.0, max_dist=None, callback=None): + def horizon_angle(self, dem, output, azimuth=0.0, max_dist=100.0, callback=None): """Calculates horizon angle (maximum upwind slope) for each grid cell in an input DEM. Keyword arguments: dem -- Input raster DEM file. output -- Output raster file. - azimuth -- Wind azimuth in degrees. - max_dist -- Optional maximum search distance (unspecified if none; in xy units). + azimuth -- Azimuth, in degrees. + max_dist -- Optional maximum search distance (unspecified if none; in xy units). Minimum value is 5 x cell size. callback -- Custom function for handling tool text outputs. """ args = [] args.append("--dem='{}'".format(dem)) args.append("--output='{}'".format(output)) args.append("--azimuth={}".format(azimuth)) - if max_dist is not None: args.append("--max_dist='{}'".format(max_dist)) + args.append("--max_dist={}".format(max_dist)) return self.run_tool('horizon_angle', args, callback) # returns 1 if error def hypsometric_analysis(self, inputs, output, watershed=None, callback=None): @@ -2681,7 +2700,7 @@ def hypsometric_analysis(self, inputs, output, watershed=None, callback=None): args.append("--output='{}'".format(output)) return self.run_tool('hypsometric_analysis', args, callback) # returns 1 if error - def hypsometrically_tinted_hillshade(self, dem, output, altitude=45.0, hs_weight=0.5, brightness=0.5, atmospheric=0.0, palette="atlas", reverse=False, zfactor=1.0, full_mode=False, callback=None): + def hypsometrically_tinted_hillshade(self, dem, output, altitude=45.0, hs_weight=0.5, brightness=0.5, atmospheric=0.0, palette="atlas", reverse=False, zfactor=None, full_mode=False, callback=None): """Creates an colour shaded relief image from an input DEM. Keyword arguments: @@ -2692,7 +2711,7 @@ def hypsometrically_tinted_hillshade(self, dem, output, altitude=45.0, hs_weight hs_weight -- Weight given to hillshade relative to relief (0.0-1.0). brightness -- Brightness factor (0.0-1.0). atmospheric -- Atmospheric effects weight (0.0-1.0). - palette -- Options include 'atlas', 'high_relief', 'arid', 'soft', 'muted', 'purple', 'viridi', 'gn_yl', 'pi_y_g', 'bl_yl_rd', and 'deep. + palette -- Options include 'atlas', 'high_relief', 'arid', 'soft', 'muted', 'purple', 'viridi', 'gn_yl', 'pi_y_g', 'bl_yl_rd', and 'deep'. reverse -- Optional flag indicating whether to use reverse the palette. zfactor -- Optional multiplier for when the vertical and horizontal units are not the same. full_mode -- Optional flag indicating whether to use full 360-degrees of illumination sources. @@ -2707,10 +2726,28 @@ def hypsometrically_tinted_hillshade(self, dem, output, altitude=45.0, hs_weight args.append("--atmospheric={}".format(atmospheric)) args.append("--palette={}".format(palette)) if reverse: args.append("--reverse") - args.append("--zfactor={}".format(zfactor)) + if zfactor is not None: args.append("--zfactor='{}'".format(zfactor)) if full_mode: args.append("--full_mode") return self.run_tool('hypsometrically_tinted_hillshade', args, callback) # returns 1 if error + def map_off_terrain_objects(self, dem, output, max_slope=40.0, min_size=1, callback=None): + """Maps off-terrain objects in a digital elevation model (DEM). + + Keyword arguments: + + dem -- Input raster DEM file. + output -- Output raster file. + max_slope -- Maximum inter-cell absolute slope. + min_size -- Minimum feature size, in grid cells. + callback -- Custom function for handling tool text outputs. + """ + args = [] + args.append("--dem='{}'".format(dem)) + args.append("--output='{}'".format(output)) + args.append("--max_slope={}".format(max_slope)) + args.append("--min_size={}".format(min_size)) + return self.run_tool('map_off_terrain_objects', args, callback) # returns 1 if error + def max_anisotropy_dev(self, dem, out_mag, out_scale, max_scale, min_scale=3, step=2, callback=None): """Calculates the maximum anisotropy (directionality) in elevation deviation over a range of spatial scales. @@ -2865,7 +2902,7 @@ def min_downslope_elev_change(self, dem, output, callback=None): args.append("--output='{}'".format(output)) return self.run_tool('min_downslope_elev_change', args, callback) # returns 1 if error - def multidirectional_hillshade(self, dem, output, altitude=45.0, zfactor=1.0, full_mode=False, callback=None): + def multidirectional_hillshade(self, dem, output, altitude=45.0, zfactor=None, full_mode=False, callback=None): """Calculates a multi-direction hillshade raster from an input DEM. Keyword arguments: @@ -2881,7 +2918,7 @@ def multidirectional_hillshade(self, dem, output, altitude=45.0, zfactor=1.0, fu args.append("--dem='{}'".format(dem)) args.append("--output='{}'".format(output)) args.append("--altitude={}".format(altitude)) - args.append("--zfactor={}".format(zfactor)) + if zfactor is not None: args.append("--zfactor='{}'".format(zfactor)) if full_mode: args.append("--full_mode") return self.run_tool('multidirectional_hillshade', args, callback) # returns 1 if error @@ -3051,7 +3088,7 @@ def num_upslope_neighbours(self, dem, output, callback=None): args.append("--output='{}'".format(output)) return self.run_tool('num_upslope_neighbours', args, callback) # returns 1 if error - def pennock_landform_class(self, dem, output, slope=3.0, prof=0.1, plan=0.0, zfactor=1.0, callback=None): + def pennock_landform_class(self, dem, output, slope=3.0, prof=0.1, plan=0.0, zfactor=None, callback=None): """Classifies hillslope zones based on slope, profile curvature, and plan curvature. Keyword arguments: @@ -3070,7 +3107,7 @@ def pennock_landform_class(self, dem, output, slope=3.0, prof=0.1, plan=0.0, zfa args.append("--slope={}".format(slope)) args.append("--prof={}".format(prof)) args.append("--plan={}".format(plan)) - args.append("--zfactor={}".format(zfactor)) + if zfactor is not None: args.append("--zfactor='{}'".format(zfactor)) return self.run_tool('pennock_landform_class', args, callback) # returns 1 if error def percent_elev_range(self, dem, output, filterx=3, filtery=3, callback=None): @@ -3091,7 +3128,7 @@ def percent_elev_range(self, dem, output, filterx=3, filtery=3, callback=None): args.append("--filtery={}".format(filtery)) return self.run_tool('percent_elev_range', args, callback) # returns 1 if error - def plan_curvature(self, dem, output, zfactor=1.0, callback=None): + def plan_curvature(self, dem, output, zfactor=None, callback=None): """Calculates a plan (contour) curvature raster from an input DEM. Keyword arguments: @@ -3104,7 +3141,7 @@ def plan_curvature(self, dem, output, zfactor=1.0, callback=None): args = [] args.append("--dem='{}'".format(dem)) args.append("--output='{}'".format(output)) - args.append("--zfactor={}".format(zfactor)) + if zfactor is not None: args.append("--zfactor='{}'".format(zfactor)) return self.run_tool('plan_curvature', args, callback) # returns 1 if error def profile(self, lines, surface, output, callback=None): @@ -3123,7 +3160,7 @@ def profile(self, lines, surface, output, callback=None): args.append("--output='{}'".format(output)) return self.run_tool('profile', args, callback) # returns 1 if error - def profile_curvature(self, dem, output, zfactor=1.0, callback=None): + def profile_curvature(self, dem, output, zfactor=None, callback=None): """Calculates a profile curvature raster from an input DEM. Keyword arguments: @@ -3136,10 +3173,10 @@ def profile_curvature(self, dem, output, zfactor=1.0, callback=None): args = [] args.append("--dem='{}'".format(dem)) args.append("--output='{}'".format(output)) - args.append("--zfactor={}".format(zfactor)) + if zfactor is not None: args.append("--zfactor='{}'".format(zfactor)) return self.run_tool('profile_curvature', args, callback) # returns 1 if error - def relative_aspect(self, dem, output, azimuth=0.0, zfactor=1.0, callback=None): + def relative_aspect(self, dem, output, azimuth=0.0, zfactor=None, callback=None): """Calculates relative aspect (relative to a user-specified direction) from an input DEM. Keyword arguments: @@ -3154,7 +3191,7 @@ def relative_aspect(self, dem, output, azimuth=0.0, zfactor=1.0, callback=None): args.append("--dem='{}'".format(dem)) args.append("--output='{}'".format(output)) args.append("--azimuth={}".format(azimuth)) - args.append("--zfactor={}".format(zfactor)) + if zfactor is not None: args.append("--zfactor='{}'".format(zfactor)) return self.run_tool('relative_aspect', args, callback) # returns 1 if error def relative_topographic_position(self, dem, output, filterx=11, filtery=11, callback=None): @@ -3193,7 +3230,7 @@ def remove_off_terrain_objects(self, dem, output, filter=11, slope=15.0, callbac args.append("--slope={}".format(slope)) return self.run_tool('remove_off_terrain_objects', args, callback) # returns 1 if error - def ruggedness_index(self, dem, output, zfactor=1.0, callback=None): + def ruggedness_index(self, dem, output, zfactor=None, callback=None): """Calculates the Riley et al.'s (1999) terrain ruggedness index from an input DEM. Keyword arguments: @@ -3206,7 +3243,7 @@ def ruggedness_index(self, dem, output, zfactor=1.0, callback=None): args = [] args.append("--dem='{}'".format(dem)) args.append("--output='{}'".format(output)) - args.append("--zfactor={}".format(zfactor)) + if zfactor is not None: args.append("--zfactor='{}'".format(zfactor)) return self.run_tool('ruggedness_index', args, callback) # returns 1 if error def sediment_transport_index(self, sca, slope, output, sca_exponent=0.4, slope_exponent=1.3, callback=None): @@ -3229,7 +3266,7 @@ def sediment_transport_index(self, sca, slope, output, sca_exponent=0.4, slope_e args.append("--slope_exponent={}".format(slope_exponent)) return self.run_tool('sediment_transport_index', args, callback) # returns 1 if error - def slope(self, dem, output, zfactor=1.0, units="degrees", callback=None): + def slope(self, dem, output, zfactor=None, units="degrees", callback=None): """Calculates a slope raster from an input DEM. Keyword arguments: @@ -3243,7 +3280,7 @@ def slope(self, dem, output, zfactor=1.0, units="degrees", callback=None): args = [] args.append("--dem='{}'".format(dem)) args.append("--output='{}'".format(output)) - args.append("--zfactor={}".format(zfactor)) + if zfactor is not None: args.append("--zfactor='{}'".format(zfactor)) args.append("--units={}".format(units)) return self.run_tool('slope', args, callback) # returns 1 if error @@ -3279,7 +3316,7 @@ def spherical_std_dev_of_normals(self, dem, output, filter=11, callback=None): args.append("--filter={}".format(filter)) return self.run_tool('spherical_std_dev_of_normals', args, callback) # returns 1 if error - def standard_deviation_of_slope(self, i, output, zfactor=1.0, filterx=11, filtery=11, callback=None): + def standard_deviation_of_slope(self, i, output, zfactor=None, filterx=11, filtery=11, callback=None): """Calculates the standard deviation of slope from an input DEM. Keyword arguments: @@ -3294,7 +3331,7 @@ def standard_deviation_of_slope(self, i, output, zfactor=1.0, filterx=11, filter args = [] args.append("--input='{}'".format(i)) args.append("--output='{}'".format(output)) - args.append("--zfactor={}".format(zfactor)) + if zfactor is not None: args.append("--zfactor='{}'".format(zfactor)) args.append("--filterx={}".format(filterx)) args.append("--filtery={}".format(filtery)) return self.run_tool('standard_deviation_of_slope', args, callback) # returns 1 if error @@ -3331,7 +3368,7 @@ def surface_area_ratio(self, dem, output, callback=None): args.append("--output='{}'".format(output)) return self.run_tool('surface_area_ratio', args, callback) # returns 1 if error - def tangential_curvature(self, dem, output, zfactor=1.0, callback=None): + def tangential_curvature(self, dem, output, zfactor=None, callback=None): """Calculates a tangential curvature raster from an input DEM. Keyword arguments: @@ -3344,10 +3381,42 @@ def tangential_curvature(self, dem, output, zfactor=1.0, callback=None): args = [] args.append("--dem='{}'".format(dem)) args.append("--output='{}'".format(output)) - args.append("--zfactor={}".format(zfactor)) + if zfactor is not None: args.append("--zfactor='{}'".format(zfactor)) return self.run_tool('tangential_curvature', args, callback) # returns 1 if error - def total_curvature(self, dem, output, zfactor=1.0, callback=None): + def time_in_daylight(self, dem, output, lat, long, az_fraction=10.0, max_dist=100.0, utc_offset="00:00", start_day=1, end_day=365, start_time="00:00:00", end_time="23:59:59", callback=None): + """Calculates the proportion of time a location is within an area of shadow. + + Keyword arguments: + + dem -- Input raster DEM file. + output -- Output raster file. + az_fraction -- Azimuth fraction in degrees. + max_dist -- Optional maximum search distance. Minimum value is 5 x cell size. + lat -- Centre point latitude. + long -- Centre point longitude. + utc_offset -- UTC time offset, in hours (e.g. -04:00, +06:00). + start_day -- Start day of the year (1-365). + end_day -- End day of the year (1-365). + start_time -- Starting hour to track shadows (e.g. 5, 5:00, 05:00:00). Assumes 24-hour time: HH:MM:SS. 'sunrise' is also a valid time. + end_time -- Starting hour to track shadows (e.g. 21, 21:00, 21:00:00). Assumes 24-hour time: HH:MM:SS. 'sunset' is also a valid time. + callback -- Custom function for handling tool text outputs. + """ + args = [] + args.append("--dem='{}'".format(dem)) + args.append("--output='{}'".format(output)) + args.append("--az_fraction={}".format(az_fraction)) + args.append("--max_dist={}".format(max_dist)) + args.append("--lat='{}'".format(lat)) + args.append("--long='{}'".format(long)) + args.append("--utc_offset={}".format(utc_offset)) + args.append("--start_day={}".format(start_day)) + args.append("--end_day={}".format(end_day)) + args.append("--start_time={}".format(start_time)) + args.append("--end_time={}".format(end_time)) + return self.run_tool('time_in_daylight', args, callback) # returns 1 if error + + def total_curvature(self, dem, output, zfactor=None, callback=None): """Calculates a total curvature raster from an input DEM. Keyword arguments: @@ -3360,7 +3429,7 @@ def total_curvature(self, dem, output, zfactor=1.0, callback=None): args = [] args.append("--dem='{}'".format(dem)) args.append("--output='{}'".format(output)) - args.append("--zfactor={}".format(zfactor)) + if zfactor is not None: args.append("--zfactor='{}'".format(zfactor)) return self.run_tool('total_curvature', args, callback) # returns 1 if error def viewshed(self, dem, stations, output, height=2.0, callback=None): @@ -4025,7 +4094,7 @@ def insert_dams(self, dem, dam_pts, output, damlength, callback=None): 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): + def isobasins(self, dem, output, size, connections=False, callback=None): """Divides a landscape into nearly equal sized drainage basins (i.e. watersheds). Keyword arguments: @@ -4033,12 +4102,14 @@ def isobasins(self, dem, output, size, callback=None): dem -- Input raster DEM file. output -- Output raster file. size -- Target basin size, in grid cells. + connections -- Output upstream-downstream flow connections among basins?. callback -- Custom function for handling tool text outputs. """ args = [] args.append("--dem='{}'".format(dem)) args.append("--output='{}'".format(output)) args.append("--size='{}'".format(size)) + if connections: args.append("--connections") return self.run_tool('isobasins', args, callback) # returns 1 if error def jenson_snap_pour_points(self, pour_pts, streams, output, snap_dist, callback=None): @@ -4611,19 +4682,23 @@ def remove_spurs(self, i, output, iterations=10, callback=None): args.append("--iterations={}".format(iterations)) return self.run_tool('remove_spurs', args, callback) # returns 1 if error - def resample(self, inputs, destination, method="cc", callback=None): + def resample(self, inputs, output, cell_size=None, base=None, method="cc", callback=None): """Resamples one or more input images into a destination image. Keyword arguments: inputs -- Input raster files. - destination -- Destination raster file. + 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. method -- Resampling method; options include 'nn' (nearest neighbour), 'bilinear', and 'cc' (cubic convolution). callback -- Custom function for handling tool text outputs. """ args = [] args.append("--inputs='{}'".format(inputs)) - args.append("--destination='{}'".format(destination)) + 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)) args.append("--method={}".format(method)) return self.run_tool('resample', args, callback) # returns 1 if error @@ -5845,6 +5920,30 @@ def lidar_colourize(self, in_lidar, in_image, output, callback=None): args.append("--output='{}'".format(output)) return self.run_tool('lidar_colourize', args, callback) # returns 1 if error + def lidar_digital_surface_model(self, i=None, output=None, resolution=1.0, radius=0.5, minz=None, maxz=None, max_triangle_edge_length=None, callback=None): + """Creates a top-surface digital surface model (DSM) from a LiDAR point cloud. + + Keyword arguments: + + i -- Input LiDAR file (including extension). + output -- Output raster file (including extension). + resolution -- Output raster's grid resolution. + radius -- Search Radius. + minz -- Optional minimum elevation for inclusion in interpolation. + maxz -- Optional maximum elevation for inclusion in interpolation. + 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 = [] + if i is not None: args.append("--input='{}'".format(i)) + if output is not None: args.append("--output='{}'".format(output)) + args.append("--resolution={}".format(resolution)) + args.append("--radius={}".format(radius)) + if minz is not None: args.append("--minz='{}'".format(minz)) + if maxz is not None: args.append("--maxz='{}'".format(maxz)) + if max_triangle_edge_length is not None: args.append("--max_triangle_edge_length='{}'".format(max_triangle_edge_length)) + return self.run_tool('lidar_digital_surface_model', args, callback) # returns 1 if error + def lidar_elevation_slice(self, i, output, minz=None, maxz=None, cls=False, inclassval=2, outclassval=1, callback=None): """Outputs all of the points within a LiDAR (LAS) point file that lie between a specified elevation range. @@ -6115,7 +6214,7 @@ def lidar_point_stats(self, i=None, resolution=1.0, num_points=True, num_pulses= if predom_class: args.append("--predom_class") return self.run_tool('lidar_point_stats', args, callback) # returns 1 if error - def lidar_ransac_planes(self, i, output, radius=2.0, num_iter=50, num_samples=5, threshold=0.35, model_size=8, max_slope=80.0, classify=False, callback=None): + def lidar_ransac_planes(self, i, output, radius=2.0, num_iter=50, num_samples=5, threshold=0.35, model_size=8, max_slope=80.0, classify=False, last_returns=False, callback=None): """Performs a RANSAC analysis to identify points within a LiDAR point cloud that belong to linear planes. Keyword arguments: @@ -6129,6 +6228,7 @@ def lidar_ransac_planes(self, i, output, radius=2.0, num_iter=50, num_samples=5, model_size -- Acceptable model size. max_slope -- Maximum planar slope. classify -- Classify points as ground (2) or off-ground (1). + last_returns -- Only include last- and only-return points. callback -- Custom function for handling tool text outputs. """ args = [] @@ -6141,6 +6241,7 @@ def lidar_ransac_planes(self, i, output, radius=2.0, num_iter=50, num_samples=5, args.append("--model_size={}".format(model_size)) args.append("--max_slope={}".format(max_slope)) if classify: args.append("--classify") + if last_returns: args.append("--last_returns") return self.run_tool('lidar_ransac_planes', args, callback) # returns 1 if error def lidar_rbf_interpolation(self, i=None, output=None, parameter="elevation", returns="all", resolution=1.0, num_points=20, exclude_cls=None, minz=None, maxz=None, func_type="ThinPlateSpline", poly_order="none", weight=5, callback=None): @@ -6226,11 +6327,11 @@ def lidar_rooftop_analysis(self, buildings, output, i=None, radius=2.0, num_iter radius -- Search Radius. num_iter -- Number of iterations. num_samples -- Number of sample points on which to build the model. - threshold -- Threshold used to determine inlier points. - model_size -- Acceptable model size. - max_slope -- Maximum planar slope. + threshold -- Threshold used to determine inlier points (in elevation units). + model_size -- Acceptable model size, in points. + max_slope -- Maximum planar slope, in degrees. norm_diff -- Maximum difference in normal vectors, in degrees. - azimuth -- Illumination source azimuth in degrees. + azimuth -- Illumination source azimuth, in degrees. altitude -- Illumination source altitude in degrees. callback -- Custom function for handling tool text outputs. """ @@ -6383,7 +6484,7 @@ def lidar_tile_footprint(self, output, i=None, hull=False, callback=None): if hull: args.append("--hull") return self.run_tool('lidar_tile_footprint', args, callback) # returns 1 if error - def lidar_tin_gridding(self, i=None, output=None, parameter="elevation", returns="all", resolution=1.0, exclude_cls=None, minz=None, maxz=None, max_triangle_edge_length=None, callback=None): + def lidar_tin_gridding(self, i=None, output=None, parameter="elevation", returns="all", resolution=1.0, exclude_cls="7,18", minz=None, maxz=None, max_triangle_edge_length=None, callback=None): """Creates a raster grid based on a Delaunay triangular irregular network (TIN) fitted to LiDAR points. Keyword arguments: @@ -6405,7 +6506,7 @@ def lidar_tin_gridding(self, i=None, output=None, parameter="elevation", returns args.append("--parameter={}".format(parameter)) args.append("--returns={}".format(returns)) args.append("--resolution={}".format(resolution)) - if exclude_cls is not None: args.append("--exclude_cls='{}'".format(exclude_cls)) + args.append("--exclude_cls={}".format(exclude_cls)) if minz is not None: args.append("--minz='{}'".format(minz)) if maxz is not None: args.append("--maxz='{}'".format(maxz)) if max_triangle_edge_length is not None: args.append("--max_triangle_edge_length='{}'".format(max_triangle_edge_length))