1
2 /****************************************************************
3 *
4 * MODULE: v.in.lidar
5 *
6 * AUTHOR(S): Markus Metz
7 * Vaclav Petras (decimation, cats, areas, zrange)
8 * based on v.in.ogr
9 *
10 * PURPOSE: Import LiDAR LAS points
11 *
12 * COPYRIGHT: (C) 2011-2015 by the GRASS Development Team
13 *
14 * This program is free software under the
15 * GNU General Public License (>=v2).
16 * Read the file COPYING that comes with GRASS
17 * for details.
18 **************************************************************/
19
20 #include <stdlib.h>
21 #include <string.h>
22 #include <unistd.h>
23 #include <grass/gis.h>
24 #include <grass/dbmi.h>
25 #include <grass/vector.h>
26 #include <grass/gprojects.h>
27 #include <grass/glocale.h>
28 #include <liblas/capi/liblas.h>
29
30 #include "count_decimation.h"
31 #include "projection.h"
32 #include "lidar.h"
33 #include "attributes.h"
34 #include "info.h"
35 #include "vector_mask.h"
36 #include "filters.h"
37
38 #ifndef MAX
39 # define MIN(a,b) ((a<b) ? a : b)
40 # define MAX(a,b) ((a>b) ? a : b)
41 #endif
42
check_layers_not_equal(int primary,int secondary,const char * primary_name,const char * secondary_name)43 static void check_layers_not_equal(int primary, int secondary,
44 const char *primary_name,
45 const char *secondary_name)
46 {
47 if (primary && primary == secondary)
48 G_fatal_error(_("Values of %s and %s are the same."
49 " All categories would be stored only"
50 " in layer number <%d>"), primary_name,
51 secondary_name, primary);
52 }
53
check_layers_in_list_not_equal(struct Option ** options,int * values,size_t size)54 static void check_layers_in_list_not_equal(struct Option **options,
55 int *values, size_t size)
56 {
57 size_t layer_index_1, layer_index_2;
58 for (layer_index_1 = 0; layer_index_1 < size; layer_index_1++) {
59 for (layer_index_2 = 0; layer_index_2 < size; layer_index_2++) {
60 if (layer_index_1 != layer_index_2) {
61 check_layers_not_equal(values[layer_index_1],
62 values[layer_index_2],
63 options[layer_index_1]->key,
64 options[layer_index_2]->key);
65 }
66 }
67 }
68 }
69
70
main(int argc,char * argv[])71 int main(int argc, char *argv[])
72 {
73 int i;
74 float xmin = 0., ymin = 0., xmax = 0., ymax = 0.;
75 struct GModule *module;
76 struct Option *in_opt, *out_opt, *spat_opt, *filter_opt, *class_opt;
77 struct Option *id_layer_opt;
78 struct Option *return_layer_opt;
79 struct Option *class_layer_opt;
80 struct Option *rgb_layer_opt;
81 struct Option *vector_mask_opt, *vector_mask_field_opt;
82 struct Option *skip_opt, *preserve_opt, *offset_opt, *limit_opt;
83 struct Option *outloc_opt, *zrange_opt;
84 struct Flag *print_flag, *notab_flag, *region_flag, *notopo_flag;
85 struct Flag *nocats_flag;
86 struct Flag *over_flag, *extend_flag, *no_import_flag;
87 struct Flag *invert_mask_flag;
88 struct Flag *only_valid_flag;
89 char buf[2000];
90 struct Key_Value *loc_proj_info = NULL, *loc_proj_units = NULL;
91 struct Key_Value *proj_info, *proj_units;
92 const char *projstr;
93 struct Cell_head cellhd, loc_wind, cur_wind;
94 char error_msg[8192];
95
96 /* Vector */
97 struct Map_info Map;
98 int cat;
99
100 /* Attributes */
101 struct field_info *Fi;
102 dbDriver *driver;
103
104 /* LAS */
105 LASReaderH LAS_reader;
106 LASHeaderH LAS_header;
107 LASSRSH LAS_srs;
108 LASPointH LAS_point;
109 double scale_x, scale_y, scale_z, offset_x, offset_y, offset_z;
110 int las_point_format;
111 int have_time, have_color;
112 int point_class;
113
114 struct line_pnts *Points;
115 struct line_cats *Cats;
116
117 int cat_max_reached = FALSE;
118
119 #ifdef HAVE_LONG_LONG_INT
120 unsigned long long n_features; /* what libLAS reports as point count */
121 unsigned long long points_imported; /* counter how much we have imported */
122 unsigned long long feature_count, n_outside, zrange_filtered,
123 n_outside_mask, n_filtered, n_class_filtered, not_valid;
124 #else
125 unsigned long n_features;
126 unsigned long points_imported;
127 unsigned long feature_count, n_outside, zrange_filtered,
128 n_outside_mask, n_filtered, n_class_filtered, not_valid;
129 #endif
130
131 int overwrite;
132
133 G_gisinit(argv[0]);
134
135 module = G_define_module();
136 G_add_keyword(_("vector"));
137 G_add_keyword(_("import"));
138 G_add_keyword(_("LIDAR"));
139 G_add_keyword(_("level1"));
140
141 module->description = _("Converts LAS LiDAR point clouds to a GRASS vector map with libLAS.");
142
143 in_opt = G_define_standard_option(G_OPT_F_BIN_INPUT);
144 in_opt->label = _("LAS input file");
145 in_opt->description = _("LiDAR input files in LAS format (*.las or *.laz)");
146
147 out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
148
149 id_layer_opt = G_define_standard_option(G_OPT_V_FIELD);
150 id_layer_opt->key = "id_layer";
151 id_layer_opt->label = _("Layer number to store generated point ID as category");
152 id_layer_opt->description = _("Set to 1 by default, use -c to not store it");
153 id_layer_opt->answer = NULL;
154 id_layer_opt->guisection = _("Categories");
155
156 return_layer_opt = G_define_standard_option(G_OPT_V_FIELD);
157 return_layer_opt->key = "return_layer";
158 return_layer_opt->label =
159 _("Layer number to store return information as category");
160 return_layer_opt->description = _("Leave empty to not store it");
161 return_layer_opt->answer = NULL;
162 return_layer_opt->guisection = _("Categories");
163
164 class_layer_opt = G_define_standard_option(G_OPT_V_FIELD);
165 class_layer_opt->key = "class_layer";
166 class_layer_opt->label =
167 _("Layer number to store class number as category");
168 class_layer_opt->description = _("Leave empty to not store it");
169 class_layer_opt->answer = NULL;
170 class_layer_opt->guisection = _("Categories");
171
172 rgb_layer_opt = G_define_standard_option(G_OPT_V_FIELD);
173 rgb_layer_opt->key = "rgb_layer";
174 rgb_layer_opt->label =
175 _("Layer number where RBG colors are stored as category");
176 rgb_layer_opt->description = _("Leave empty to not store it");
177 rgb_layer_opt->answer = NULL;
178 rgb_layer_opt->guisection = _("Categories");
179
180 spat_opt = G_define_option();
181 spat_opt->key = "spatial";
182 spat_opt->type = TYPE_DOUBLE;
183 spat_opt->multiple = YES;
184 spat_opt->required = NO;
185 spat_opt->key_desc = "xmin,ymin,xmax,ymax";
186 spat_opt->label = _("Import subregion only");
187 spat_opt->guisection = _("Selection");
188 spat_opt->description =
189 _("Format: xmin,ymin,xmax,ymax - usually W,S,E,N");
190
191 zrange_opt = G_define_option();
192 zrange_opt->key = "zrange";
193 zrange_opt->type = TYPE_DOUBLE;
194 zrange_opt->required = NO;
195 zrange_opt->key_desc = "min,max";
196 zrange_opt->description = _("Filter range for z data (min,max)");
197 zrange_opt->guisection = _("Selection");
198
199 filter_opt = G_define_option();
200 filter_opt->key = "return_filter";
201 filter_opt->type = TYPE_STRING;
202 filter_opt->required = NO;
203 filter_opt->label = _("Only import points of selected return type");
204 filter_opt->description = _("If not specified, all points are imported");
205 filter_opt->options = "first,last,mid";
206 filter_opt->guisection = _("Selection");
207
208 class_opt = G_define_option();
209 class_opt->key = "class_filter";
210 class_opt->type = TYPE_INTEGER;
211 class_opt->multiple = YES;
212 class_opt->required = NO;
213 class_opt->label = _("Only import points of selected class(es)");
214 class_opt->description = _("Input is comma separated integers. "
215 "If not specified, all points are imported.");
216 class_opt->guisection = _("Selection");
217
218 vector_mask_opt = G_define_standard_option(G_OPT_V_INPUT);
219 vector_mask_opt->key = "mask";
220 vector_mask_opt->required = NO;
221 vector_mask_opt->label = _("Areas where to import points");
222 vector_mask_opt->description = _("Name of vector map with areas where the points should be imported");
223 vector_mask_opt->guisection = _("Selection");
224
225 vector_mask_field_opt = G_define_standard_option(G_OPT_V_FIELD);
226 vector_mask_field_opt->key = "mask_layer";
227 vector_mask_field_opt->label = _("Layer number or name for mask option");
228 vector_mask_field_opt->guisection = _("Selection");
229
230 skip_opt = G_define_option();
231 skip_opt->key = "skip";
232 skip_opt->type = TYPE_INTEGER;
233 skip_opt->multiple = NO;
234 skip_opt->required = NO;
235 skip_opt->label = _("Do not import every n-th point");
236 skip_opt->description = _("For example, 5 will import 80 percent of points. "
237 "If not specified, all points are imported");
238 skip_opt->guisection = _("Decimation");
239
240 preserve_opt = G_define_option();
241 preserve_opt->key = "preserve";
242 preserve_opt->type = TYPE_INTEGER;
243 preserve_opt->multiple = NO;
244 preserve_opt->required = NO;
245 preserve_opt->label = _("Import only every n-th point");
246 preserve_opt->description = _("For example, 4 will import 25 percent of points. "
247 "If not specified, all points are imported");
248 preserve_opt->guisection = _("Decimation");
249
250 offset_opt = G_define_option();
251 offset_opt->key = "offset";
252 offset_opt->type = TYPE_INTEGER;
253 offset_opt->multiple = NO;
254 offset_opt->required = NO;
255 offset_opt->label = _("Skip first n points");
256 offset_opt->description = _("Skips the given number of points at the beginning.");
257 offset_opt->guisection = _("Decimation");
258
259 limit_opt = G_define_option();
260 limit_opt->key = "limit";
261 limit_opt->type = TYPE_INTEGER;
262 limit_opt->multiple = NO;
263 limit_opt->required = NO;
264 limit_opt->label = _("Import only n points");
265 limit_opt->description = _("Imports only the given number of points");
266 limit_opt->guisection = _("Decimation");
267
268 outloc_opt = G_define_option();
269 outloc_opt->key = "location";
270 outloc_opt->type = TYPE_STRING;
271 outloc_opt->required = NO;
272 outloc_opt->description = _("Name for new location to create");
273 outloc_opt->key_desc = "name";
274
275 print_flag = G_define_flag();
276 print_flag->key = 'p';
277 print_flag->description =
278 _("Print LAS file info and exit");
279 print_flag->suppress_required = YES;
280
281 region_flag = G_define_flag();
282 region_flag->key = 'r';
283 region_flag->guisection = _("Selection");
284 region_flag->description = _("Limit import to the current region");
285
286 invert_mask_flag = G_define_flag();
287 invert_mask_flag->key = 'u';
288 invert_mask_flag->description = _("Invert mask when selecting points");
289 invert_mask_flag->guisection = _("Selection");
290
291 only_valid_flag = G_define_flag();
292 only_valid_flag->key = 'v';
293 only_valid_flag->label = _("Use only valid points");
294 only_valid_flag->description =
295 _("Points invalid according to APSRS LAS specification will be"
296 " filtered out");
297 only_valid_flag->guisection = _("Selection");
298
299 extend_flag = G_define_flag();
300 extend_flag->key = 'e';
301 extend_flag->description =
302 _("Extend region extents based on new dataset");
303
304 notab_flag = G_define_standard_flag(G_FLG_V_TABLE);
305 notab_flag->guisection = _("Speed");
306
307 nocats_flag = G_define_flag();
308 nocats_flag->key = 'c';
309 nocats_flag->label =
310 _("Do not automatically add unique ID as category to each point");
311 nocats_flag->description =
312 _("Create only requested layers and categories");
313 nocats_flag->guisection = _("Speed");
314
315 notopo_flag = G_define_standard_flag(G_FLG_V_TOPO);
316 notopo_flag->guisection = _("Speed");
317
318 over_flag = G_define_flag();
319 over_flag->key = 'o';
320 over_flag->label =
321 _("Override projection check (use current location's projection)");
322 over_flag->description =
323 _("Assume that the dataset has same projection as the current location");
324
325 no_import_flag = G_define_flag();
326 no_import_flag->key = 'i';
327 no_import_flag->description =
328 _("Create the location specified by the \"location\" parameter and exit."
329 " Do not import the vector data.");
330 no_import_flag->suppress_required = YES;
331
332 G_option_exclusive(skip_opt, preserve_opt, NULL);
333 G_option_requires(nocats_flag, notab_flag, NULL);
334 G_option_exclusive(nocats_flag, id_layer_opt, NULL);
335 G_option_requires(return_layer_opt, id_layer_opt, nocats_flag, NULL);
336 G_option_requires(class_layer_opt, id_layer_opt, nocats_flag, NULL);
337 G_option_requires(rgb_layer_opt, id_layer_opt, nocats_flag, NULL);
338
339 /* The parser checks if the map already exists in current mapset, this is
340 * wrong if location options is used, so we switch out the check and do it
341 * in the module after the parser */
342 overwrite = G_check_overwrite(argc, argv);
343
344 if (G_parser(argc, argv))
345 exit(EXIT_FAILURE);
346
347 /* Don't crash on cmd line if file not found */
348 if (access(in_opt->answer, F_OK) != 0) {
349 G_fatal_error(_("Input file <%s> does not exist"), in_opt->answer);
350 }
351 /* Open LAS file*/
352 LAS_reader = LASReader_Create(in_opt->answer);
353 if (LAS_reader == NULL)
354 G_fatal_error(_("Unable to open file <%s> as a LiDAR point cloud. %s"),
355 in_opt->answer, LASError_GetLastErrorMsg());
356 LAS_header = LASReader_GetHeader(LAS_reader);
357 if (LAS_header == NULL) {
358 G_fatal_error(_("Unable to read LAS header of <%s>"), in_opt->answer);
359 }
360
361 LAS_srs = LASHeader_GetSRS(LAS_header);
362
363 scale_x = LASHeader_GetScaleX(LAS_header);
364 scale_y = LASHeader_GetScaleY(LAS_header);
365 scale_z = LASHeader_GetScaleZ(LAS_header);
366
367 offset_x = LASHeader_GetOffsetX(LAS_header);
368 offset_y = LASHeader_GetOffsetY(LAS_header);
369 offset_z = LASHeader_GetOffsetZ(LAS_header);
370
371 xmin = LASHeader_GetMinX(LAS_header);
372 xmax = LASHeader_GetMaxX(LAS_header);
373 ymin = LASHeader_GetMinY(LAS_header);
374 ymax = LASHeader_GetMaxY(LAS_header);
375
376 /* Print LAS header */
377 if (print_flag->answer) {
378 /* print... */
379 print_lasinfo(LAS_header, LAS_srs);
380
381 LASSRS_Destroy(LAS_srs);
382 LASHeader_Destroy(LAS_header);
383 LASReader_Destroy(LAS_reader);
384
385 exit(EXIT_SUCCESS);
386 }
387
388 int only_valid = FALSE;
389 if (only_valid_flag->answer)
390 only_valid = TRUE;
391
392 struct ReturnFilter return_filter_struct;
393 return_filter_create_from_string(&return_filter_struct, filter_opt->answer);
394 struct ClassFilter class_filter;
395 class_filter_create_from_strings(&class_filter, class_opt->answers);
396
397 int id_layer = 1;
398 int return_layer = 0;
399 int class_layer = 0;
400 int rgb_layer = 0;
401 if (id_layer_opt->answer)
402 id_layer = atoi(id_layer_opt->answer);
403 if (return_layer_opt->answer)
404 return_layer = atoi(return_layer_opt->answer);
405 if (class_layer_opt->answer)
406 class_layer = atoi(class_layer_opt->answer);
407 if (rgb_layer_opt->answer)
408 rgb_layer = atoi(rgb_layer_opt->answer);
409
410 if (nocats_flag->answer) {
411 id_layer = 0;
412 }
413 /* no cats forces no table earlier */
414 if (!notab_flag->answer && !id_layer) {
415 G_message(_("-%c flag is not set but ID layer is not specified"), notab_flag->key);
416 G_fatal_error(_("ID layer is required to store attribute table"));
417 }
418
419 struct Option *layer_options[4] = {id_layer_opt, return_layer_opt,
420 class_layer_opt, rgb_layer_opt};
421 int layer_values[4] = {id_layer, return_layer, class_layer, rgb_layer};
422 check_layers_in_list_not_equal(layer_options, layer_values, 4);
423
424 if (id_layer)
425 G_verbose_message(_("Storing generated point IDs as categories"
426 " in the layer <%d>, consequently no more"
427 " than %d points can be imported"),
428 id_layer, GV_CAT_MAX);
429
430 double zrange_min, zrange_max;
431 int use_zrange = FALSE;
432
433 if (zrange_opt->answer != NULL) {
434 if (zrange_opt->answers[0] == NULL || zrange_opt->answers[1] == NULL)
435 G_fatal_error(_("Invalid zrange <%s>"), zrange_opt->answer);
436 sscanf(zrange_opt->answers[0], "%lf", &zrange_min);
437 sscanf(zrange_opt->answers[1], "%lf", &zrange_max);
438 /* for convenience, switch order to make valid input */
439 if (zrange_min > zrange_max) {
440 double tmp = zrange_max;
441
442 zrange_max = zrange_min;
443 zrange_min = tmp;
444 }
445 use_zrange = TRUE;
446 }
447
448 if (region_flag->answer) {
449 if (spat_opt->answer)
450 G_fatal_error(_("Select either the current region flag or the spatial option, not both"));
451
452 G_get_window(&cur_wind);
453 xmin = cur_wind.west;
454 xmax = cur_wind.east;
455 ymin = cur_wind.south;
456 ymax = cur_wind.north;
457 }
458 if (spat_opt->answer) {
459 /* See as reference: gdal/ogr/ogr_capi_test.c */
460
461 /* cut out a piece of the map */
462 /* order: xmin,ymin,xmax,ymax */
463 int arg_s_num = 0;
464 i = 0;
465 while (spat_opt->answers[i]) {
466 if (i == 0)
467 xmin = atof(spat_opt->answers[i]);
468 if (i == 1)
469 ymin = atof(spat_opt->answers[i]);
470 if (i == 2)
471 xmax = atof(spat_opt->answers[i]);
472 if (i == 3)
473 ymax = atof(spat_opt->answers[i]);
474 arg_s_num++;
475 i++;
476 }
477 if (arg_s_num != 4)
478 G_fatal_error(_("4 parameters required for 'spatial' parameter"));
479 }
480 if (spat_opt->answer || region_flag->answer) {
481 G_debug(2, "cut out with boundaries: xmin:%f ymin:%f xmax:%f ymax:%f",
482 xmin, ymin, xmax, ymax);
483 }
484
485 /* fetch boundaries */
486 G_get_window(&cellhd);
487 cellhd.north = ymax;
488 cellhd.south = ymin;
489 cellhd.west = xmin;
490 cellhd.east = xmax;
491 cellhd.rows = 20; /* TODO - calculate useful values */
492 cellhd.cols = 20;
493 cellhd.ns_res = (cellhd.north - cellhd.south) / cellhd.rows;
494 cellhd.ew_res = (cellhd.east - cellhd.west) / cellhd.cols;
495
496 /* Fetch input map projection in GRASS form. */
497 proj_info = NULL;
498 proj_units = NULL;
499 projstr = LASSRS_GetWKT_CompoundOK(LAS_srs);
500
501 /* Do we need to create a new location? */
502 if (outloc_opt->answer != NULL) {
503 /* Convert projection information non-interactively as we can't
504 * assume the user has a terminal open */
505 if (GPJ_wkt_to_grass(&cellhd, &proj_info,
506 &proj_units, projstr, 0) < 0) {
507 G_fatal_error(_("Unable to convert input map projection to GRASS "
508 "format; cannot create new location."));
509 }
510 else {
511 if (0 != G_make_location(outloc_opt->answer, &cellhd,
512 proj_info, proj_units)) {
513 G_fatal_error(_("Unable to create new location <%s>"),
514 outloc_opt->answer);
515 }
516 G_message(_("Location <%s> created"), outloc_opt->answer);
517 }
518
519 /* If the i flag is set, clean up? and exit here */
520 if(no_import_flag->answer)
521 exit(EXIT_SUCCESS);
522
523 /* TODO: */
524 G_warning("Import into new location not yet implemented");
525 /* at this point the module should be using G_create_alt_env()
526 to change context to the newly created location; once done
527 it should switch back with G_switch_env(). See r.in.gdal */
528 }
529 else {
530 /* Does the projection of the current location match the dataset? */
531 /* G_get_window seems to be unreliable if the location has been changed */
532 G_get_default_window(&loc_wind);
533 projstr = LASSRS_GetWKT_CompoundOK(LAS_srs);
534 /* we are printing the non-warning messages only for first file */
535 projection_check_wkt(cellhd, loc_wind, projstr, over_flag->answer,
536 TRUE);
537 }
538
539 if (!outloc_opt->answer) { /* Check if the map exists */
540 if (G_find_vector2(out_opt->answer, G_mapset())) {
541 if (overwrite)
542 G_warning(_("Vector map <%s> already exists and will be overwritten"),
543 out_opt->answer);
544 else
545 G_fatal_error(_("Vector map <%s> already exists"),
546 out_opt->answer);
547 }
548 }
549
550 /* open output vector */
551 sprintf(buf, "%s", out_opt->answer);
552 /* strip any @mapset from vector output name */
553 G_find_vector(buf, G_mapset());
554 if (Vect_open_new(&Map, out_opt->answer, 1) < 0)
555 G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
556
557 Vect_hist_command(&Map);
558
559 /* libLAS uses uint32_t according to the source code
560 * or unsigned int according to the online doc,
561 * so just storing in long doesn't help.
562 * Thus, we use this just for the messages and percents.
563 */
564 n_features = LASHeader_GetPointRecordsCount(LAS_header);
565 las_point_format = LASHeader_GetDataFormatId(LAS_header);
566
567 have_time = (las_point_format == 1 || las_point_format == 3 ||
568 las_point_format == 4 || las_point_format == 5);
569
570 have_color = (las_point_format == 2 || las_point_format == 3 ||
571 las_point_format == 5);
572
573 /* Add DB link */
574 if (!notab_flag->answer) {
575 create_table_for_lidar(&Map, out_opt->answer, id_layer, &driver,
576 &Fi, have_time, have_color);
577 }
578
579 struct VectorMask vector_mask;
580 if (vector_mask_opt->answer) {
581 VectorMask_init(&vector_mask, vector_mask_opt->answer,
582 vector_mask_field_opt->answer, (int)invert_mask_flag->answer);
583 }
584
585 /* Import feature */
586 points_imported = 0;
587 cat = 1;
588 not_valid = 0;
589 feature_count = 0;
590 n_outside = 0;
591 n_filtered = 0;
592 n_class_filtered = 0;
593 n_outside_mask = 0;
594 zrange_filtered = 0;
595
596 Points = Vect_new_line_struct();
597 Cats = Vect_new_cats_struct();
598
599 struct CountDecimationControl count_decimation_control;
600
601 count_decimation_init_from_str(&count_decimation_control,
602 skip_opt->answer, preserve_opt->answer,
603 offset_opt->answer, limit_opt->answer);
604 if (!count_decimation_is_valid(&count_decimation_control))
605 G_fatal_error(_("Settings for count-based decimation are not valid"));
606 /* we don't check if the decimation is noop */
607
608 #ifdef HAVE_LONG_LONG_INT
609 G_important_message(_("Scanning %llu points..."), n_features);
610 #else
611 G_important_message(_("Scanning %lu points..."), n_features);
612 #endif
613 while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
614 double x, y, z;
615
616 G_percent(feature_count++, n_features, 1); /* show something happens */
617
618 /* We always count them and report because r.in.lidar behavior
619 * changed in between 7.0 and 7.2 from undefined (but skipping
620 * invalid points) to filtering them out only when requested. */
621 if (!LASPoint_IsValid(LAS_point)) {
622 not_valid++;
623 if (only_valid)
624 continue;
625 }
626
627 Vect_reset_line(Points);
628 Vect_reset_cats(Cats);
629
630 x = LASPoint_GetX(LAS_point);
631 y = LASPoint_GetY(LAS_point);
632 z = LASPoint_GetZ(LAS_point);
633
634 if (spat_opt->answer || region_flag->answer) {
635 if (x < xmin || x > xmax || y < ymin || y > ymax) {
636 n_outside++;
637 continue;
638 }
639 }
640 if (use_zrange) {
641 if (z < zrange_min || z > zrange_max) {
642 zrange_filtered++;
643 continue;
644 }
645 }
646 int return_n = LASPoint_GetReturnNumber(LAS_point);
647 int n_returns = LASPoint_GetNumberOfReturns(LAS_point);
648 if (return_filter_is_out(&return_filter_struct, return_n, n_returns)) {
649 n_filtered++;
650 continue;
651 }
652 point_class = (int) LASPoint_GetClassification(LAS_point);
653 if (class_filter_is_out(&class_filter, point_class)) {
654 n_class_filtered++;
655 continue;
656 }
657 if (vector_mask_opt->answer) {
658 if (!VectorMask_point_in(&vector_mask, x, y)) {
659 n_outside_mask++;
660 continue;
661 }
662 }
663 if (count_decimation_is_out(&count_decimation_control))
664 continue;
665
666 Vect_append_point(Points, x, y, z);
667 if (id_layer)
668 Vect_cat_set(Cats, id_layer, cat);
669 if (return_layer) {
670 int return_c = return_to_cat(return_n, n_returns);
671 Vect_cat_set(Cats, return_layer, return_c);
672 }
673 if (class_layer) {
674 /* 0 is not a valid category and
675 * classes 0 and 1 as practically the same */
676 if (point_class == 0)
677 Vect_cat_set(Cats, class_layer, 1);
678 else
679 Vect_cat_set(Cats, class_layer, point_class);
680 }
681 if (have_color && rgb_layer) {
682 /* TODO: if attr table, acquired again, performance difference? */
683 LASColorH LAS_color = LASPoint_GetColor(LAS_point);
684 if (rgb_layer) {
685 int red = LASColor_GetRed(LAS_color);
686 int green = LASColor_GetGreen(LAS_color);
687 int blue = LASColor_GetBlue(LAS_color);
688 int rgb = red;
689 rgb = (rgb << 8) + green;
690 rgb = (rgb << 8) + blue;
691 rgb++; /* cat 0 is not valid, add one */
692 Vect_cat_set(Cats, rgb_layer, rgb);
693 }
694 }
695 Vect_write_line(&Map, GV_POINT, Points, Cats);
696
697 /* Attributes */
698 if (!notab_flag->answer) {
699 las_point_to_attributes(Fi, driver, cat, LAS_point, x, y, z,
700 have_time, have_color);
701 }
702
703 if (count_decimation_is_end(&count_decimation_control))
704 break;
705 if (id_layer && cat == GV_CAT_MAX) {
706 cat_max_reached = TRUE;
707 break;
708 }
709 cat++;
710 points_imported++;
711 }
712 G_percent(n_features, n_features, 1); /* finish it */
713
714 if (!notab_flag->answer) {
715 db_commit_transaction(driver);
716 db_close_database_shutdown_driver(driver);
717 }
718
719 if (vector_mask_opt->answer) {
720 VectorMask_destroy(&vector_mask);
721 }
722
723 LASSRS_Destroy(LAS_srs);
724 LASHeader_Destroy(LAS_header);
725 LASReader_Destroy(LAS_reader);
726
727 /* close map */
728 if (!notopo_flag->answer)
729 Vect_build(&Map);
730 Vect_close(&Map);
731
732 #ifdef HAVE_LONG_LONG_INT
733 unsigned long long not_valid_filtered = 0;
734 #else
735 unsigned long not_valid_filtered = 0;
736 #endif
737 if (only_valid)
738 not_valid_filtered = not_valid;
739
740 /* can be easily determined only when iterated over all points */
741 if (!count_decimation_control.limit_n && !cat_max_reached
742 && points_imported != n_features
743 - not_valid_filtered - n_outside - n_filtered - n_class_filtered
744 - n_outside_mask - count_decimation_control.offset_n_counter
745 - count_decimation_control.n_count_filtered - zrange_filtered)
746 G_warning(_("The underlying libLAS library is at its limits."
747 " Previously reported counts might have been distorted."
748 " However, the import itself should be unaffected."));
749
750 #ifdef HAVE_LONG_LONG_INT
751 if (count_decimation_control.limit_n) {
752 G_message(_("%llu points imported (limit was %llu)"),
753 count_decimation_control.limit_n_counter,
754 count_decimation_control.limit_n);
755 }
756 else {
757 G_message(_("%llu points imported"), points_imported);
758 }
759 if (not_valid && only_valid)
760 G_message(_("%llu input points were not valid and filtered out"), not_valid);
761 if (n_outside)
762 G_message(_("%llu input points were outside of the selected area"), n_outside);
763 if (n_outside_mask)
764 G_message(_("%llu input points were outside of the area specified by mask"), n_outside_mask);
765 if (n_filtered)
766 G_message(_("%llu input points were filtered out by return number"), n_filtered);
767 if (n_class_filtered)
768 G_message(_("%llu input points were filtered out by class number"), n_class_filtered);
769 if (zrange_filtered)
770 G_message(_("%llu input points were filtered outsite the range for z coordinate"), zrange_filtered);
771 if (count_decimation_control.offset_n_counter)
772 G_message(_("%llu input points were skipped at the begging using offset"),
773 count_decimation_control.offset_n_counter);
774 if (count_decimation_control.n_count_filtered)
775 G_message(_("%llu input points were skipped by count-based decimation"),
776 count_decimation_control.n_count_filtered);
777 #else
778 if (count_decimation_control.limit_n)
779 G_message(_("%lu points imported (limit was %d)"),
780 count_decimation_control.limit_n_counter,
781 count_decimation_control.limit_n);
782 else
783 G_message(_("%lu points imported"), points_imported);
784 if (not_valid && only_valid)
785 G_message(_("%lu input points were not valid and filtered out"), not_valid);
786 if (n_outside)
787 G_message(_("%lu input points were outside of the selected area"), n_outside);
788 if (n_outside_mask)
789 G_message(_("%lu input points were outside of the area specified by mask"), n_outside_mask);
790 if (n_filtered)
791 G_message(_("%lu input points were filtered out by return number"), n_filtered);
792 if (n_class_filtered)
793 G_message(_("%lu input points were filtered out by class number"), n_class_filtered);
794 if (zrange_filtered)
795 G_message(_("%lu input points were filtered outsite the range for z coordinate"), zrange_filtered);
796 if (count_decimation_control.offset_n_counter)
797 G_message(_("%lu input points were skipped at the begging using offset"),
798 count_decimation_control.offset_n_counter);
799 if (count_decimation_control.n_count_filtered)
800 G_message(_("%lu input points were skipped by count-based decimation"),
801 count_decimation_control.n_count_filtered);
802 G_message(_("Accuracy of the printed point counts might be limited by your computer architecture."));
803 #endif
804 if (count_decimation_control.limit_n)
805 G_message(_("The rest of points was ignored"));
806
807 #ifdef HAVE_LONG_LONG_INT
808 if (not_valid && !only_valid)
809 G_message(_("%llu input points were not valid, use -%c flag to filter"
810 " them out"), not_valid, only_valid_flag->key);
811 #else
812 if (not_valid && !only_valid)
813 G_message(_("%lu input points were not valid, use -%c flag to filter"
814 " them out"), not_valid, only_valid_flag->key);
815 #endif
816
817 if (cat_max_reached)
818 G_warning(_("Maximum number of categories reached (%d). Import ended prematurely."
819 " Try to import without using category as an ID."), GV_CAT_MAX);
820
821 /* -------------------------------------------------------------------- */
822 /* Extend current window based on dataset. */
823 /* -------------------------------------------------------------------- */
824 if (extend_flag->answer) {
825 G_get_set_window(&loc_wind);
826
827 loc_wind.north = MAX(loc_wind.north, cellhd.north);
828 loc_wind.south = MIN(loc_wind.south, cellhd.south);
829 loc_wind.west = MIN(loc_wind.west, cellhd.west);
830 loc_wind.east = MAX(loc_wind.east, cellhd.east);
831
832 loc_wind.rows = (int)ceil((loc_wind.north - loc_wind.south)
833 / loc_wind.ns_res);
834 loc_wind.south = loc_wind.north - loc_wind.rows * loc_wind.ns_res;
835
836 loc_wind.cols = (int)ceil((loc_wind.east - loc_wind.west)
837 / loc_wind.ew_res);
838 loc_wind.east = loc_wind.west + loc_wind.cols * loc_wind.ew_res;
839
840 G_put_window(&loc_wind);
841 }
842
843 exit(EXIT_SUCCESS);
844 }
845