-
Notifications
You must be signed in to change notification settings - Fork 161
/
stream_link_id.rs
executable file
·438 lines (400 loc) · 16.7 KB
/
stream_link_id.rs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
/*
This tool is part of the WhiteboxTools geospatial analysis library.
Authors: Dr. John Lindsay
Created: 25/06/2017
Last Modified: 03/09/2020
License: MIT
*/
use whitebox_raster::*;
use whitebox_common::structures::Array2D;
use crate::tools::*;
use std::env;
use std::f64;
use std::io::{Error, ErrorKind};
use std::path;
/// This tool can be used to assign each link in a stream network a unique numeric identifier. This grid is
/// used by a number of other stream network analysis tools.
///
/// The input streams file (`--streams`) is used to designate which grid cells contain a stream and the pointer
/// image is used to traverse the stream network. Stream cells are designated in the streams image as all values
/// greater than zero. Thus, all non-stream or background grid cells are commonly assigned either zeros or NoData
/// values. Background cells will be assigned the NoData value in the output image, unless the `--zero_background`
/// parameter is used, in which case non-stream cells will be assigned zero values in the output.
///
/// The user must specify the name of a flow pointer (flow direction) raster (`--d8_pntr`), a streams raster
/// (`--streams`), and the output raster (`--output`). The flow pointer and streams rasters should be generated
/// using the `D8Pointer` algorithm. This will require a depressionless DEM, processed using either the
/// `BreachDepressions` or `FillDepressions` tool.
///
/// By default, the pointer raster is assumed to use the clockwise indexing method used by WhiteboxTools.
/// If the pointer file contains ESRI flow direction values instead, the `--esri_pntr` parameter must be specified.
///
/// # See Also
/// `D8Pointer`, `TributaryIdentifier`, `BreachDepressions`, `FillDepressions`
pub struct StreamLinkIdentifier {
name: String,
description: String,
toolbox: String,
parameters: Vec<ToolParameter>,
example_usage: String,
}
impl StreamLinkIdentifier {
pub fn new() -> StreamLinkIdentifier {
// public constructor
let name = "StreamLinkIdentifier".to_string();
let toolbox = "Stream Network Analysis".to_string();
let description =
"Assigns a unique identifier to each link in a stream network.".to_string();
let mut parameters = vec![];
parameters.push(ToolParameter {
name: "Input D8 Pointer File".to_owned(),
flags: vec!["--d8_pntr".to_owned()],
description: "Input raster D8 pointer file.".to_owned(),
parameter_type: ParameterType::ExistingFile(ParameterFileType::Raster),
default_value: None,
optional: false,
});
parameters.push(ToolParameter {
name: "Input Streams File".to_owned(),
flags: vec!["--streams".to_owned()],
description: "Input raster streams file.".to_owned(),
parameter_type: ParameterType::ExistingFile(ParameterFileType::Raster),
default_value: None,
optional: false,
});
parameters.push(ToolParameter {
name: "Output File".to_owned(),
flags: vec!["-o".to_owned(), "--output".to_owned()],
description: "Output raster file.".to_owned(),
parameter_type: ParameterType::NewFile(ParameterFileType::Raster),
default_value: None,
optional: false,
});
parameters.push(ToolParameter {
name: "Does the pointer file use the ESRI pointer scheme?".to_owned(),
flags: vec!["--esri_pntr".to_owned()],
description: "D8 pointer uses the ESRI style scheme.".to_owned(),
parameter_type: ParameterType::Boolean,
default_value: Some("false".to_owned()),
optional: true,
});
parameters.push(ToolParameter {
name: "Should a background value of zero be used?".to_owned(),
flags: vec!["--zero_background".to_owned()],
description: "Flag indicating whether a background value of zero should be used."
.to_owned(),
parameter_type: ParameterType::Boolean,
default_value: None,
optional: true,
});
let sep: String = path::MAIN_SEPARATOR.to_string();
let e = format!("{}", env::current_exe().unwrap().display());
let mut parent = env::current_exe().unwrap();
parent.pop();
let p = format!("{}", parent.display());
let mut short_exe = e
.replace(&p, "")
.replace(".exe", "")
.replace(".", "")
.replace(&sep, "");
if e.contains(".exe") {
short_exe += ".exe";
}
let usage = format!(">>.*{0} -r={1} -v --wd=\"*path*to*data*\" --d8_pntr=D8.tif --streams=streams.tif -o=output.tif
>>.*{0} -r={1} -v --wd=\"*path*to*data*\" --d8_pntr=D8.tif --streams=streams.tif -o=output.tif --esri_pntr --zero_background", short_exe, name).replace("*", &sep);
StreamLinkIdentifier {
name: name,
description: description,
toolbox: toolbox,
parameters: parameters,
example_usage: usage,
}
}
}
impl WhiteboxTool for StreamLinkIdentifier {
fn get_source_file(&self) -> String {
String::from(file!())
}
fn get_tool_name(&self) -> String {
self.name.clone()
}
fn get_tool_description(&self) -> String {
self.description.clone()
}
fn get_tool_parameters(&self) -> String {
let mut s = String::from("{\"parameters\": [");
for i in 0..self.parameters.len() {
if i < self.parameters.len() - 1 {
s.push_str(&(self.parameters[i].to_string()));
s.push_str(",");
} else {
s.push_str(&(self.parameters[i].to_string()));
}
}
s.push_str("]}");
s
}
fn get_example_usage(&self) -> String {
self.example_usage.clone()
}
fn get_toolbox(&self) -> String {
self.toolbox.clone()
}
fn run<'a>(
&self,
args: Vec<String>,
working_directory: &'a str,
verbose: bool,
) -> Result<(), Error> {
let mut d8_file = String::new();
let mut streams_file = String::new();
let mut output_file = String::new();
let mut esri_style = false;
let mut background_val = f64::NEG_INFINITY;
if args.len() == 0 {
return Err(Error::new(
ErrorKind::InvalidInput,
"Tool run with no parameters.",
));
}
for i in 0..args.len() {
let mut arg = args[i].replace("\"", "");
arg = arg.replace("\'", "");
let cmd = arg.split("="); // in case an equals sign was used
let vec = cmd.collect::<Vec<&str>>();
let mut keyval = false;
if vec.len() > 1 {
keyval = true;
}
if vec[0].to_lowercase() == "-d8_pntr" || vec[0].to_lowercase() == "--d8_pntr" {
if keyval {
d8_file = vec[1].to_string();
} else {
d8_file = args[i + 1].to_string();
}
} else if vec[0].to_lowercase() == "-streams" || vec[0].to_lowercase() == "--streams" {
if keyval {
streams_file = vec[1].to_string();
} else {
streams_file = args[i + 1].to_string();
}
} else if vec[0].to_lowercase() == "-o" || vec[0].to_lowercase() == "--output" {
if keyval {
output_file = vec[1].to_string();
} else {
output_file = args[i + 1].to_string();
}
} else if vec[0].to_lowercase() == "-esri_pntr"
|| vec[0].to_lowercase() == "--esri_pntr"
|| vec[0].to_lowercase() == "--esri_style"
{
if vec.len() == 1 || !vec[1].to_string().to_lowercase().contains("false") {
esri_style = true;
}
} else if vec[0].to_lowercase() == "-zero_background"
|| vec[0].to_lowercase() == "--zero_background"
{
if vec.len() == 1 || !vec[1].to_string().to_lowercase().contains("false") {
background_val = 0f64;
}
}
}
if verbose {
let tool_name = self.get_tool_name();
let welcome_len = format!("* Welcome to {} *", tool_name).len().max(28);
// 28 = length of the 'Powered by' by statement.
println!("{}", "*".repeat(welcome_len));
println!("* Welcome to {} {}*", tool_name, " ".repeat(welcome_len - 15 - tool_name.len()));
println!("* Powered by WhiteboxTools {}*", " ".repeat(welcome_len - 28));
println!("* www.whiteboxgeo.com {}*", " ".repeat(welcome_len - 23));
println!("{}", "*".repeat(welcome_len));
}
let sep: String = path::MAIN_SEPARATOR.to_string();
let mut progress: usize;
let mut old_progress: usize = 1;
if !d8_file.contains(&sep) && !d8_file.contains("/") {
d8_file = format!("{}{}", working_directory, d8_file);
}
if !streams_file.contains(&sep) && !streams_file.contains("/") {
streams_file = format!("{}{}", working_directory, streams_file);
}
if !output_file.contains(&sep) && !output_file.contains("/") {
output_file = format!("{}{}", working_directory, output_file);
}
if verbose {
println!("Reading pointer data...")
};
let pntr = Raster::new(&d8_file, "r")?;
let pntr_nodata = pntr.configs.nodata;
if verbose {
println!("Reading streams data...")
};
let streams = Raster::new(&streams_file, "r")?;
let start = Instant::now();
let rows = pntr.configs.rows as isize;
let columns = pntr.configs.columns as isize;
let num_cells = pntr.num_cells();
let nodata = streams.configs.nodata;
if background_val == f64::NEG_INFINITY {
background_val = nodata;
}
// make sure the input files have the same size
if streams.configs.rows != pntr.configs.rows
|| streams.configs.columns != pntr.configs.columns
{
return Err(Error::new(
ErrorKind::InvalidInput,
"The input files must have the same number of rows and columns and spatial extent.",
));
}
let mut output = Raster::initialize_using_file(&output_file, &streams);
let mut stack = Vec::with_capacity((rows * columns) as usize);
// calculate the number of inflowing cells
let mut num_inflowing: Array2D<i8> = Array2D::new(rows, columns, -1, -1)?;
let dx = [1, 1, 1, 0, -1, -1, -1, 0];
let dy = [-1, 0, 1, 1, 1, 0, -1, -1];
let mut inflowing_vals = [16f64, 32f64, 64f64, 128f64, 1f64, 2f64, 4f64, 8f64];
if esri_style {
inflowing_vals = [8f64, 16f64, 32f64, 64f64, 128f64, 1f64, 2f64, 4f64];
}
let mut num_solved_cells = 0;
let mut count: i8;
let mut current_id = 1f64;
for row in 0..rows {
for col in 0..columns {
if streams[(row, col)] > 0.0 {
count = 0i8;
for i in 0..8 {
if streams[(row + dy[i], col + dx[i])] > 0.0
&& pntr[(row + dy[i], col + dx[i])] == inflowing_vals[i]
{
count += 1;
}
}
num_inflowing[(row, col)] = count;
if count == 0 {
// It's a headwater; add it to the stack
stack.push((row, col));
output[(row, col)] = current_id;
current_id += 1f64;
}
} else {
if pntr[(row, col)] != pntr_nodata {
output[(row, col)] = background_val;
} else {
output[(row, col)] = nodata;
}
num_solved_cells += 1;
}
}
if verbose {
progress = (100.0_f64 * num_solved_cells as f64 / (num_cells - 1) as f64) as usize;
if progress != old_progress {
println!("Progress: {}%", progress);
old_progress = progress;
}
}
}
// Create a mapping from the pointer values to cells offsets.
// This may seem wasteful, using only 8 of 129 values in the array,
// but the mapping method is far faster than calculating z.ln() / ln(2.0).
// It's also a good way of allowing for different point styles.
let mut pntr_matches: [usize; 129] = [999usize; 129];
if !esri_style {
// This maps Whitebox-style D8 pointer values
// onto the cell offsets in dx and dy.
pntr_matches[1] = 0usize;
pntr_matches[2] = 1usize;
pntr_matches[4] = 2usize;
pntr_matches[8] = 3usize;
pntr_matches[16] = 4usize;
pntr_matches[32] = 5usize;
pntr_matches[64] = 6usize;
pntr_matches[128] = 7usize;
} else {
// This maps Esri-style D8 pointer values
// onto the cell offsets in dx and dy.
pntr_matches[1] = 1usize;
pntr_matches[2] = 2usize;
pntr_matches[4] = 3usize;
pntr_matches[8] = 4usize;
pntr_matches[16] = 5usize;
pntr_matches[32] = 6usize;
pntr_matches[64] = 7usize;
pntr_matches[128] = 0usize;
}
let (mut row, mut col): (isize, isize);
let (mut row_n, mut col_n): (isize, isize);
// let mut cell: (isize, isize);
let mut dir: usize;
let mut val: f64;
let mut c: usize;
while !stack.is_empty() {
let cell = stack.pop().expect("Error during pop operation.");
row = cell.0;
col = cell.1;
val = output[(row, col)];
// find the downstream cell
dir = pntr[(row, col)] as usize;
if dir > 0 {
if dir > 128 || pntr_matches[dir] == 999 {
return Err(Error::new(ErrorKind::InvalidInput,
"An unexpected value has been identified in the pointer image. This tool requires a pointer grid that has been created using either the D8 or Rho8 tools."));
}
c = pntr_matches[dir];
row_n = row + dy[c];
col_n = col + dx[c];
if num_inflowing[(row_n, col_n)] > 1 {
current_id += 1f64;
output[(row_n, col_n)] = current_id;
} else if output[(row_n, col_n)] == nodata {
output[(row_n, col_n)] = val;
}
num_inflowing.decrement(row_n, col_n, 1);
if num_inflowing[(row_n, col_n)] == 0 {
stack.push((row_n, col_n));
}
}
if verbose {
progress = (100.0_f64 * num_solved_cells as f64 / (num_cells - 1) as f64) as usize;
if progress != old_progress {
println!("Progress: {}%", progress);
old_progress = progress;
}
}
}
let elapsed_time = get_formatted_elapsed_time(start);
if current_id < 32767f64 {
output.configs.data_type = DataType::I16;
} else {
output.configs.data_type = DataType::I32;
}
output.configs.palette = "qual.plt".to_string();
output.configs.photometric_interp = PhotometricInterpretation::Categorical;
output.add_metadata_entry(format!(
"Created by whitebox_tools\' {} tool",
self.get_tool_name()
));
output.add_metadata_entry(format!("Input d8 pointer file: {}", d8_file));
output.add_metadata_entry(format!("Input streams file: {}", streams_file));
output.add_metadata_entry(format!("Elapsed Time (excluding I/O): {}", elapsed_time));
if verbose {
println!("Saving data...")
};
let _ = match output.write() {
Ok(_) => {
if verbose {
println!("Output file written")
}
}
Err(e) => return Err(e),
};
if verbose {
println!(
"{}",
&format!("Elapsed Time (excluding I/O): {}", elapsed_time)
);
}
Ok(())
}
}