[GRASS-SVN] r67705 - in grass/trunk/raster3d: . r3.in.lidar
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Jan 29 13:17:46 PST 2016
Author: wenzeslaus
Date: 2016-01-29 13:17:46 -0800 (Fri, 29 Jan 2016)
New Revision: 67705
Added:
grass/trunk/raster3d/r3.in.lidar/
grass/trunk/raster3d/r3.in.lidar/Makefile
grass/trunk/raster3d/r3.in.lidar/main.c
grass/trunk/raster3d/r3.in.lidar/r3.in.lidar.html
grass/trunk/raster3d/r3.in.lidar/r3_in_lidar.png
Modified:
grass/trunk/raster3d/Makefile
grass/trunk/raster3d/raster3dintro.html
Log:
r3.in.lidar: basic implementation [news]
Experimental, major limitations described in manual.
Modified: grass/trunk/raster3d/Makefile
===================================================================
--- grass/trunk/raster3d/Makefile 2016-01-29 09:56:07 UTC (rev 67704)
+++ grass/trunk/raster3d/Makefile 2016-01-29 21:17:46 UTC (rev 67705)
@@ -7,6 +7,7 @@
r3.gwflow \
r3.in.ascii \
r3.in.bin \
+ r3.in.lidar \
r3.in.v5d \
r3.info \
r3.mask \
Added: grass/trunk/raster3d/r3.in.lidar/Makefile
===================================================================
--- grass/trunk/raster3d/r3.in.lidar/Makefile (rev 0)
+++ grass/trunk/raster3d/r3.in.lidar/Makefile 2016-01-29 21:17:46 UTC (rev 67705)
@@ -0,0 +1,15 @@
+MODULE_TOPDIR = ../..
+
+PGM = r3.in.lidar
+
+LIBES = $(RASTERLIB) $(RASTER3DLIB) $(GISLIB) $(MATHLIB) $(GPROJLIB) $(LASLIBS) $(SEGMENTLIB) $(VECTORLIB)
+DEPENDENCIES = $(RASTERDEP) $(RASTER3DDEP) $(GISDEP) $(SEGMENTDEP) $(VECTORDEP)
+
+EXTRA_INC = $(LASINC) $(VECT_INC) $(PROJINC)
+EXTRA_CFLAGS = $(VECT_CFLAGS) $(GDALCFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+ifneq ($(USE_LIBLAS),)
+default: cmd
+endif
Property changes on: grass/trunk/raster3d/r3.in.lidar/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass/trunk/raster3d/r3.in.lidar/main.c
===================================================================
--- grass/trunk/raster3d/r3.in.lidar/main.c (rev 0)
+++ grass/trunk/raster3d/r3.in.lidar/main.c 2016-01-29 21:17:46 UTC (rev 67705)
@@ -0,0 +1,264 @@
+/****************************************************************************
+*
+* MODULE: r3.in.Lidar
+*
+* AUTHOR(S): Vaclav Petras
+*
+* PURPOSE: Imports LAS LiDAR point clouds to a 3D raster map using
+* aggregate statistics.
+*
+* COPYRIGHT: (C) 2016 Vaclav Petras and the The GRASS Development Team
+*
+* This program is free software under the GNU General Public
+* License (>=v2). Read the file COPYING that comes with GRASS
+* for details.
+*
+*****************************************************************************/
+
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/raster3d.h>
+#include <grass/glocale.h>
+#include <liblas/capi/liblas.h>
+
+static void raster3d_set_value_float(RASTER3D_Map *raster, RASTER3D_Region *region, float value)
+{
+ int col, row, depth;
+
+ for (depth = 0; depth < region->depths; depth++)
+ for (row = 0; row < region->rows; row++)
+ for (col = 0; col < region->cols; col++)
+ Rast3d_put_float(raster, col, row, depth, value);
+}
+
+/* c = a / b */
+static void raster3d_divide(RASTER3D_Map *a, RASTER3D_Map *b, RASTER3D_Map *c, RASTER3D_Region *region)
+{
+ int col, row, depth;
+ double tmp;
+
+ for (depth = 0; depth < region->depths; depth++)
+ for (row = 0; row < region->rows; row++)
+ for (col = 0; col < region->cols; col++) {
+ tmp = Rast3d_get_double(b, col, row, 0);
+ /* TODO: compare to epsilon */
+ if (tmp > 0) {
+ tmp = Rast3d_get_double(a, col, row, depth) / tmp;
+ Rast3d_put_double(c, col, row, depth, tmp);
+ }
+ else {
+ /* TODO: check this implementation */
+ Rast3d_set_null_value(&tmp, 1, DCELL_TYPE);
+ Rast3d_put_double(c, col, row, depth, tmp);
+ }
+ }
+}
+
+/* c = a / b where b has depth 1 */
+static void raster3d_divide_by_flat(RASTER3D_Map *a, RASTER3D_Map *b, RASTER3D_Map *c, RASTER3D_Region *region)
+{
+ int col, row, depth;
+ double tmp;
+
+ for (depth = 0; depth < region->depths; depth++)
+ for (row = 0; row < region->rows; row++)
+ for (col = 0; col < region->cols; col++) {
+ tmp = Rast3d_get_double(b, col, row, 0);
+ /* since it is count, using cast to integer to check
+ againts zero, limits the value to max of CELL */
+ if (((CELL) tmp) > 0) {
+ tmp = Rast3d_get_double(a, col, row, depth) / tmp;
+ Rast3d_put_double(c, col, row, depth, tmp);
+ }
+ else {
+ /* TODO: check this implementation */
+ Rast3d_set_null_value(&tmp, 1, DCELL_TYPE);
+ Rast3d_put_double(c, col, row, depth, tmp);
+ }
+ }
+}
+
+int main(int argc, char *argv[])
+{
+ struct GModule *module;
+ struct Option *input_opt;
+ struct Option *count_output_opt, *sum_output_opt, *mean_output_opt;
+ struct Option *prop_count_output_opt, *prop_sum_output_opt;
+
+ G_gisinit(argv[0]);
+ module = G_define_module();
+ G_add_keyword(_("3D raster"));
+ G_add_keyword(_("import"));
+ G_add_keyword(_("LIDAR"));
+ module->description =
+ _("Creates a 3D raster map from LAS LiDAR points");
+
+ input_opt = G_define_standard_option(G_OPT_F_BIN_INPUT);
+ input_opt->required = YES;
+ input_opt->label = _("LAS input file");
+ input_opt->description = _("LiDAR input file in LAS format (*.las or *.laz)");
+ input_opt->guisection = _("Input");
+
+ count_output_opt = G_define_standard_option(G_OPT_R3_OUTPUT);
+ count_output_opt->key = "n";
+ count_output_opt->required = YES;
+ count_output_opt->label = _("Count of points per cell");
+ count_output_opt->guisection = _("Output");
+
+ sum_output_opt = G_define_standard_option(G_OPT_R3_OUTPUT);
+ sum_output_opt->key = "sum";
+ sum_output_opt->required = YES;
+ sum_output_opt->label = _("Sum of values of point intensities per cell");
+ sum_output_opt->guisection = _("Output");
+
+ mean_output_opt = G_define_standard_option(G_OPT_R3_OUTPUT);
+ mean_output_opt->key = "mean";
+ mean_output_opt->required = YES;
+ mean_output_opt->label = _("Mean of point intensities per cell");
+ mean_output_opt->guisection = _("Output");
+
+ /* TODO: proportional versus relative in naming */
+ prop_count_output_opt = G_define_standard_option(G_OPT_R3_OUTPUT);
+ prop_count_output_opt->key = "proportional_n";
+ prop_count_output_opt->required = YES;
+ prop_count_output_opt->label =
+ _("3D raster map of proportional point count");
+ prop_count_output_opt->description =
+ _("Point count per 3D cell divided by point count per vertical"
+ " column");
+ prop_count_output_opt->guisection = _("Proportional output");
+
+ prop_sum_output_opt = G_define_standard_option(G_OPT_R3_OUTPUT);
+ prop_sum_output_opt->key = "proportional_sum";
+ prop_sum_output_opt->required = YES;
+ prop_sum_output_opt->label =
+ _("3D raster map of proportional sum of values");
+ prop_sum_output_opt->description =
+ _("Sum of values per 3D cell divided by sum of values per"
+ " vertical column");
+ prop_sum_output_opt->guisection = _("Proportional output");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ LASReaderH LAS_reader;
+ LAS_reader = LASReader_Create(input_opt->answer);
+ if (LAS_reader == NULL)
+ G_fatal_error(_("Unable to open file <%s>"), input_opt->answer);
+
+ Rast3d_init_defaults();
+ Rast3d_set_error_fun(Rast3d_fatal_error_noargs);
+
+ RASTER3D_Region region, flat_region;
+ RASTER3D_Map *count_raster, *sum_raster, *mean_raster;
+ RASTER3D_Map *count_flat_raster, *sum_flat_raster;
+ RASTER3D_Map *prop_count_raster, *prop_sum_raster;
+
+ Rast3d_get_window(®ion);
+ Rast3d_get_window(&flat_region);
+ flat_region.depths = 1;
+ Rast3d_adjust_region(&flat_region);
+
+ count_raster = Rast3d_open_new_opt_tile_size(count_output_opt->answer,
+ RASTER3D_USE_CACHE_DEFAULT,
+ ®ion, FCELL_TYPE, 32);
+ if (!count_raster)
+ Rast3d_fatal_error(_("Unable to create 3D raster map <%s>"),
+ count_output_opt->answer);
+ sum_raster = Rast3d_open_new_opt_tile_size(sum_output_opt->answer,
+ RASTER3D_USE_CACHE_DEFAULT,
+ ®ion, FCELL_TYPE, 32);
+ if (!sum_raster)
+ Rast3d_fatal_error(_("Unable to create 3D raster map <%s>"),
+ sum_output_opt->answer);
+ mean_raster = Rast3d_open_new_opt_tile_size(mean_output_opt->answer,
+ RASTER3D_USE_CACHE_DEFAULT,
+ ®ion, FCELL_TYPE, 32);
+ if (!mean_raster)
+ Rast3d_fatal_error(_("Unable to create 3D raster map <%s>"),
+ mean_output_opt->answer);
+ count_flat_raster = Rast3d_open_new_opt_tile_size("r3_in_lidar_tmp_sum_flat",
+ RASTER3D_USE_CACHE_DEFAULT,
+ &flat_region, FCELL_TYPE, 32);
+ if (!count_flat_raster)
+ Rast3d_fatal_error(_("Unable to create 3D raster map <%s>"),
+ count_output_opt->answer);
+ sum_flat_raster = Rast3d_open_new_opt_tile_size("r3_in_lidar_tmp_count_flat",
+ RASTER3D_USE_CACHE_DEFAULT,
+ &flat_region, FCELL_TYPE, 32);
+ if (!sum_flat_raster)
+ Rast3d_fatal_error(_("Unable to create 3D raster map <%s>"),
+ count_output_opt->answer);
+ prop_count_raster = Rast3d_open_new_opt_tile_size(prop_count_output_opt->answer,
+ RASTER3D_USE_CACHE_DEFAULT,
+ ®ion, FCELL_TYPE, 32);
+ if (!prop_count_raster)
+ Rast3d_fatal_error(_("Unable to create 3D raster map <%s>"),
+ prop_count_output_opt->answer);
+ prop_sum_raster = Rast3d_open_new_opt_tile_size(prop_sum_output_opt->answer,
+ RASTER3D_USE_CACHE_DEFAULT,
+ ®ion, FCELL_TYPE, 32);
+ if (!prop_sum_raster)
+ Rast3d_fatal_error(_("Unable to create 3D raster map <%s>"),
+ prop_sum_output_opt->answer);
+
+ raster3d_set_value_float(count_raster, ®ion, 0);
+ raster3d_set_value_float(sum_raster, ®ion, 0);
+ raster3d_set_value_float(count_flat_raster, &flat_region, 0);
+ raster3d_set_value_float(sum_flat_raster, &flat_region, 0);
+
+ LASPointH LAS_point;
+ double east, north, top;
+ int row, col, depth;
+ double value;
+ double tmp;
+
+ /* TODO: use long long */
+ long unsigned inside = 0;
+ long unsigned outside = 0;
+
+ while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
+ if (!LASPoint_IsValid(LAS_point))
+ continue;
+
+ east = LASPoint_GetX(LAS_point);
+ north = LASPoint_GetY(LAS_point);
+ top = LASPoint_GetZ(LAS_point);
+
+ if (!Rast3d_is_valid_location(®ion, north, east, top)) {
+ outside += 1;
+ continue;
+ }
+ Rast3d_location2coord(®ion, north, east, top, &col, &row, &depth);
+ value = LASPoint_GetIntensity(LAS_point);
+
+ tmp = Rast3d_get_double(count_raster, col, row, depth);
+ Rast3d_put_double(count_raster, col, row, depth, tmp + 1);
+ tmp = Rast3d_get_double(count_flat_raster, col, row, 0);
+ Rast3d_put_double(count_flat_raster, col, row, 0, tmp + 1);
+ tmp = Rast3d_get_double(sum_raster, col, row, depth);
+ Rast3d_put_double(sum_raster, col, row, depth, tmp + value);
+ tmp = Rast3d_get_double(sum_flat_raster, col, row, 0);
+ Rast3d_put_double(sum_flat_raster, col, row, 0, tmp + value);
+
+ inside += 1;
+ }
+
+ raster3d_divide_by_flat(count_raster, count_flat_raster, prop_count_raster, ®ion);
+ raster3d_divide_by_flat(sum_raster, sum_flat_raster, prop_sum_raster, ®ion);
+
+ raster3d_divide(sum_raster, count_raster, mean_raster, ®ion);
+
+ G_message("Number of point inside: %lu", inside);
+ G_message("Number of point outside: %lu", outside);
+
+ Rast3d_close(prop_sum_raster);
+ Rast3d_close(prop_count_raster);
+ Rast3d_close(sum_flat_raster); /* TODO: delete */
+ Rast3d_close(count_flat_raster); /* TODO: delete */
+ Rast3d_close(mean_raster);
+ Rast3d_close(sum_raster);
+ Rast3d_close(count_raster);
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass/trunk/raster3d/r3.in.lidar/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass/trunk/raster3d/r3.in.lidar/r3.in.lidar.html
===================================================================
--- grass/trunk/raster3d/r3.in.lidar/r3.in.lidar.html (rev 0)
+++ grass/trunk/raster3d/r3.in.lidar/r3.in.lidar.html 2016-01-29 21:17:46 UTC (rev 67705)
@@ -0,0 +1,74 @@
+<h2>DESCRIPTION</h2>
+
+<p>
+<img src="r3_in_lidar.png">
+<p>
+<em>
+ Figure: Proportional count of points per 3D cell. When 50% of all
+ points in a vertical column fall into a given 3D cell, the value
+ is 0.5. Here, the green color was assigned to 0.5, red to 1 and
+ yellow to 0. The figure shows vertical slices and green color
+ indicates high vegetation while red color indicates bare ground.
+</em>
+<!--
+0% 255:255:100
+50% green
+100% red
+-->
+
+
+<h2>NOTES</h2>
+
+<ul>
+ <li>
+ This module is highly experimental. Don't rely on its
+ functionality or interface. Please report issues on the mailing
+ list or in the bug tracker.
+ <li>
+ No projection check or reprojection is performed, make sure you
+ are using data in the right projection for your GRASS Location.
+ <li>
+ Selection of points according to return or class is not yet
+ supported.
+ <li>All outputs are currently mandatory.
+ <li>Some temporary maps are created but not cleaned up.
+ <li>
+ Expects points to have intensity (random result for related
+ outputs when they don't).
+</ul>
+
+
+<h2>EXAMPLES</h2>
+
+Set the region according to a 2D raster and adding 3D minimum
+(bottom), maximum (top) and vertical (top-bottom) resolution.
+
+<div class="code"><pre>
+g.region rast=secref b=80 t=160 tbres=5 -p3
+</pre></div>
+
+Now, <em>r3.in.lidar</em> will create the 3D raster of the size
+given by the computation region:
+
+<div class="code"><pre>
+r3.in.lidar input=points.las n=points_n sum=points_sum \
+ mean=points_mean proportional_n=points_n_prop \
+ proportional_sum=points_sum_prop
+</pre></div>
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r3.in.xyz.html">r3.in.xyz</a>,
+<a href="r.in.lidar.html">r.in.lidar</a>,
+<a href="v.in.lidar.html">v.in.lidar</a>,
+<a href="r.to.rast3.html">r.to.rast3</a>,
+<a href="r3.mapcalc.html">r3.mapcalc</a>,
+<a href="g.region.html">g.region</a>
+</em>
+
+
+<h2>AUTHOR</h2>
+
+Vaclav Petras, <a href="http://geospatial.ncsu.edu/osgeorel/">NCSU OSGeoREL</a>
Property changes on: grass/trunk/raster3d/r3.in.lidar/r3.in.lidar.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass/trunk/raster3d/r3.in.lidar/r3_in_lidar.png
===================================================================
(Binary files differ)
Property changes on: grass/trunk/raster3d/r3.in.lidar/r3_in_lidar.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Modified: grass/trunk/raster3d/raster3dintro.html
===================================================================
--- grass/trunk/raster3d/raster3dintro.html 2016-01-29 09:56:07 UTC (rev 67704)
+++ grass/trunk/raster3d/raster3dintro.html 2016-01-29 21:17:46 UTC (rev 67705)
@@ -85,6 +85,13 @@
2D rasters are considered as slices in this case and
merged into one 3D raster map (<a href="r.to.rast3.html">r.to.rast3</a>).
+<p>
+Import of 3D points and their statistics can be done using
+<a href="r3.in.lidar.html">r3.in.lidar</a> for LiDAR data and
+<a href="r3.in.xyz.html">r3.in.xyz</a> for CSV and other ASCII text
+formats.
+
+
<h3>3D region settings and 3D MASK</h3>
GRASS GIS 3D raster map processing is always performed in the current 3D region
More information about the grass-commit
mailing list