1 /****************************************************************************
2 *
3 * MODULE: r.in.Lidar
4 *
5 * AUTHOR(S): Markus Metz
6 * Vaclav Petras (base raster addition and refactoring)
7 * Based on r.in.xyz by Hamish Bowman, Volker Wichmann
8 *
9 * PURPOSE: Imports LAS LiDAR point clouds to a raster map using
10 * aggregate statistics.
11 *
12 * COPYRIGHT: (C) 2011-2015 Markus Metz and the The GRASS Development Team
13 *
14 * This program is free software under the GNU General Public
15 * License (>=v2). Read the file COPYING that comes with GRASS
16 * for details.
17 *
18 *****************************************************************************/
19
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <string.h>
23 #include <unistd.h>
24 #include <math.h>
25 #include <sys/types.h>
26 #include <grass/gis.h>
27 #include <grass/raster.h>
28 #include <grass/segment.h>
29 #include <grass/gprojects.h>
30 #include <grass/glocale.h>
31 #include <liblas/capi/liblas.h>
32
33 #include "local_proto.h"
34 #include "rast_segment.h"
35 #include "point_binning.h"
36 #include "filters.h"
37
38
main(int argc,char * argv[])39 int main(int argc, char *argv[])
40 {
41 int out_fd, base_raster;
42 char *infile, *outmap;
43 int percent;
44 double zrange_min, zrange_max, d_tmp;
45 double irange_min, irange_max;
46 unsigned long estimated_lines;
47
48 RASTER_MAP_TYPE rtype, base_raster_data_type;
49 struct History history;
50 char title[64];
51 SEGMENT base_segment;
52 struct PointBinning point_binning;
53 void *base_array;
54 void *raster_row;
55 struct Cell_head region;
56 struct Cell_head input_region;
57 int rows, last_rows, row0, cols; /* scan box size */
58 int row; /* counters */
59
60 int pass, npasses;
61 unsigned long line, line_total;
62 unsigned int counter;
63 unsigned long n_invalid;
64 char buff[BUFFSIZE];
65 double x, y, z;
66 double intensity;
67 int arr_row, arr_col;
68 unsigned long count, count_total;
69 int point_class;
70
71 double zscale = 1.0;
72 double iscale = 1.0;
73 double res = 0.0;
74
75 struct BinIndex bin_index_nodes;
76 bin_index_nodes.num_nodes = 0;
77 bin_index_nodes.max_nodes = 0;
78 bin_index_nodes.nodes = 0;
79
80 struct GModule *module;
81 struct Option *input_opt, *output_opt, *percent_opt, *type_opt, *filter_opt, *class_opt;
82 struct Option *method_opt, *base_raster_opt;
83 struct Option *zrange_opt, *zscale_opt;
84 struct Option *irange_opt, *iscale_opt;
85 struct Option *trim_opt, *pth_opt, *res_opt;
86 struct Option *file_list_opt;
87 struct Flag *print_flag, *scan_flag, *shell_style, *over_flag, *extents_flag;
88 struct Flag *intens_flag, *intens_import_flag;
89 struct Flag *set_region_flag;
90 struct Flag *base_rast_res_flag;
91 struct Flag *only_valid_flag;
92
93 /* LAS */
94 LASReaderH LAS_reader;
95 LASHeaderH LAS_header;
96 LASSRSH LAS_srs;
97 LASPointH LAS_point;
98 int return_filter;
99
100 const char *projstr;
101 struct Cell_head cellhd, loc_wind;
102
103 unsigned int n_filtered;
104
105 G_gisinit(argv[0]);
106
107 module = G_define_module();
108 G_add_keyword(_("raster"));
109 G_add_keyword(_("import"));
110 G_add_keyword(_("LIDAR"));
111 G_add_keyword(_("statistics"));
112 G_add_keyword(_("conversion"));
113 G_add_keyword(_("aggregation"));
114 G_add_keyword(_("binning"));
115 module->description =
116 _("Creates a raster map from LAS LiDAR points using univariate statistics.");
117
118 input_opt = G_define_standard_option(G_OPT_F_BIN_INPUT);
119 input_opt->required = NO;
120 input_opt->label = _("LAS input file");
121 input_opt->description = _("LiDAR input files in LAS format (*.las or *.laz)");
122 input_opt->guisection = _("Input");
123
124 output_opt = G_define_standard_option(G_OPT_R_OUTPUT);
125 output_opt->required = NO;
126 output_opt->guisection = _("Output");
127
128 file_list_opt = G_define_standard_option(G_OPT_F_INPUT);
129 file_list_opt->key = "file";
130 file_list_opt->label = _("File containing names of LAS input files");
131 file_list_opt->description = _("LiDAR input files in LAS format (*.las or *.laz)");
132 file_list_opt->required = NO;
133 file_list_opt->guisection = _("Input");
134
135 method_opt = G_define_option();
136 method_opt->key = "method";
137 method_opt->type = TYPE_STRING;
138 method_opt->required = NO;
139 method_opt->description = _("Statistic to use for raster values");
140 method_opt->options =
141 "n,min,max,range,sum,mean,stddev,variance,coeff_var,median,percentile,skewness,trimmean";
142 method_opt->answer = "mean";
143 method_opt->guisection = _("Statistic");
144 G_asprintf((char **)&(method_opt->descriptions),
145 "n;%s;"
146 "min;%s;"
147 "max;%s;"
148 "range;%s;"
149 "sum;%s;"
150 "mean;%s;"
151 "stddev;%s;"
152 "variance;%s;"
153 "coeff_var;%s;"
154 "median;%s;"
155 "percentile;%s;"
156 "skewness;%s;"
157 "trimmean;%s",
158 _("Number of points in cell"),
159 _("Minimum value of point values in cell"),
160 _("Maximum value of point values in cell"),
161 _("Range of point values in cell"),
162 _("Sum of point values in cell"),
163 _("Mean (average) value of point values in cell"),
164 _("Standard deviation of point values in cell"),
165 _("Variance of point values in cell"),
166 _("Coefficient of variance of point values in cell"),
167 _("Median value of point values in cell"),
168 _("pth (nth) percentile of point values in cell"),
169 _("Skewness of point values in cell"),
170 _("Trimmed mean of point values in cell"));
171
172 type_opt = G_define_standard_option(G_OPT_R_TYPE);
173 type_opt->required = NO;
174 type_opt->answer = "FCELL";
175
176 base_raster_opt = G_define_standard_option(G_OPT_R_INPUT);
177 base_raster_opt->key = "base_raster";
178 base_raster_opt->required = NO;
179 base_raster_opt->label =
180 _("Subtract raster values from the Z coordinates");
181 base_raster_opt->description =
182 _("The scale for Z is applied beforehand, the range filter for"
183 " Z afterwards");
184 base_raster_opt->guisection = _("Transform");
185
186 zrange_opt = G_define_option();
187 zrange_opt->key = "zrange";
188 zrange_opt->type = TYPE_DOUBLE;
189 zrange_opt->required = NO;
190 zrange_opt->key_desc = "min,max";
191 zrange_opt->label = _("Filter range for Z data (min,max)");
192 zrange_opt->description = _("Applied after base_raster transformation step");
193 zrange_opt->guisection = _("Selection");
194
195 zscale_opt = G_define_option();
196 zscale_opt->key = "zscale";
197 zscale_opt->type = TYPE_DOUBLE;
198 zscale_opt->required = NO;
199 zscale_opt->answer = "1.0";
200 zscale_opt->description = _("Scale to apply to Z data");
201 zscale_opt->guisection = _("Transform");
202
203 irange_opt = G_define_option();
204 irange_opt->key = "intensity_range";
205 irange_opt->type = TYPE_DOUBLE;
206 irange_opt->required = NO;
207 irange_opt->key_desc = "min,max";
208 irange_opt->description = _("Filter range for intensity values (min,max)");
209 irange_opt->guisection = _("Selection");
210
211 iscale_opt = G_define_option();
212 iscale_opt->key = "intensity_scale";
213 iscale_opt->type = TYPE_DOUBLE;
214 iscale_opt->required = NO;
215 iscale_opt->answer = "1.0";
216 iscale_opt->description = _("Scale to apply to intensity values");
217 iscale_opt->guisection = _("Transform");
218
219 percent_opt = G_define_option();
220 percent_opt->key = "percent";
221 percent_opt->type = TYPE_INTEGER;
222 percent_opt->required = NO;
223 percent_opt->answer = "100";
224 percent_opt->options = "1-100";
225 percent_opt->description = _("Percent of map to keep in memory");
226
227 /* I would prefer to call the following "percentile", but that has too
228 * much namespace overlap with the "percent" option above */
229 pth_opt = G_define_option();
230 pth_opt->key = "pth";
231 pth_opt->type = TYPE_INTEGER;
232 pth_opt->required = NO;
233 pth_opt->options = "1-100";
234 pth_opt->description = _("pth percentile of the values");
235 pth_opt->guisection = _("Statistic");
236
237 trim_opt = G_define_option();
238 trim_opt->key = "trim";
239 trim_opt->type = TYPE_DOUBLE;
240 trim_opt->required = NO;
241 trim_opt->options = "0-50";
242 trim_opt->label = _("Discard given percentage of the smallest and largest values");
243 trim_opt->description =
244 _("Discard <trim> percent of the smallest and <trim> percent of the largest observations");
245 trim_opt->guisection = _("Statistic");
246
247 res_opt = G_define_option();
248 res_opt->key = "resolution";
249 res_opt->type = TYPE_DOUBLE;
250 res_opt->required = NO;
251 res_opt->description =
252 _("Output raster resolution");
253 res_opt->guisection = _("Output");
254
255 filter_opt = G_define_option();
256 filter_opt->key = "return_filter";
257 filter_opt->type = TYPE_STRING;
258 filter_opt->required = NO;
259 filter_opt->label = _("Only import points of selected return type");
260 filter_opt->description = _("If not specified, all points are imported");
261 filter_opt->options = "first,last,mid";
262 filter_opt->guisection = _("Selection");
263
264 class_opt = G_define_option();
265 class_opt->key = "class_filter";
266 class_opt->type = TYPE_INTEGER;
267 class_opt->multiple = YES;
268 class_opt->required = NO;
269 class_opt->label = _("Only import points of selected class(es)");
270 class_opt->description = _("Input is comma separated integers. "
271 "If not specified, all points are imported.");
272 class_opt->guisection = _("Selection");
273
274 print_flag = G_define_flag();
275 print_flag->key = 'p';
276 print_flag->description =
277 _("Print LAS file info and exit");
278
279 extents_flag = G_define_flag();
280 extents_flag->key = 'e';
281 extents_flag->label =
282 _("Use the extent of the input for the raster extent");
283 extents_flag->description =
284 _("Set internally computational region extents based on the"
285 " point cloud");
286 extents_flag->guisection = _("Output");
287
288 set_region_flag = G_define_flag();
289 set_region_flag->key = 'n';
290 set_region_flag->label =
291 _("Set computation region to match the new raster map");
292 set_region_flag->description =
293 _("Set computation region to match the 2D extent and resolution"
294 " of the newly created new raster map");
295 set_region_flag->guisection = _("Output");
296
297 over_flag = G_define_flag();
298 over_flag->key = 'o';
299 over_flag->label =
300 _("Override projection check (use current location's projection)");
301 over_flag->description =
302 _("Assume that the dataset has same projection as the current location");
303
304 scan_flag = G_define_flag();
305 scan_flag->key = 's';
306 scan_flag->description = _("Scan data file for extent then exit");
307
308 shell_style = G_define_flag();
309 shell_style->key = 'g';
310 shell_style->description =
311 _("In scan mode, print using shell script style");
312
313 intens_flag = G_define_flag();
314 intens_flag->key = 'i';
315 intens_flag->label =
316 _("Use intensity values rather than Z values");
317 intens_flag->description =
318 _("Uses intensity values everywhere as if they would be Z"
319 " coordinates");
320
321 intens_import_flag = G_define_flag();
322 intens_import_flag->key = 'j';
323 intens_import_flag->description =
324 _("Use Z values for filtering, but intensity values for statistics");
325
326 base_rast_res_flag = G_define_flag();
327 base_rast_res_flag->key = 'd';
328 base_rast_res_flag->label =
329 _("Use base raster resolution instead of computational region");
330 base_rast_res_flag->description =
331 _("For getting values from base raster, use its actual"
332 " resolution instead of computational region resolution");
333 base_rast_res_flag->guisection = _("Transform");
334
335 only_valid_flag = G_define_flag();
336 only_valid_flag->key = 'v';
337 only_valid_flag->label = _("Use only valid points");
338 only_valid_flag->description =
339 _("Points invalid according to APSRS LAS specification will be"
340 " filtered out");
341 only_valid_flag->guisection = _("Selection");
342
343 G_option_required(input_opt, file_list_opt, NULL);
344 G_option_exclusive(input_opt, file_list_opt, NULL);
345 G_option_required(output_opt, print_flag, scan_flag, shell_style, NULL);
346 G_option_exclusive(intens_flag, intens_import_flag, NULL);
347 G_option_requires(base_rast_res_flag, base_raster_opt, NULL);
348
349 if (G_parser(argc, argv))
350 exit(EXIT_FAILURE);
351
352 int only_valid = FALSE;
353 n_invalid = 0;
354 if (only_valid_flag->answer)
355 only_valid = TRUE;
356
357 /* we could use rules but this gives more info and allows continuing */
358 if (set_region_flag->answer && !(extents_flag->answer || res_opt->answer)) {
359 G_warning(_("Flag %c makes sense only with %s option or -%c flag"),
360 set_region_flag->key, res_opt->key, extents_flag->key);
361 /* avoid the call later on */
362 set_region_flag->answer = '\0';
363 }
364
365 /* Trim option is used only for trimmean method */
366 if (trim_opt->answer != NULL && strcmp(method_opt->answer, "trimmean") != 0) {
367 G_fatal_error(_("Trim option can be used only with trimmean method"));
368 }
369
370 struct StringList infiles;
371
372 if (file_list_opt->answer) {
373 if (access(file_list_opt->answer, F_OK) != 0)
374 G_fatal_error(_("File <%s> does not exist"), file_list_opt->answer);
375 string_list_from_file(&infiles, file_list_opt->answer);
376 }
377 else {
378 string_list_from_one_item(&infiles, input_opt->answer);
379 }
380
381 /* parse input values */
382 outmap = output_opt->answer;
383
384 if (shell_style->answer && !scan_flag->answer) {
385 scan_flag->answer = 1; /* pointer not int, so set = shell_style->answer ? */
386 }
387
388 /* check zrange and extent relation */
389 if (scan_flag->answer || extents_flag->answer) {
390 if (zrange_opt->answer)
391 G_warning(_("zrange will not be taken into account during scan"));
392 }
393
394 Rast_get_window(®ion);
395 /* G_get_window seems to be unreliable if the location has been changed */
396 G_get_set_window(&loc_wind); /* TODO: v.in.lidar uses G_get_default_window() */
397
398 estimated_lines = 0;
399 int i;
400 for (i = 0; i < infiles.num_items; i++) {
401 infile = infiles.items[i];
402 /* don't if file not found */
403 if (access(infile, F_OK) != 0)
404 G_fatal_error(_("Input file <%s> does not exist"), infile);
405 /* Open LAS file*/
406 LAS_reader = LASReader_Create(infile);
407 if (LAS_reader == NULL)
408 G_fatal_error(_("Unable to open file <%s> as a LiDAR point cloud"),
409 infile);
410 LAS_header = LASReader_GetHeader(LAS_reader);
411 if (LAS_header == NULL) {
412 G_fatal_error(_("Unable to read LAS header of <%s>"), infile);
413 }
414
415 LAS_srs = LASHeader_GetSRS(LAS_header);
416
417 /* print info or check projection if we are actually importing */
418 if (print_flag->answer) {
419 /* print filename when there is more than one file */
420 if (infiles.num_items > 1)
421 fprintf(stdout, "File: %s\n", infile);
422 /* Print LAS header */
423 print_lasinfo(LAS_header, LAS_srs);
424 }
425 else {
426 /* report that we are checking more files */
427 if (i == 1)
428 G_message(_("First file's projection checked,"
429 " checking projection of the other files..."));
430 /* Fetch input map projection in GRASS form. */
431 projstr = LASSRS_GetWKT_CompoundOK(LAS_srs);
432 /* we are printing the non-warning messages only for first file */
433 projection_check_wkt(cellhd, loc_wind, projstr, over_flag->answer,
434 shell_style->answer || i);
435 /* if there is a problem in some other file, first OK message
436 * is printed but than a warning, this is not ideal but hopefully
437 * not so confusing when importing multiple files */
438 }
439 if (scan_flag->answer || extents_flag->answer) {
440 /* we assign to the first one (i==0) but update for the rest */
441 scan_bounds(LAS_reader, shell_style->answer, extents_flag->answer, i,
442 zscale, ®ion);
443 }
444 /* number of estimated point across all files */
445 /* TODO: this should be ull which won't work with percent report */
446 estimated_lines += LASHeader_GetPointRecordsCount(LAS_header);
447 /* We are closing all again and we will be opening them later,
448 * so we don't have to worry about limit for open files. */
449 LASSRS_Destroy(LAS_srs);
450 LASHeader_Destroy(LAS_header);
451 LASReader_Destroy(LAS_reader);
452 }
453 /* if we are not importing, end */
454 if (print_flag->answer || scan_flag->answer)
455 exit(EXIT_SUCCESS);
456
457 return_filter = LAS_ALL;
458 if (filter_opt->answer) {
459 if (strcmp(filter_opt->answer, "first") == 0)
460 return_filter = LAS_FIRST;
461 else if (strcmp(filter_opt->answer, "last") == 0)
462 return_filter = LAS_LAST;
463 else if (strcmp(filter_opt->answer, "mid") == 0)
464 return_filter = LAS_MID;
465 else
466 G_fatal_error(_("Unknown filter option <%s>"), filter_opt->answer);
467 }
468 struct ReturnFilter return_filter_struct;
469 return_filter_struct.filter = return_filter;
470 struct ClassFilter class_filter;
471 class_filter_create_from_strings(&class_filter, class_opt->answers);
472
473 percent = atoi(percent_opt->answer);
474 /* TODO: we already used zscale */
475 /* TODO: we don't report intensity range */
476 if (zscale_opt->answer)
477 zscale = atof(zscale_opt->answer);
478 if (iscale_opt->answer)
479 iscale = atof(iscale_opt->answer);
480
481 /* parse zrange */
482 if (zrange_opt->answer != NULL) {
483 if (zrange_opt->answers[0] == NULL)
484 G_fatal_error(_("Invalid zrange"));
485
486 sscanf(zrange_opt->answers[0], "%lf", &zrange_min);
487 sscanf(zrange_opt->answers[1], "%lf", &zrange_max);
488
489 if (zrange_min > zrange_max) {
490 d_tmp = zrange_max;
491 zrange_max = zrange_min;
492 zrange_min = d_tmp;
493 }
494 }
495 /* parse irange */
496 if (irange_opt->answer != NULL) {
497 if (irange_opt->answers[0] == NULL)
498 G_fatal_error(_("Invalid %s"), irange_opt->key);
499
500 sscanf(irange_opt->answers[0], "%lf", &irange_min);
501 sscanf(irange_opt->answers[1], "%lf", &irange_max);
502
503 if (irange_min > irange_max) {
504 d_tmp = irange_max;
505 irange_max = irange_min;
506 irange_min = d_tmp;
507 }
508 }
509
510 point_binning_set(&point_binning, method_opt->answer, pth_opt->answer,
511 trim_opt->answer, FALSE);
512
513 base_array = NULL;
514
515 if (strcmp("CELL", type_opt->answer) == 0)
516 rtype = CELL_TYPE;
517 else if (strcmp("DCELL", type_opt->answer) == 0)
518 rtype = DCELL_TYPE;
519 else
520 rtype = FCELL_TYPE;
521
522 if (point_binning.method == METHOD_N)
523 rtype = CELL_TYPE;
524
525 if (res_opt->answer) {
526 /* align to resolution */
527 res = atof(res_opt->answer);
528
529 if (!G_scan_resolution(res_opt->answer, &res, region.proj))
530 G_fatal_error(_("Invalid input <%s=%s>"), res_opt->key, res_opt->answer);
531
532 if (res <= 0)
533 G_fatal_error(_("Option '%s' must be > 0.0"), res_opt->key);
534
535 region.ns_res = region.ew_res = res;
536
537 region.north = ceil(region.north / res) * res;
538 region.south = floor(region.south / res) * res;
539 region.east = ceil(region.east / res) * res;
540 region.west = floor(region.west / res) * res;
541
542 G_adjust_Cell_head(®ion, 0, 0);
543 }
544 else if (extents_flag->answer) {
545 /* align to current region */
546 Rast_align_window(®ion, &loc_wind);
547 }
548 Rast_set_output_window(®ion);
549
550 rows = last_rows = region.rows;
551 npasses = 1;
552 if (percent < 100) {
553 rows = (int)(region.rows * (percent / 100.0));
554 npasses = region.rows / rows;
555 last_rows = region.rows - npasses * rows;
556 if (last_rows)
557 npasses++;
558 else
559 last_rows = rows;
560
561 }
562 cols = region.cols;
563
564 G_debug(2, "region.n=%f region.s=%f region.ns_res=%f", region.north,
565 region.south, region.ns_res);
566 G_debug(2, "region.rows=%d [box_rows=%d] region.cols=%d", region.rows,
567 rows, region.cols);
568
569 /* using row-based chunks (used for output) when input and output
570 * region matches and using segment library when they don't */
571 int use_segment = 0;
572 int use_base_raster_res = 0;
573 /* TODO: see if the input region extent is smaller than the raster
574 * if yes, the we need to load the whole base raster if the -e
575 * flag was defined (alternatively clip the regions) */
576 if (base_rast_res_flag->answer)
577 use_base_raster_res = 1;
578 if (base_raster_opt->answer && (res_opt->answer || use_base_raster_res
579 || extents_flag->answer))
580 use_segment = 1;
581 if (base_raster_opt->answer && !use_segment) {
582 /* TODO: do we need to test existence first? mapset? */
583 base_raster = Rast_open_old(base_raster_opt->answer, "");
584 base_raster_data_type = Rast_get_map_type(base_raster);
585 base_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(base_raster_data_type));
586 }
587 if (base_raster_opt->answer && use_segment) {
588 if (use_base_raster_res) {
589 /* read raster actual extent and resolution */
590 Rast_get_cellhd(base_raster_opt->answer, "", &input_region);
591 /* TODO: make it only as small as the output is or points are */
592 Rast_set_input_window(&input_region); /* we have split window */
593 } else {
594 Rast_get_input_window(&input_region);
595 }
596 rast_segment_open(&base_segment, base_raster_opt->answer, &base_raster_data_type);
597 }
598
599 if (!scan_flag->answer) {
600 if (!check_rows_cols_fit_to_size_t(rows, cols))
601 G_fatal_error(_("Unable to process the whole map at once. "
602 "Please set the '%s' option to some value lower than 100."),
603 percent_opt->key);
604 point_binning_memory_test(&point_binning, rows, cols, rtype);
605 }
606
607 /* open output map */
608 out_fd = Rast_open_new(outmap, rtype);
609
610 /* allocate memory for a single row of output data */
611 raster_row = Rast_allocate_output_buf(rtype);
612
613 G_message(_("Reading data..."));
614
615 count_total = line_total = 0;
616
617 /* main binning loop(s) */
618 for (pass = 1; pass <= npasses; pass++) {
619
620 if (npasses > 1)
621 G_message(_("Pass #%d (of %d)..."), pass, npasses);
622
623 /* figure out segmentation */
624 row0 = (pass - 1) * rows;
625 if (pass == npasses) {
626 rows = last_rows;
627 }
628
629 if (base_array) {
630 G_debug(2, "filling base raster array");
631 for (row = 0; row < rows; row++) {
632 Rast_get_row(base_raster, base_array + ((size_t) row * cols * Rast_cell_size(base_raster_data_type)), row, base_raster_data_type);
633 }
634 }
635
636 G_debug(2, "pass=%d/%d rows=%d", pass, npasses, rows);
637
638 point_binning_allocate(&point_binning, rows, cols, rtype);
639
640 line = 0;
641 count = 0;
642 counter = 0;
643 G_percent_reset();
644
645 /* loop of input files */
646 for (i = 0; i < infiles.num_items; i++) {
647 infile = infiles.items[i];
648 /* we already know file is there, so just do basic checks */
649 LAS_reader = LASReader_Create(infile);
650 if (LAS_reader == NULL)
651 G_fatal_error(_("Unable to open file <%s>"), infile);
652
653 while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
654 line++;
655 counter++;
656
657 if (counter == 100000) { /* speed */
658 if (line < estimated_lines)
659 G_percent(line, estimated_lines, 3);
660 counter = 0;
661 }
662
663 /* We always count them and report because behavior
664 * changed in between 7.0 and 7.2 from undefined (but skipping
665 * invalid points) to filtering them out only when requested. */
666 if (!LASPoint_IsValid(LAS_point)) {
667 n_invalid++;
668 if (only_valid)
669 continue;
670 }
671
672 x = LASPoint_GetX(LAS_point);
673 y = LASPoint_GetY(LAS_point);
674 if (intens_flag->answer)
675 /* use intensity as z here to allow all filters (and
676 * modifications) below to be applied for intensity */
677 z = LASPoint_GetIntensity(LAS_point);
678 else
679 z = LASPoint_GetZ(LAS_point);
680
681 int return_n = LASPoint_GetReturnNumber(LAS_point);
682 int n_returns = LASPoint_GetNumberOfReturns(LAS_point);
683 if (return_filter_is_out(&return_filter_struct, return_n, n_returns)) {
684 n_filtered++;
685 continue;
686 }
687 point_class = (int) LASPoint_GetClassification(LAS_point);
688 if (class_filter_is_out(&class_filter, point_class))
689 continue;
690
691 if (y <= region.south || y > region.north) {
692 continue;
693 }
694 if (x < region.west || x >= region.east) {
695 continue;
696 }
697
698 /* find the bin in the current array box */
699 arr_row = (int)((region.north - y) / region.ns_res) - row0;
700 if (arr_row < 0 || arr_row >= rows)
701 continue;
702 arr_col = (int)((x - region.west) / region.ew_res);
703
704 z = z * zscale;
705
706 if (base_array) {
707 double base_z;
708 if (row_array_get_value_row_col(base_array, arr_row, arr_col,
709 cols, base_raster_data_type,
710 &base_z))
711 z -= base_z;
712 else
713 continue;
714 }
715 else if (use_segment) {
716 double base_z;
717 if (rast_segment_get_value_xy(&base_segment, &input_region,
718 base_raster_data_type, x, y,
719 &base_z))
720 z -= base_z;
721 else
722 continue;
723 }
724
725 if (zrange_opt->answer) {
726 if (z < zrange_min || z > zrange_max) {
727 continue;
728 }
729 }
730
731 if (intens_import_flag->answer || irange_opt->answer) {
732 intensity = LASPoint_GetIntensity(LAS_point);
733 intensity *= iscale;
734 if (irange_opt->answer) {
735 if (intensity < irange_min || intensity > irange_max) {
736 continue;
737 }
738 }
739 /* use intensity for statistics */
740 if (intens_import_flag->answer)
741 z = intensity;
742 }
743
744 count++;
745 /* G_debug(5, "x: %f, y: %f, z: %f", x, y, z); */
746
747 update_value(&point_binning, &bin_index_nodes, cols,
748 arr_row, arr_col, rtype, x, y, z);
749 } /* while !EOF of one input file */
750 /* close input LAS file */
751 LASReader_Destroy(LAS_reader);
752 } /* end of loop for all input files files */
753
754 G_percent(1, 1, 1); /* flush */
755 G_debug(2, "pass %d finished, %lu coordinates in box", pass, count);
756 count_total += count;
757 line_total += line;
758
759 /* calc stats and output */
760 G_message(_("Writing output raster map..."));
761 for (row = 0; row < rows; row++) {
762 /* potentially vector writing can be independent on the binning */
763 write_values(&point_binning, &bin_index_nodes, raster_row, row,
764 cols, rtype, NULL);
765
766 G_percent(row, rows, 10);
767
768 /* write out line of raster data */
769 Rast_put_row(out_fd, raster_row, rtype);
770 }
771
772 /* free memory */
773 point_binning_free(&point_binning, &bin_index_nodes);
774 } /* passes loop */
775 if (base_array)
776 Rast_close(base_raster);
777 if (use_segment)
778 Segment_close(&base_segment);
779
780 G_percent(1, 1, 1); /* flush */
781 G_free(raster_row);
782
783 G_message(_("%lu points found in input file(s)"), line_total);
784
785 /* close raster file & write history */
786 Rast_close(out_fd);
787
788 sprintf(title, "Raw X,Y,Z data binned into a raster grid by cell %s",
789 method_opt->answer);
790 Rast_put_cell_title(outmap, title);
791
792 Rast_short_history(outmap, "raster", &history);
793 Rast_command_history(&history);
794 Rast_set_history(&history, HIST_DATSRC_1, infile);
795 Rast_write_history(outmap, &history);
796
797 /* set computation region to the new raster map */
798 /* TODO: should be in the done message */
799 if (set_region_flag->answer)
800 G_put_window(®ion);
801
802 if (n_invalid && only_valid)
803 G_message(_("%lu input points were invalid and filtered out"),
804 n_invalid);
805 if (n_invalid && !only_valid)
806 G_message(_("%lu input points were invalid, use -%c flag to filter"
807 " them out"), n_invalid, only_valid_flag->key);
808 if (infiles.num_items > 1) {
809 sprintf(buff, _("Raster map <%s> created."
810 " %lu points from %d files found in region."),
811 outmap, count_total, infiles.num_items);
812 }
813 else {
814 sprintf(buff, _("Raster map <%s> created."
815 " %lu points found in region."),
816 outmap, count_total);
817 }
818
819 G_done_msg("%s", buff);
820 G_debug(1, "Processed %lu points.", line_total);
821
822 string_list_free(&infiles);
823
824 exit(EXIT_SUCCESS);
825
826 }
827