[GRASS-SVN] r66601 - grass/trunk/raster/r.in.lidar
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Oct 25 21:18:55 PDT 2015
Author: wenzeslaus
Date: 2015-10-25 21:18:55 -0700 (Sun, 25 Oct 2015)
New Revision: 66601
Added:
grass/trunk/raster/r.in.lidar/string_list.c
Modified:
grass/trunk/raster/r.in.lidar/info.c
grass/trunk/raster/r.in.lidar/local_proto.h
grass/trunk/raster/r.in.lidar/main.c
grass/trunk/raster/r.in.lidar/projection.c
Log:
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 r66593, r66594, r66595 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.
Modified: grass/trunk/raster/r.in.lidar/info.c
===================================================================
--- grass/trunk/raster/r.in.lidar/info.c 2015-10-25 17:51:25 UTC (rev 66600)
+++ grass/trunk/raster/r.in.lidar/info.c 2015-10-26 04:18:55 UTC (rev 66601)
@@ -90,7 +90,7 @@
}
-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 @@
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;
Modified: grass/trunk/raster/r.in.lidar/local_proto.h
===================================================================
--- grass/trunk/raster/r.in.lidar/local_proto.h 2015-10-25 17:51:25 UTC (rev 66600)
+++ grass/trunk/raster/r.in.lidar/local_proto.h 2015-10-26 04:18:55 UTC (rev 66601)
@@ -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 @@
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__ */
Modified: grass/trunk/raster/r.in.lidar/main.c
===================================================================
--- grass/trunk/raster/r.in.lidar/main.c 2015-10-25 17:51:25 UTC (rev 66600)
+++ grass/trunk/raster/r.in.lidar/main.c 2015-10-26 04:18:55 UTC (rev 66601)
@@ -78,6 +78,7 @@
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,12 +103,22 @@
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";
method_opt->type = TYPE_STRING;
@@ -183,6 +194,7 @@
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 @@
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 @@
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);
+ /* 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"));
}
- /* 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);
- }
- LAS_srs = LASHeader_GetSRS(LAS_header);
+ Rast_get_window(®ion);
+ /* 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() */
- /* Print LAS header */
- if (print_flag->answer) {
- /* print... */
- print_lasinfo(LAS_header, LAS_srs);
+ 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, ®ion);
+ }
+ /* 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 @@
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 @@
if (point_binning.method == METHOD_N)
rtype = CELL_TYPE;
-
- Rast_get_window(®ion);
-
-
- 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, ®ion);
-
- 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 @@
* 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 @@
}
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 @@
/* 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 @@
/* 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,82 +509,93 @@
counter = 0;
G_percent_reset();
- while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
- line++;
- counter++;
+ /* 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);
- if (counter == 100000) { /* speed */
- if (line < estimated_lines)
- G_percent(line, estimated_lines, 3);
- counter = 0;
- }
+ while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
+ line++;
+ counter++;
- if (!LASPoint_IsValid(LAS_point)) {
- continue;
- }
+ if (counter == 100000) { /* speed */
+ if (line < estimated_lines)
+ G_percent(line, estimated_lines, 3);
+ counter = 0;
+ }
- 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);
+ if (!LASPoint_IsValid(LAS_point)) {
+ continue;
+ }
- 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;
+ 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);
- 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
+ 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;
- }
- 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
+ }
+ 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;
+ }
- count++;
- /* G_debug(5, "x: %f, y: %f, z: %f", x, y, z); */
+ z = z * zscale;
- update_value(&point_binning, &bin_index_nodes, cols, arr_row, arr_col, rtype, z);
- } /* while !EOF */
+ /* 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
+ 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 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);
count_total += count;
@@ -567,6 +611,7 @@
/* 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 @@
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);
Modified: grass/trunk/raster/r.in.lidar/projection.c
===================================================================
--- grass/trunk/raster/r.in.lidar/projection.c 2015-10-25 17:51:25 UTC (rev 66600)
+++ grass/trunk/raster/r.in.lidar/projection.c 2015-10-26 04:18:55 UTC (rev 66601)
@@ -99,7 +99,7 @@
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 @@
"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 @@
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 @@
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"));
}
Added: grass/trunk/raster/r.in.lidar/string_list.c
===================================================================
--- grass/trunk/raster/r.in.lidar/string_list.c (rev 0)
+++ grass/trunk/raster/r.in.lidar/string_list.c 2015-10-26 04:18:55 UTC (rev 66601)
@@ -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;
+}
Property changes on: grass/trunk/raster/r.in.lidar/string_list.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list