[GRASS-SVN] r69379 - grass/trunk/raster3d/r3.in.lidar
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Sep 5 17:50:50 PDT 2016
Author: wenzeslaus
Date: 2016-09-05 17:50:49 -0700 (Mon, 05 Sep 2016)
New Revision: 69379
Added:
grass/trunk/raster3d/r3.in.lidar/info.c
grass/trunk/raster3d/r3.in.lidar/info.h
grass/trunk/raster3d/r3.in.lidar/projection.c
grass/trunk/raster3d/r3.in.lidar/projection.h
grass/trunk/raster3d/r3.in.lidar/string_list.c
grass/trunk/raster3d/r3.in.lidar/string_list.h
Modified:
grass/trunk/raster3d/r3.in.lidar/filters.c
grass/trunk/raster3d/r3.in.lidar/filters.h
grass/trunk/raster3d/r3.in.lidar/main.c
grass/trunk/raster3d/r3.in.lidar/rast_segment.c
grass/trunk/raster3d/r3.in.lidar/rast_segment.h
Log:
r3.in.lidar: multiple files input, projection check, scan modes, irange, scales (based on r.in.lidar) [news]
Modified: grass/trunk/raster3d/r3.in.lidar/filters.c
===================================================================
--- grass/trunk/raster3d/r3.in.lidar/filters.c 2016-09-05 16:16:36 UTC (rev 69378)
+++ grass/trunk/raster3d/r3.in.lidar/filters.c 2016-09-06 00:50:49 UTC (rev 69379)
@@ -1,11 +1,12 @@
/*
- * v.in.lidar filtering functions
+ * lidar-related filtering functions
*
- * Copyright 2011-2015 by Markus Metz, and The GRASS Development Team
* Authors:
* Markus Metz (r.in.lidar)
- * Vaclav Petras (move code to a separate files)
+ * Vaclav Petras (refactoring and various additions)
*
+ * Copyright 2011-2016 by Markus Metz, and The GRASS Development Team
+ *
* This program is free software licensed under the GPL (>=v2).
* Read the COPYING file that comes with GRASS for details.
*
@@ -19,7 +20,25 @@
#include <grass/gis.h>
#include <grass/glocale.h>
+int range_filter_from_option(struct Option *option, double *min, double *max)
+{
+ if (option->answer != NULL) {
+ if (option->answers[0] == NULL || option->answers[1] == NULL)
+ G_fatal_error(_("Invalid range <%s> for option %s"), option->answer, option->key);
+ sscanf(option->answers[0], "%lf", min);
+ sscanf(option->answers[1], "%lf", max);
+ /* for convenience, switch order to make valid input */
+ if (*min > *max) {
+ double tmp = *max;
+ *max = *min;
+ *min = tmp;
+ }
+ return TRUE;
+ }
+ return FALSE;
+}
+
int return_filter_create_from_string(struct ReturnFilter *return_filter,
const char *name)
{
Modified: grass/trunk/raster3d/r3.in.lidar/filters.h
===================================================================
--- grass/trunk/raster3d/r3.in.lidar/filters.h 2016-09-05 16:16:36 UTC (rev 69378)
+++ grass/trunk/raster3d/r3.in.lidar/filters.h 2016-09-06 00:50:49 UTC (rev 69379)
@@ -1,11 +1,12 @@
/*
- * v.in.lidar filtering functions
+ * lidar-related filtering functions
*
- * Copyright 2011-2015 by Markus Metz, and The GRASS Development Team
* Authors:
* Markus Metz (r.in.lidar)
- * Vaclav Petras (move code to a separate files)
+ * Vaclav Petras (refactoring and various additions)
*
+ * Copyright 2011-2016 by Markus Metz, and The GRASS Development Team
+ *
* This program is free software licensed under the GPL (>=v2).
* Read the COPYING file that comes with GRASS for details.
*
@@ -33,14 +34,8 @@
struct Option;
-int spatial_filter_from_option(struct Option *option, double *xmin,
- double *ymin, double *xmax, double *ymax);
-int spatial_filter_from_current_region(double *xmin, double *ymin,
- double *xmax, double *ymax);
+int range_filter_from_option(struct Option *option, double *min, double *max);
-int zrange_filter_from_option(struct Option *option,
- double *zmin, double *zmax);
-
int return_filter_create_from_string(struct ReturnFilter *return_filter,
const char *name);
int return_filter_is_out(struct ReturnFilter *return_filter, int return_n,
Copied: grass/trunk/raster3d/r3.in.lidar/info.c (from rev 69373, grass/trunk/raster/r.in.lidar/info.c)
===================================================================
--- grass/trunk/raster3d/r3.in.lidar/info.c (rev 0)
+++ grass/trunk/raster3d/r3.in.lidar/info.c 2016-09-06 00:50:49 UTC (rev 69379)
@@ -0,0 +1,180 @@
+/*
+ * lidar-related printing and scanning functions
+ *
+ * Authors:
+ * Markus Metz (r.in.lidar)
+ * Vaclav Petras (refactoring and various additions)
+ *
+ * Copyright 2011-2016 by Markus Metz, and The GRASS Development Team
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+#include <string.h>
+
+#include <grass/glocale.h>
+
+#include <liblas/capi/liblas.h>
+
+#include "local_proto.h"
+
+
+void print_lasinfo(LASHeaderH LAS_header, LASSRSH LAS_srs)
+{
+ char *las_srs_proj4 = LASSRS_GetProj4(LAS_srs);
+ int las_point_format = LASHeader_GetDataFormatId(LAS_header);
+
+ fprintf(stdout, "\nUsing LAS Library Version '%s'\n\n",
+ LAS_GetFullVersion());
+ fprintf(stdout, "LAS File Version: %d.%d\n",
+ LASHeader_GetVersionMajor(LAS_header),
+ LASHeader_GetVersionMinor(LAS_header));
+ fprintf(stdout, "System ID: '%s'\n",
+ LASHeader_GetSystemId(LAS_header));
+ fprintf(stdout, "Generating Software: '%s'\n",
+ LASHeader_GetSoftwareId(LAS_header));
+ fprintf(stdout, "File Creation Day/Year: %d/%d\n",
+ LASHeader_GetCreationDOY(LAS_header),
+ LASHeader_GetCreationYear(LAS_header));
+ fprintf(stdout, "Point Data Format: %d\n",
+ las_point_format);
+ fprintf(stdout, "Number of Point Records: %d\n",
+ LASHeader_GetPointRecordsCount(LAS_header));
+ fprintf(stdout, "Number of Points by Return: %d %d %d %d %d\n",
+ LASHeader_GetPointRecordsByReturnCount(LAS_header, 0),
+ LASHeader_GetPointRecordsByReturnCount(LAS_header, 1),
+ LASHeader_GetPointRecordsByReturnCount(LAS_header, 2),
+ LASHeader_GetPointRecordsByReturnCount(LAS_header, 3),
+ LASHeader_GetPointRecordsByReturnCount(LAS_header, 4));
+ fprintf(stdout, "Scale Factor X Y Z: %g %g %g\n",
+ LASHeader_GetScaleX(LAS_header),
+ LASHeader_GetScaleY(LAS_header), LASHeader_GetScaleZ(LAS_header));
+ fprintf(stdout, "Offset X Y Z: %g %g %g\n",
+ LASHeader_GetOffsetX(LAS_header),
+ LASHeader_GetOffsetY(LAS_header),
+ LASHeader_GetOffsetZ(LAS_header));
+ fprintf(stdout, "Min X Y Z: %g %g %g\n",
+ LASHeader_GetMinX(LAS_header),
+ LASHeader_GetMinY(LAS_header), LASHeader_GetMinZ(LAS_header));
+ fprintf(stdout, "Max X Y Z: %g %g %g\n",
+ LASHeader_GetMaxX(LAS_header),
+ LASHeader_GetMaxY(LAS_header), LASHeader_GetMaxZ(LAS_header));
+ if (las_srs_proj4 && strlen(las_srs_proj4) > 0) {
+ fprintf(stdout, "Spatial Reference:\n");
+ fprintf(stdout, "%s\n", las_srs_proj4);
+ }
+ else {
+ fprintf(stdout, "Spatial Reference: None\n");
+ }
+
+ fprintf(stdout, "\nData Fields:\n");
+ fprintf(stdout,
+ " 'X'\n 'Y'\n 'Z'\n 'Intensity'\n 'Return Number'\n");
+ fprintf(stdout, " 'Number of Returns'\n 'Scan Direction'\n");
+ fprintf(stdout,
+ " 'Flighline Edge'\n 'Classification'\n 'Scan Angle Rank'\n");
+ fprintf(stdout, " 'User Data'\n 'Point Source ID'\n");
+ if (las_point_format == 1 || las_point_format == 3 ||
+ las_point_format == 4 || las_point_format == 5) {
+ fprintf(stdout, " 'GPS Time'\n");
+ }
+ if (las_point_format == 2 || las_point_format == 3 ||
+ las_point_format == 5) {
+ fprintf(stdout, " 'Red'\n 'Green'\n 'Blue'\n");
+ }
+ fprintf(stdout, "\n");
+ fflush(stdout);
+
+ return;
+}
+
+
+int scan_bounds(LASReaderH LAS_reader, int shell_style, int extents, int update,
+ double zscale, struct Cell_head *region)
+{
+ unsigned long line;
+ int first;
+ double min_x, max_x, min_y, max_y, min_z, max_z;
+ double x, y, z;
+ LASPointH LAS_point;
+
+ line = 0;
+ first = TRUE;
+
+ /* init to nan in case no points are found */
+ min_x = max_x = min_y = max_y = min_z = max_z = 0.0 / 0.0;
+
+ G_verbose_message(_("Scanning data ..."));
+
+ LASReader_Seek(LAS_reader, 0);
+
+ while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
+ line++;
+
+ /* we don't do any filtering here */
+
+ x = LASPoint_GetX(LAS_point);
+ y = LASPoint_GetY(LAS_point);
+ z = LASPoint_GetZ(LAS_point);
+
+ if (first) {
+ min_x = x;
+ max_x = x;
+ min_y = y;
+ max_y = y;
+ min_z = z;
+ max_z = z;
+ first = FALSE;
+ }
+ else {
+ if (x < min_x)
+ min_x = x;
+ if (x > max_x)
+ max_x = x;
+ if (y < min_y)
+ min_y = y;
+ if (y > max_y)
+ max_y = y;
+ if (z < min_z)
+ min_z = z;
+ if (z > max_z)
+ max_z = z;
+ }
+ }
+
+ if (!extents) {
+ if (!shell_style) {
+ fprintf(stderr, _("Range: min max\n"));
+ fprintf(stdout, "x: %11f %11f\n", min_x, max_x);
+ fprintf(stdout, "y: %11f %11f\n", min_y, max_y);
+ fprintf(stdout, "z: %11f %11f\n", min_z * zscale, max_z * zscale);
+ }
+ else
+ fprintf(stdout, "n=%f s=%f e=%f w=%f b=%f t=%f\n",
+ max_y, min_y, max_x, min_x, min_z * zscale,
+ max_z * zscale);
+ }
+ 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;
+ region->north = max_y;
+ region->south = min_y;
+ }
+
+ G_debug(1, "Processed %lu points.", line);
+ G_debug(1, "region template: g.region n=%f s=%f e=%f w=%f",
+ max_y, min_y, max_x, min_x);
+ return 0;
+}
Added: grass/trunk/raster3d/r3.in.lidar/info.h
===================================================================
--- grass/trunk/raster3d/r3.in.lidar/info.h (rev 0)
+++ grass/trunk/raster3d/r3.in.lidar/info.h 2016-09-06 00:50:49 UTC (rev 69379)
@@ -0,0 +1,25 @@
+/*
+ * lidar-related printing and scanning functions
+ *
+ * Authors:
+ * Markus Metz (r.in.lidar)
+ * Vaclav Petras (refactoring and various additions)
+ *
+ * Copyright 2011-2016 by Vaclav Petras, and The GRASS Development Team
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+#ifndef __INFO_H__
+#define __INFO_H__
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <liblas/capi/liblas.h>
+
+void print_lasinfo(LASHeaderH, LASSRSH);
+int scan_bounds(LASReaderH, int, int, int, double, struct Cell_head *);
+
+#endif /* __INFO_H__ */
Property changes on: grass/trunk/raster3d/r3.in.lidar/info.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Modified: grass/trunk/raster3d/r3.in.lidar/main.c
===================================================================
--- grass/trunk/raster3d/r3.in.lidar/main.c 2016-09-05 16:16:36 UTC (rev 69378)
+++ grass/trunk/raster3d/r3.in.lidar/main.c 2016-09-06 00:50:49 UTC (rev 69379)
@@ -22,6 +22,9 @@
#include <grass/glocale.h>
#include <liblas/capi/liblas.h>
+#include "info.h"
+#include "string_list.h"
+#include "projection.h"
#include "rast_segment.h"
#include "filters.h"
@@ -124,12 +127,18 @@
{
struct GModule *module;
struct Option *input_opt;
+ struct Option *file_list_opt;
struct Option *count_output_opt, *sum_output_opt, *mean_output_opt;
struct Option *prop_count_output_opt, *prop_sum_output_opt;
struct Option *filter_opt, *class_opt;
+ struct Option *zscale_opt;
+ struct Option *iscale_opt;
+ struct Option *irange_opt;
struct Option *base_raster_opt;
struct Flag *base_rast_res_flag;
struct Flag *only_valid_flag;
+ struct Flag *print_flag, *scan_flag, *shell_style;
+ struct Flag *over_flag;
G_gisinit(argv[0]);
@@ -142,15 +151,22 @@
G_add_keyword(_("aggregation"));
G_add_keyword(_("binning"));
module->description =
- _("Creates a 3D raster map from LAS LiDAR points using univariate statistics.");
+ _("Creates a 3D raster map from LAS LiDAR points using univariate statistics.");
input_opt = G_define_standard_option(G_OPT_F_BIN_INPUT);
- input_opt->required = YES;
+ input_opt->required = NO;
input_opt->label = _("LAS input file");
input_opt->description =
_("LiDAR input file in LAS format (*.las or *.laz)");
input_opt->guisection = _("Input");
+ 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");
+
count_output_opt = G_define_standard_option(G_OPT_R3_OUTPUT);
count_output_opt->key = "n";
count_output_opt->required = YES;
@@ -223,6 +239,30 @@
base_rast_res_flag->description =
_("Use base raster actual resolution instead of computational region");
+ zscale_opt = G_define_option();
+ zscale_opt->key = "zscale";
+ zscale_opt->type = TYPE_DOUBLE;
+ zscale_opt->required = NO;
+ zscale_opt->answer = "1.0";
+ zscale_opt->description = _("Scale to apply to Z data");
+ zscale_opt->guisection = _("Transform");
+
+ irange_opt = G_define_option();
+ irange_opt->key = "intensity_range";
+ irange_opt->type = TYPE_DOUBLE;
+ irange_opt->required = NO;
+ irange_opt->key_desc = "min,max";
+ irange_opt->description = _("Filter range for intensity values (min,max)");
+ irange_opt->guisection = _("Selection");
+
+ iscale_opt = G_define_option();
+ iscale_opt->key = "intensity_scale";
+ iscale_opt->type = TYPE_DOUBLE;
+ iscale_opt->required = NO;
+ iscale_opt->answer = "1.0";
+ iscale_opt->description = _("Scale to apply to intensity values");
+ iscale_opt->guisection = _("Transform");
+
only_valid_flag = G_define_flag();
only_valid_flag->key = 'v';
only_valid_flag->label = _("Use only valid points");
@@ -231,29 +271,148 @@
" filtered out");
only_valid_flag->guisection = _("Selection");
+ over_flag = G_define_flag();
+ over_flag->key = 'o';
+ over_flag->label =
+ _("Override projection check (use current location's projection)");
+ over_flag->description =
+ _("Assume that the dataset has same projection as the current location");
+
+ print_flag = G_define_flag();
+ print_flag->key = 'p';
+ print_flag->description = _("Print LAS file info and exit");
+ print_flag->suppress_required = YES;
+
+ scan_flag = G_define_flag();
+ scan_flag->key = 's';
+ scan_flag->description = _("Scan data file for extent then exit");
+ scan_flag->suppress_required = YES;
+
+ shell_style = G_define_flag();
+ shell_style->key = 'g';
+ shell_style->description = _("In scan mode, print using shell script style");
+ shell_style->suppress_required = YES;
+
+ G_option_required(input_opt, file_list_opt, NULL);
+ G_option_exclusive(input_opt, file_list_opt, NULL);
G_option_requires(base_rast_res_flag, base_raster_opt, NULL);
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
+ if (shell_style->answer && !scan_flag->answer) {
+ scan_flag->answer = TRUE;
+ }
+
+ /* check intensity range and extent relation */
+ if (scan_flag->answer) {
+ if (irange_opt->answer)
+ G_warning(_("%s will not be taken into account during scan"), irange_opt->key);
+ }
+
int only_valid = FALSE;
if (only_valid_flag->answer)
only_valid = TRUE;
+ 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);
+ }
+
+ double zscale = 1.0;
+ double iscale = 1.0;
+ if (zscale_opt->answer)
+ zscale = atof(zscale_opt->answer);
+ if (iscale_opt->answer)
+ iscale = atof(iscale_opt->answer);
+
LASReaderH LAS_reader;
LASHeaderH LAS_header;
+ LASSRSH LAS_srs;
+ char *infile;
- LAS_reader = LASReader_Create(input_opt->answer);
- if (LAS_reader == NULL)
- G_fatal_error(_("Unable to open file <%s>"), input_opt->answer);
- LAS_header = LASReader_GetHeader(LAS_reader);
- if (LAS_header == NULL) {
+ /* for the CRS info */
+ const char *projstr;
+ struct Cell_head current_region;
+ struct Cell_head file_region;
+ G_get_set_window(¤t_region);
+
+ /* extent for all data */
+ struct Cell_head data_region;
+
+ long unsigned header_count = 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"),
- input_opt->answer);
+ infile);
+ }
+
+ LAS_srs = LASHeader_GetSRS(LAS_header);
+
+ /* 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(file_region, current_region, 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) {
+ /* we assign to the first one (i==0) but update for the rest */
+ scan_bounds(LAS_reader, shell_style->answer, FALSE, i,
+ zscale, &data_region);
+ }
+ /* number of estimated point across all files */
+ /* TODO: this should be ull which won't work with percent report */
+ header_count += 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);
+
Rast3d_init_defaults();
Rast3d_set_error_fun(Rast3d_fatal_error_noargs);
+ double irange_min;
+ double irange_max;
+ int use_irange =
+ range_filter_from_option(irange_opt, &irange_min, &irange_max);
+
struct ReturnFilter return_filter_struct;
int use_return_filter =
return_filter_create_from_string(&return_filter_struct,
@@ -291,6 +450,7 @@
int type = FCELL_TYPE;
int max_tile_size = 32;
+ /* TODO: this should probably happen before we change 2D region just to be sure */
Rast3d_get_window(&binning.region);
Rast3d_get_window(&binning.flat_region);
binning.flat_region.depths = 1;
@@ -368,73 +528,88 @@
long unsigned in_nulls = 0; /* or outside */
long unsigned n_return_filtered = 0;
long unsigned n_class_filtered = 0;
+ long unsigned irange_filtered = 0;
long unsigned n_invalid = 0;
long unsigned counter = 0;
- long unsigned header_count = LASHeader_GetPointRecordsCount(LAS_header);
-
G_verbose_message(_("Reading points..."));
G_percent_reset();
- while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
- if (counter == 100000) { /* report only some for speed */
- if (inside < header_count) /* TODO: inside can greatly underestimate */
- G_percent(inside, header_count, 3);
- counter = 0;
- }
- counter++;
- /* We always count them and report because r.in.lidar behavior
- * changed in between 7.0 and 7.2 from undefined (but skipping
- * invalid points) to filtering them out only when requested. */
- if (!LASPoint_IsValid(LAS_point)) {
- n_invalid++;
- if (only_valid)
- continue;
- }
- if (use_return_filter) {
- 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_return_filtered++;
- continue;
+ /* 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);
+ while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
+ if (counter == 100000) { /* report only some for speed */
+ if (inside < header_count) /* TODO: inside can greatly underestimate */
+ G_percent(inside, header_count, 3);
+ counter = 0;
}
- }
- if (use_class_filter) {
- int point_class = (int) LASPoint_GetClassification(LAS_point);
- if (class_filter_is_out(&class_filter, point_class)) {
- n_class_filtered++;
+ counter++;
+ /* We always count them and report because r.in.lidar behavior
+ * changed in between 7.0 and 7.2 from undefined (but skipping
+ * invalid points) to filtering them out only when requested. */
+ if (!LASPoint_IsValid(LAS_point)) {
+ n_invalid++;
+ if (only_valid)
+ continue;
+ }
+ if (use_return_filter) {
+ 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_return_filtered++;
+ continue;
+ }
+ }
+ if (use_class_filter) {
+ int point_class = (int) LASPoint_GetClassification(LAS_point);
+ if (class_filter_is_out(&class_filter, point_class)) {
+ n_class_filtered++;
+ continue;
+ }
+ }
+ value = LASPoint_GetIntensity(LAS_point);
+ value *= iscale;
+ if (use_irange && (value < irange_min || value > irange_max)) {
+ irange_filtered++;
continue;
}
- }
- east = LASPoint_GetX(LAS_point);
- north = LASPoint_GetY(LAS_point);
- top = LASPoint_GetZ(LAS_point);
+ east = LASPoint_GetX(LAS_point);
+ north = LASPoint_GetY(LAS_point);
+ top = LASPoint_GetZ(LAS_point);
+ top *= zscale;
- if (use_segment) {
- if (rast_segment_get_value_xy(&base_segment, &input_region,
- base_raster_data_type, east, north,
- &base_z)) {
- top -= base_z;
+ if (use_segment) {
+ if (rast_segment_get_value_xy(&base_segment, &input_region,
+ base_raster_data_type, east, north,
+ &base_z)) {
+ top -= base_z;
+ }
+ else {
+ /* TODO: separate message for this */
+ in_nulls += 1;
+ continue;
+ }
}
- else {
- in_nulls += 1;
+
+ Rast3d_location2coord(&binning.region, north, east, top, &col, &row,
+ &depth);
+ if (col >= binning.region.cols || row >= binning.region.rows ||
+ depth >= binning.region.depths || col < 0 || row < 0 ||
+ depth < 0) {
+ outside += 1;
continue;
}
+ binning_add_point(binning, row, col, depth, value);
+ inside += 1;
}
-
- Rast3d_location2coord(&binning.region, north, east, top, &col, &row,
- &depth);
- if (col >= binning.region.cols || row >= binning.region.rows ||
- depth >= binning.region.depths || col < 0 || row < 0 ||
- depth < 0) {
- outside += 1;
- continue;
- }
- value = LASPoint_GetIntensity(LAS_point);
- binning_add_point(binning, row, col, depth, value);
- inside += 1;
+ /* close input LAS file */
+ LASReader_Destroy(LAS_reader);
}
+ /* end of loop for all input files */
G_percent(1, 1, 1); /* flush */
@@ -483,6 +658,8 @@
G_message(_("%lu input points were filtered out by return number"), n_return_filtered);
if (n_class_filtered)
G_message(_("%lu input points were filtered out by class number"), n_class_filtered);
+ if (irange_filtered)
+ G_message(_("%lu input points had intensity out of range"), irange_filtered);
if (n_invalid && !only_valid)
G_message(_("%lu input points were not valid, use -%c flag to filter"
" them out"), n_invalid, only_valid_flag->key);
Copied: grass/trunk/raster3d/r3.in.lidar/projection.c (from rev 69373, grass/trunk/raster/r.in.lidar/projection.c)
===================================================================
--- grass/trunk/raster3d/r3.in.lidar/projection.c (rev 0)
+++ grass/trunk/raster3d/r3.in.lidar/projection.c 2016-09-06 00:50:49 UTC (rev 69379)
@@ -0,0 +1,144 @@
+/*
+ * lidar-related projection functions
+ *
+ * Authors:
+ * Markus Metz (r.in.lidar)
+ * Vaclav Petras (refactoring and various additions)
+ *
+ * Copyright 2011-2016 by Markus Metz, and The GRASS Development Team
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+
+#include <string.h>
+
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/gprojects.h>
+
+#include "local_proto.h"
+
+
+void projection_mismatch_report(struct Cell_head cellhd,
+ struct Cell_head loc_wind,
+ struct Key_Value *loc_proj_info,
+ struct Key_Value *loc_proj_units,
+ struct Key_Value *proj_info,
+ struct Key_Value *proj_units, int err)
+{
+ int i_value;
+ char error_msg[8192];
+
+ strcpy(error_msg,
+ _("Projection of dataset does not"
+ " appear to match current location.\n\n"));
+
+ /* TODO: output this info sorted by key: */
+ if (loc_wind.proj != cellhd.proj || err != -2) {
+ if (loc_proj_info != NULL) {
+ strcat(error_msg, _("GRASS LOCATION PROJ_INFO is:\n"));
+ for (i_value = 0; i_value < loc_proj_info->nitems; i_value++)
+ sprintf(error_msg + strlen(error_msg), "%s: %s\n",
+ loc_proj_info->key[i_value],
+ loc_proj_info->value[i_value]);
+ strcat(error_msg, "\n");
+ }
+
+ if (proj_info != NULL) {
+ strcat(error_msg, _("Import dataset PROJ_INFO is:\n"));
+ for (i_value = 0; i_value < proj_info->nitems; i_value++)
+ sprintf(error_msg + strlen(error_msg), "%s: %s\n",
+ proj_info->key[i_value], proj_info->value[i_value]);
+ }
+ else {
+ strcat(error_msg, _("Import dataset PROJ_INFO is:\n"));
+ if (cellhd.proj == PROJECTION_XY)
+ sprintf(error_msg + strlen(error_msg),
+ "Dataset proj = %d (unreferenced/unknown)\n",
+ cellhd.proj);
+ else if (cellhd.proj == PROJECTION_LL)
+ sprintf(error_msg + strlen(error_msg),
+ "Dataset proj = %d (lat/long)\n", cellhd.proj);
+ else if (cellhd.proj == PROJECTION_UTM)
+ sprintf(error_msg + strlen(error_msg),
+ "Dataset proj = %d (UTM), zone = %d\n",
+ cellhd.proj, cellhd.zone);
+ else
+ sprintf(error_msg + strlen(error_msg),
+ "Dataset proj = %d (unknown), zone = %d\n",
+ cellhd.proj, cellhd.zone);
+ }
+ }
+ else {
+ if (loc_proj_units != NULL) {
+ strcat(error_msg, "GRASS LOCATION PROJ_UNITS is:\n");
+ for (i_value = 0; i_value < loc_proj_units->nitems; i_value++)
+ sprintf(error_msg + strlen(error_msg), "%s: %s\n",
+ loc_proj_units->key[i_value],
+ loc_proj_units->value[i_value]);
+ strcat(error_msg, "\n");
+ }
+
+ if (proj_units != NULL) {
+ strcat(error_msg, "Import dataset PROJ_UNITS is:\n");
+ for (i_value = 0; i_value < proj_units->nitems; i_value++)
+ sprintf(error_msg + strlen(error_msg), "%s: %s\n",
+ proj_units->key[i_value], proj_units->value[i_value]);
+ }
+ }
+ sprintf(error_msg + strlen(error_msg),
+ _("\nIn case of no significant differences in the projection definitions,"
+ " use the -o flag to ignore them and use"
+ " current location definition.\n"));
+ strcat(error_msg,
+ _("Consider generating a new location with 'location' parameter"
+ " from input data set.\n"));
+ G_fatal_error("%s", error_msg);
+}
+
+void projection_check_wkt(struct Cell_head cellhd,
+ struct Cell_head loc_wind,
+ 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;
+ int err = 0;
+
+ proj_info = NULL;
+ proj_units = NULL;
+
+ /* Projection only required for checking so convert non-interactively */
+ if (GPJ_wkt_to_grass(&cellhd, &proj_info, &proj_units, projstr, 0) < 0)
+ G_warning(_("Unable to convert input map projection information to "
+ "GRASS format for checking"));
+
+ /* Does the projection of the current location match the dataset? */
+
+ /* fetch LOCATION PROJ info */
+ if (loc_wind.proj != PROJECTION_XY) {
+ loc_proj_info = G_get_projinfo();
+ loc_proj_units = G_get_projunits();
+ }
+
+ if (override) {
+ cellhd.proj = loc_wind.proj;
+ cellhd.zone = loc_wind.zone;
+ if (verbose)
+ G_message(_("Over-riding projection check"));
+ }
+ else if (loc_wind.proj != cellhd.proj
+ || (err =
+ G_compare_projections(loc_proj_info, loc_proj_units,
+ proj_info, proj_units)) != TRUE) {
+ projection_mismatch_report(cellhd, loc_wind, loc_proj_info,
+ loc_proj_units,
+ proj_info, proj_units, err);
+ }
+ else if (verbose) {
+ G_message(_("Projection of input dataset and current location "
+ "appear to match"));
+ }
+}
Added: grass/trunk/raster3d/r3.in.lidar/projection.h
===================================================================
--- grass/trunk/raster3d/r3.in.lidar/projection.h (rev 0)
+++ grass/trunk/raster3d/r3.in.lidar/projection.h 2016-09-06 00:50:49 UTC (rev 69379)
@@ -0,0 +1,34 @@
+/*
+ * lidar-related projection functions
+ *
+ * Authors:
+ * Markus Metz (r.in.lidar)
+ * Vaclav Petras (refactoring and various additions)
+ *
+ * Copyright 2011-2016 by Vaclav Petras, and The GRASS Development Team
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+#ifndef __PROJECTION_H__
+#define __PROJECTION_H__
+
+#include <grass/gis.h>
+
+
+void projection_mismatch_report(struct Cell_head cellhd,
+ struct Cell_head loc_wind,
+ struct Key_Value *loc_proj_info,
+ struct Key_Value *loc_proj_units,
+ struct Key_Value *proj_info,
+ struct Key_Value *proj_units,
+ int err);
+void projection_check_wkt(struct Cell_head cellhd,
+ struct Cell_head loc_wind,
+ const char *projstr,
+ int override,
+ int verbose);
+
+#endif /* __PROJECTION_H__ */
Property changes on: grass/trunk/raster3d/r3.in.lidar/projection.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Modified: grass/trunk/raster3d/r3.in.lidar/rast_segment.c
===================================================================
--- grass/trunk/raster3d/r3.in.lidar/rast_segment.c 2016-09-05 16:16:36 UTC (rev 69378)
+++ grass/trunk/raster3d/r3.in.lidar/rast_segment.c 2016-09-06 00:50:49 UTC (rev 69379)
@@ -1,4 +1,17 @@
+/*
+ * segment library convenient wrapper functions
+ *
+ * Authors:
+ * Vaclav Petras
+ *
+ * Copyright 2015 by Vaclav Petras, and the GRASS Development Team
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
#include <grass/gis.h>
#include <grass/glocale.h>
#include <grass/raster.h>
Modified: grass/trunk/raster3d/r3.in.lidar/rast_segment.h
===================================================================
--- grass/trunk/raster3d/r3.in.lidar/rast_segment.h 2016-09-05 16:16:36 UTC (rev 69378)
+++ grass/trunk/raster3d/r3.in.lidar/rast_segment.h 2016-09-06 00:50:49 UTC (rev 69379)
@@ -1,3 +1,17 @@
+/*
+ * segment library convenient wrapper functions
+ *
+ * Authors:
+ * Vaclav Petras
+ *
+ * Copyright 2015 by Vaclav Petras, and the GRASS Development Team
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+
#ifndef GRASS_RAST_SEGMENT_H
#define GRASS_RAST_SEGMENT_H
Copied: grass/trunk/raster3d/r3.in.lidar/string_list.c (from rev 69373, grass/trunk/raster/r.in.lidar/string_list.c)
===================================================================
--- grass/trunk/raster3d/r3.in.lidar/string_list.c (rev 0)
+++ grass/trunk/raster3d/r3.in.lidar/string_list.c 2016-09-06 00:50:49 UTC (rev 69379)
@@ -0,0 +1,78 @@
+/*
+ * Functionality to handle list of strings
+ *
+ * Authors:
+ * Vaclav Petras
+ *
+ * Copyright 2015 by Vaclav Petras, and the GRASS Development Team
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+
+#include <stdio.h>
+#include <string.h>
+
+#include <grass/gis.h>
+#include <grass/glocale.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 */
+ if (!file)
+ G_fatal_error(_("Cannot open file %s for reading"), filename);
+ char *line = G_malloc(GPATH_MAX * sizeof(char));
+
+ while (G_getl2(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;
+ string_list_add_item(string_list, strdup(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;
+}
Added: grass/trunk/raster3d/r3.in.lidar/string_list.h
===================================================================
--- grass/trunk/raster3d/r3.in.lidar/string_list.h (rev 0)
+++ grass/trunk/raster3d/r3.in.lidar/string_list.h 2016-09-06 00:50:49 UTC (rev 69379)
@@ -0,0 +1,31 @@
+/*
+ * Functionality to handle list of strings
+ *
+ * Authors:
+ * Vaclav Petras
+ *
+ * Copyright 2015 by Vaclav Petras, and the GRASS Development Team
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+
+#ifndef __STRING_LIST_H__
+#define __STRING_LIST_H__
+
+/* 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 /* __STRING_LIST_H__ */
Property changes on: grass/trunk/raster3d/r3.in.lidar/string_list.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list