Fixes for #107 and #109

Cargo.toml

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

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.

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::<Vec<_>>();
+
+                            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::<Vec<_>>();
+                            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;
                                 }
                             }
                         }

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,

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,