r.in.lidar: read multiple LAS files specified in a text file [news]

Printing and projection check is simply repeated for multiple files.

guisection Required and part of Optional moved to Input and Output.
input option preserved as single file input and made binary.

Related code in main moved together. Files open (and closed) as needed.
Files iterated on the level of points, so still writting the raster
just once.

Custom resolution and data-based extent supported. Base raster
segment reading fixed for the case of using current region.
Data-based extent is still limited by the current region
when base raster is used unless base raster resolution flag
is used.

Follows refactoring in https://trac.osgeo.org/grass/changeset/66593, https://trac.osgeo.org/grass/changeset/66594, https://trac.osgeo.org/grass/changeset/66595 and r66596.

Individual file extents could be preserved to avoid reading the
files for some row groups. Alternativelly, segement-based
reading and writting could be used.


git-svn-id: https://svn.osgeo.org/grass/grass/trunk@66601 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.in.lidar/info.c

@@ -90,7 +90,7 @@ void print_lasinfo(LASHeaderH LAS_header, LASSRSH LAS_srs)
 }
 
 
-int scan_bounds(LASReaderH LAS_reader, int shell_style, int extents,
+int scan_bounds(LASReaderH LAS_reader, int shell_style, int extents, int update,
                 double zscale, struct Cell_head *region)
 {
     unsigned long line;
@@ -161,6 +161,16 @@ int scan_bounds(LASReaderH LAS_reader, int shell_style, int extents,
         G_debug(1, "region template: g.region n=%f s=%f e=%f w=%f",
                 max_y, min_y, max_x, min_x);
     }
+    else if (update) {
+        if (min_x < region->west)
+            region->west = min_x;
+        if (max_x > region->east)
+            region->east = max_x;
+        if (min_y < region->south)
+            region->south = min_y;
+        if (max_y > region->north)
+            region->north = max_y;
+    }
     else {
         region->east = max_x;
         region->west = min_x;

raster/r.in.lidar/local_proto.h

@@ -49,7 +49,7 @@
 
 /* info.c */
 void print_lasinfo(LASHeaderH, LASSRSH);
-int scan_bounds(LASReaderH, int, int, double, struct Cell_head *);
+int scan_bounds(LASReaderH, int, int, int, double, struct Cell_head *);
 
 /* support.c */
 int blank_array(void *, int, int, RASTER_MAP_TYPE, int);
@@ -71,9 +71,22 @@ void projection_check_wkt(struct Cell_head cellhd,
                           struct Cell_head loc_wind,
                           const char *projstr,
                           int override,
-                          int shellstyle);
+                          int verbose);
 /* raster reading */
 int row_array_get_value_row_col(void *array, int arr_row, int arr_col,
                                 int cols, RASTER_MAP_TYPE rtype, double *value);
 
+/* multiple files */
+
+struct StringList
+{
+    int num_items;
+    int max_items;
+    char **items;
+};
+
+void string_list_from_file(struct StringList *string_list, char *filename);
+void string_list_from_one_item(struct StringList *string_list, char *item);
+void string_list_free(struct StringList *string_list);
+
 #endif /* __LOCAL_PROTO_H__ */

raster/r.in.lidar/main.c

@@ -78,6 +78,7 @@ int main(int argc, char *argv[])
     struct Option *input_opt, *output_opt, *percent_opt, *type_opt, *filter_opt, *class_opt;
     struct Option *method_opt, *base_raster_opt, *zrange_opt, *zscale_opt;
     struct Option *trim_opt, *pth_opt, *res_opt;
+    struct Option *file_list_opt;
     struct Flag *print_flag, *scan_flag, *shell_style, *over_flag, *extents_flag, *intens_flag;
     struct Flag *base_rast_res_flag;
 
@@ -102,11 +103,21 @@ int main(int argc, char *argv[])
     module->description =
 	_("Creates a raster map from LAS LiDAR points using univariate statistics.");
 
-    input_opt = G_define_standard_option(G_OPT_F_INPUT);
+    input_opt = G_define_standard_option(G_OPT_F_BIN_INPUT);
+    input_opt->required = NO;
     input_opt->label = _("LAS input file");
     input_opt->description = _("LiDAR input files in LAS format (*.las or *.laz)");
+    input_opt->guisection = _("Input");
 
     output_opt = G_define_standard_option(G_OPT_R_OUTPUT);
+    output_opt->guisection = _("Output");
+
+    file_list_opt = G_define_standard_option(G_OPT_F_INPUT);
+    file_list_opt->key = "file";
+    file_list_opt->label = _("File containing names of LAS input files");
+    file_list_opt->description = _("LiDAR input files in LAS format (*.las or *.laz)");
+    file_list_opt->required = NO;
+    file_list_opt->guisection = _("Input");
 
     method_opt = G_define_option();
     method_opt->key = "method";
@@ -183,6 +194,7 @@ int main(int argc, char *argv[])
     res_opt->required = NO;
     res_opt->description =
 	_("Output raster resolution");
+    res_opt->guisection = _("Output");
 
     filter_opt = G_define_option();
     filter_opt->key = "return_filter";
@@ -213,6 +225,7 @@ int main(int argc, char *argv[])
     extents_flag->key = 'e';
     extents_flag->description =
 	_("Extend region extents based on new dataset");
+    extents_flag->guisection = _("Output");
 
     over_flag = G_define_flag();
     over_flag->key = 'o';
@@ -244,44 +257,92 @@ int main(int argc, char *argv[])
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
+    struct StringList infiles;
+
+    if (file_list_opt->answer) {
+        if (access(file_list_opt->answer, F_OK) != 0)
+            G_fatal_error(_("File <%s> does not exist"), file_list_opt->answer);
+        string_list_from_file(&infiles, file_list_opt->answer);
+    }
+    else {
+        string_list_from_one_item(&infiles, input_opt->answer);
+    }
 
     /* parse input values */
-    infile = input_opt->answer;
     outmap = output_opt->answer;
 
     if (shell_style->answer && !scan_flag->answer) {
 	scan_flag->answer = 1; /* pointer not int, so set = shell_style->answer ? */
     }
 
-    /* Don't crash on cmd line if file not found */
-    if (access(infile, F_OK) != 0) {
-	G_fatal_error(_("Input file <%s> does not exist"), infile);
-    }
-    /* Open LAS file*/
-    LAS_reader = LASReader_Create(infile);
-    if (LAS_reader == NULL) {
-	G_fatal_error(_("Unable to open file <%s>"), infile);
-    }
-    
-    LAS_header = LASReader_GetHeader(LAS_reader);
-    if  (LAS_header == NULL) {
-	G_fatal_error(_("Input file <%s> is not a LAS LiDAR point cloud"),
-	                infile);
+    /* check zrange and extent relation */
+    if (scan_flag->answer || extents_flag->answer) {
+        if (zrange_opt->answer)
+            G_warning(_("zrange will not be taken into account during scan"));
     }
 
-    LAS_srs = LASHeader_GetSRS(LAS_header);
-
-    /* Print LAS header */
-    if (print_flag->answer) {
-	/* print... */
-	print_lasinfo(LAS_header, LAS_srs);
+    Rast_get_window(&region);
+    /* G_get_window seems to be unreliable if the location has been changed */
+    G_get_set_window(&loc_wind);        /* TODO: v.in.lidar uses G_get_default_window() */
+
+    estimated_lines = 0;
+    int i;
+    for (i = 0; i < infiles.num_items; i++) {
+        infile = infiles.items[i];
+        /* don't if file not found */
+        if (access(infile, F_OK) != 0)
+            G_fatal_error(_("Input file <%s> does not exist"), infile);
+        /* Open LAS file*/
+        LAS_reader = LASReader_Create(infile);
+        if (LAS_reader == NULL)
+            G_fatal_error(_("Unable to open file <%s>"), infile);
+        LAS_header = LASReader_GetHeader(LAS_reader);
+        if  (LAS_header == NULL) {
+            G_fatal_error(_("Input file <%s> is not a LAS LiDAR point cloud"),
+                            infile);
+        }
 
-	LASSRS_Destroy(LAS_srs);
-	LASHeader_Destroy(LAS_header);
-	LASReader_Destroy(LAS_reader);
+        LAS_srs = LASHeader_GetSRS(LAS_header);
 
-	exit(EXIT_SUCCESS);
+        /* print info or check projection if we are actually importing */
+        if (print_flag->answer) {
+            /* print filename when there is more than one file */
+            if (infiles.num_items > 1)
+                fprintf(stdout, "File: %s\n", infile);
+            /* Print LAS header */
+            print_lasinfo(LAS_header, LAS_srs);
+        }
+        else {
+            /* report that we are checking more files */
+            if (i == 1)
+                G_message(_("First file's projection checked,"
+                            " checking projection of the other files..."));
+            /* Fetch input map projection in GRASS form. */
+            projstr = LASSRS_GetWKT_CompoundOK(LAS_srs);
+            /* we are printing the non-warning messages only for first file */
+            projection_check_wkt(cellhd, loc_wind, projstr, over_flag->answer,
+                                 shell_style->answer || i);
+            /* if there is a problem in some other file, first OK message
+             * is printed but than a warning, this is not ideal but hopefully
+             * not so confusing when importing multiple files */
+        }
+        if (scan_flag->answer || extents_flag->answer) {
+            /* we assign to the first one (i==0) but update for the rest */
+            scan_bounds(LAS_reader, shell_style->answer, extents_flag->answer, i,
+                        zscale, &region);
+        }
+        /* number of estimated point across all files */
+        /* TODO: this should be ull which won't work with percent report */
+        estimated_lines += LASHeader_GetPointRecordsCount(LAS_header);
+        /* We are closing all again and we will be opening them later,
+         * so we don't have to worry about limit for open files. */
+        LASSRS_Destroy(LAS_srs);
+        LASHeader_Destroy(LAS_header);
+        LASReader_Destroy(LAS_reader);
     }
+    /* if we are not importing, end */
+    if (print_flag->answer || scan_flag->answer)
+        exit(EXIT_SUCCESS);
 
     return_filter = LAS_ALL;
     if (filter_opt->answer) {
@@ -299,13 +360,6 @@ int main(int argc, char *argv[])
     struct ClassFilter class_filter;
     class_filter_create_from_strings(&class_filter, class_opt->answers);
 
-    /* Fetch input map projection in GRASS form. */
-    projstr = LASSRS_GetWKT_CompoundOK(LAS_srs);
-
-    if (TRUE) {
-        projection_check_wkt(cellhd, loc_wind, projstr, over_flag->answer, shell_style->answer);
-    }
-
     percent = atoi(percent_opt->answer);
     zscale = atof(zscale_opt->answer);
 
@@ -338,25 +392,6 @@ int main(int argc, char *argv[])
     if (point_binning.method == METHOD_N)
 	rtype = CELL_TYPE;
 
-
-    Rast_get_window(&region);
-
-
-    if (scan_flag->answer || extents_flag->answer) {
-	if (zrange_opt->answer)
-	    G_warning(_("zrange will not be taken into account during scan"));
-
-	scan_bounds(LAS_reader, shell_style->answer, extents_flag->answer,
-	            zscale, &region);
-
-	if (!extents_flag->answer) {
-	    LASHeader_Destroy(LAS_header);
-	    LASReader_Destroy(LAS_reader);
-
-	    exit(EXIT_SUCCESS);
-	}
-    }
-    
     if (res_opt->answer) {
 	/* align to resolution */
 	res = atof(res_opt->answer);
@@ -396,9 +431,13 @@ int main(int argc, char *argv[])
      * region matches and using segment library when they don't */
     int use_segment = 0;
     int use_base_raster_res = 0;
+    /* TODO: see if the input region extent is smaller than the raster
+     * if yes, the we need to load the whole base raster if the -e
+     * flag was defined (alternatively clip the regions) */
     if (base_rast_res_flag->answer)
         use_base_raster_res = 1;
-    if (base_raster_opt->answer && (res_opt->answer || use_base_raster_res))
+    if (base_raster_opt->answer && (res_opt->answer || use_base_raster_res
+                                    || extents_flag->answer))
         use_segment = 1;
     if (base_raster_opt->answer && !use_segment) {
         /* TODO: do we need to test existence first? mapset? */
@@ -408,10 +447,12 @@ int main(int argc, char *argv[])
     }
     if (base_raster_opt->answer && use_segment) {
         if (use_base_raster_res) {
-            /* TODO: how to get cellhd already stored in the open map? */
+            /* read raster actual extent and resolution */
             Rast_get_cellhd(base_raster_opt->answer, "", &input_region);
             /* TODO: make it only as small as the output is or points are */
             Rast_set_input_window(&input_region);  /* we have split window */
+        } else {
+            Rast_get_input_window(&input_region);
         }
         rast_segment_open(&base_segment, base_raster_opt->answer, &base_raster_data_type);
     }
@@ -427,8 +468,6 @@ int main(int argc, char *argv[])
     /* open output map */
     out_fd = Rast_open_new(outmap, rtype);
 
-    estimated_lines = LASHeader_GetPointRecordsCount(LAS_header);
-
     /* allocate memory for a single row of output data */
     raster_row = Rast_allocate_output_buf(rtype);
 
@@ -441,16 +480,10 @@ int main(int argc, char *argv[])
 
     /* main binning loop(s) */
     for (pass = 1; pass <= npasses; pass++) {
-	LASError LAS_error;
 
 	if (npasses > 1)
 	    G_message(_("Pass #%d (of %d) ..."), pass, npasses);
 
-	LAS_error = LASReader_Seek(LAS_reader, 0);
-	
-	if (LAS_error != LE_None)
-	    G_fatal_error(_("Could not rewind input file"));
-
 	/* figure out segmentation */
 	pass_north = pass_south;  /* exact copy to avoid fp errors */
 	pass_south = region.north - pass * rows * region.ns_res;
@@ -476,81 +509,92 @@ int main(int argc, char *argv[])
 	counter = 0;
 	G_percent_reset();
 
-	while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
-	    line++;
-	    counter++;
-
-	    if (counter == 100000) {	/* speed */
-		if (line < estimated_lines)
-		    G_percent(line, estimated_lines, 3);
-		counter = 0;
-	    }
-
-	    if (!LASPoint_IsValid(LAS_point)) {
-		continue;
-	    }
-
-	    x = LASPoint_GetX(LAS_point);
-	    y = LASPoint_GetY(LAS_point);
-	    if (intens_flag->answer)
-		/* use z variable here to allow for scaling of intensity below */
-		z = LASPoint_GetIntensity(LAS_point);
-	    else
-		z = LASPoint_GetZ(LAS_point);
-
-        int return_n = LASPoint_GetReturnNumber(LAS_point);
-        int n_returns = LASPoint_GetNumberOfReturns(LAS_point);
-        if (return_filter_is_out(&return_filter_struct, return_n, n_returns)) {
-            n_filtered++;
-            continue;
-        }
-        point_class = (int) LASPoint_GetClassification(LAS_point);
-        if (class_filter_is_out(&class_filter, point_class))
-            continue;
-
-	    if (y <= pass_south || y > pass_north) {
-		continue;
-	    }
-	    if (x < region.west || x >= region.east) {
-		continue;
-	    }
-
-	    z = z * zscale;
-
-            /* find the bin in the current array box */
-            arr_row = (int)((pass_north - y) / region.ns_res);
-            arr_col = (int)((x - region.west) / region.ew_res);
-
-            if (base_array) {
-                double base_z;
-                if (row_array_get_value_row_col(base_array, arr_row, arr_col,
-                                                cols, base_raster_data_type,
-                                                &base_z))
-                    z -= base_z;
-                else
+        /* loop of input files */
+        for (i = 0; i < infiles.num_items; i++) {
+            infile = infiles.items[i];
+            /* we already know file is there, so just do basic checks */
+            LAS_reader = LASReader_Create(infile);
+            if (LAS_reader == NULL)
+                G_fatal_error(_("Unable to open file <%s>"), infile);
+
+            while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
+                line++;
+                counter++;
+
+                if (counter == 100000) {        /* speed */
+                    if (line < estimated_lines)
+                        G_percent(line, estimated_lines, 3);
+                    counter = 0;
+                }
+
+                if (!LASPoint_IsValid(LAS_point)) {
                     continue;
-            }
-            else if (use_segment) {
-                double base_z;
-                if (rast_segment_get_value_xy(&base_segment, &input_region,
-                                              base_raster_data_type, x, y,
-                                              &base_z))
-                    z -= base_z;
+                }
+
+                x = LASPoint_GetX(LAS_point);
+                y = LASPoint_GetY(LAS_point);
+                if (intens_flag->answer)
+                    /* use z variable here to allow for scaling of intensity below */
+                    z = LASPoint_GetIntensity(LAS_point);
                 else
+                    z = LASPoint_GetZ(LAS_point);
+
+                int return_n = LASPoint_GetReturnNumber(LAS_point);
+                int n_returns = LASPoint_GetNumberOfReturns(LAS_point);
+                if (return_filter_is_out(&return_filter_struct, return_n, n_returns)) {
+                    n_filtered++;
+                    continue;
+                }
+                point_class = (int) LASPoint_GetClassification(LAS_point);
+                if (class_filter_is_out(&class_filter, point_class))
                     continue;
-            }
 
-            if (zrange_opt->answer) {
-                if (z < zrange_min || z > zrange_max) {
+                if (y <= pass_south || y > pass_north) {
                     continue;
                 }
-            }
+                if (x < region.west || x >= region.east) {
+                    continue;
+                }
+
+                z = z * zscale;
+
+                /* find the bin in the current array box */
+                arr_row = (int)((pass_north - y) / region.ns_res);
+                arr_col = (int)((x - region.west) / region.ew_res);
 
-	    count++;
-	    /*          G_debug(5, "x: %f, y: %f, z: %f", x, y, z); */
+                if (base_array) {
+                    double base_z;
+                    if (row_array_get_value_row_col(base_array, arr_row, arr_col,
+                                                    cols, base_raster_data_type,
+                                                    &base_z))
+                        z -= base_z;
+                    else
+                        continue;
+                }
+                else if (use_segment) {
+                    double base_z;
+                    if (rast_segment_get_value_xy(&base_segment, &input_region,
+                                                  base_raster_data_type, x, y,
+                                                  &base_z))
+                        z -= base_z;
+                    else
+                        continue;
+                }
+
+                if (zrange_opt->answer) {
+                    if (z < zrange_min || z > zrange_max) {
+                        continue;
+                    }
+                }
+
+                count++;
+                /*          G_debug(5, "x: %f, y: %f, z: %f", x, y, z); */
 
-	    update_value(&point_binning, &bin_index_nodes, cols, arr_row, arr_col, rtype, z);
-	}			/* while !EOF */
+                update_value(&point_binning, &bin_index_nodes, cols, arr_row, arr_col, rtype, z);
+            }                        /* while !EOF of one input file */
+            /* close input LAS file */
+            LASReader_Destroy(LAS_reader);
+        }           /* end of loop for all input files files */
 
 	G_percent(1, 1, 1);	/* flush */
 	G_debug(2, "pass %d finished, %lu coordinates in box", pass, count);
@@ -567,6 +611,7 @@ int main(int argc, char *argv[])
 
 	/* free memory */
 	point_binning_free(&point_binning, &bin_index_nodes);
+    string_list_free(&infiles);
     }				/* passes loop */
     if (base_array)
         Rast_close(base_raster);
@@ -576,10 +621,6 @@ int main(int argc, char *argv[])
     G_percent(1, 1, 1);		/* flush */
     G_free(raster_row);
 
-    /* close input LAS file */
-    LASHeader_Destroy(LAS_header);
-    LASReader_Destroy(LAS_reader);
-
     /* close raster file & write history */
     Rast_close(out_fd);
 

raster/r.in.lidar/projection.c

@@ -99,7 +99,7 @@ void projection_mismatch_report(struct Cell_head cellhd,
 
 void projection_check_wkt(struct Cell_head cellhd,
                           struct Cell_head loc_wind,
-                          const char *projstr, int override, int shellstyle)
+                          const char *projstr, int override, int verbose)
 {
     struct Key_Value *loc_proj_info = NULL, *loc_proj_units = NULL;
     struct Key_Value *proj_info, *proj_units;
@@ -114,8 +114,7 @@ void projection_check_wkt(struct Cell_head cellhd,
                     "GRASS format for checking"));
 
     /* Does the projection of the current location match the dataset? */
-    /* G_get_window seems to be unreliable if the location has been changed */
-    G_get_set_window(&loc_wind);        /* TODO: v.in.lidar uses G_get_default_window() */
+
     /* fetch LOCATION PROJ info */
     if (loc_wind.proj != PROJECTION_XY) {
         loc_proj_info = G_get_projinfo();
@@ -125,7 +124,8 @@ void projection_check_wkt(struct Cell_head cellhd,
     if (override) {
         cellhd.proj = loc_wind.proj;
         cellhd.zone = loc_wind.zone;
-        G_message(_("Over-riding projection check"));
+        if (verbose)
+            G_message(_("Over-riding projection check"));
     }
     else if (loc_wind.proj != cellhd.proj
              || (err =
@@ -135,7 +135,7 @@ void projection_check_wkt(struct Cell_head cellhd,
                                    loc_proj_units,
                                    proj_info, proj_units, err);
     }
-    else if (!shellstyle) {
+    else if (verbose) {
         G_message(_("Projection of input dataset and current location "
                     "appear to match"));
     }

raster/r.in.lidar/string_list.c

@@ -0,0 +1,65 @@
+
+#include <stdio.h>
+#include <string.h>
+
+#include <grass/gis.h>
+
+#include "local_proto.h"
+
+#define SIZE_INCREMENT 10
+
+static int string_list_add_item(struct StringList *string_list, char *item)
+{
+    int n = string_list->num_items++;
+
+    if (string_list->num_items >= string_list->max_items) {
+        string_list->max_items += SIZE_INCREMENT;
+        string_list->items = G_realloc(string_list->items,
+                                       (size_t) string_list->max_items *
+                                       sizeof(char *));
+    }
+    /* n contains the index */
+    string_list->items[n] = item;
+    return n;
+}
+
+void string_list_from_file(struct StringList *string_list, char *filename)
+{
+    string_list->num_items = 0;
+    string_list->max_items = 0;
+    string_list->items = NULL;
+    FILE *file = fopen(filename, "r");  /* should check the result */
+    char *line = G_malloc(GPATH_MAX * sizeof(char));
+
+    while (G_getl(line, GPATH_MAX * sizeof(char), file)) {
+        G_debug(5, "line content from file %s: %s\n", filename, line);
+        string_list_add_item(string_list, line);
+        line = G_malloc(GPATH_MAX);
+    }
+    /* last allocation was not necessary */
+    G_free(line);
+    fclose(file);
+}
+
+void string_list_from_one_item(struct StringList *string_list, char *item)
+{
+    string_list->num_items = 0;
+    string_list->max_items = 0;
+    string_list->items = NULL;
+    char *new_item = G_malloc(strlen(item) * sizeof(char));
+
+    strcpy(new_item, item);
+    string_list_add_item(string_list, new_item);
+}
+
+void string_list_free(struct StringList *string_list)
+{
+    int i;
+
+    for (i = 0; i < string_list->num_items; i++)
+        G_free(string_list->items[i]);
+    G_free(string_list->items);
+    string_list->num_items = 0;
+    string_list->max_items = 0;
+    string_list->items = NULL;
+}