浏览代码

add LAS LiDAR version of r.in.xyz

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@46570 15284696-431f-4ddb-bdfa-cd5b030d7da7
Markus Metz 14 年之前
父节点
当前提交
65ad7ef4ac
共有 5 个文件被更改,包括 1665 次插入0 次删除
  1. 13 0
      raster/r.in.lidar/Makefile
  2. 55 0
      raster/r.in.lidar/local_proto.h
  3. 1192 0
      raster/r.in.lidar/main.c
  4. 275 0
      raster/r.in.lidar/r.in.lidar.html
  5. 130 0
      raster/r.in.lidar/support.c

+ 13 - 0
raster/r.in.lidar/Makefile

@@ -0,0 +1,13 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.in.lidar
+
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB) $(GPROJLIB) $(LASLIBS)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+ifneq ($(USE_LIBLAS),)
+default: cmd
+endif
+

+ 55 - 0
raster/r.in.lidar/local_proto.h

@@ -0,0 +1,55 @@
+/****************************************************************************
+ *
+ * MODULE:       r.in.Lidar
+ *               
+ * AUTHOR(S):    Markus Metz
+ *               Based on r.in.xyz by Hamish Bowman, Volker Wichmann
+ *
+ * PURPOSE:      Imports LAS LiDAR point clouds to a raster map using 
+ *               aggregate statistics.
+ *
+ * COPYRIGHT:    (C) 2011 Markus Metz 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.
+ *
+ *****************************************************************************/
+
+#ifndef __LOCAL_PROTO_H__
+#define __LOCAL_PROTO_H__
+
+
+#include <grass/gis.h>
+#include <liblas/capi/liblas.h>
+
+
+#define BUFFSIZE 256
+
+#define METHOD_N           1
+#define METHOD_MIN         2
+#define METHOD_MAX         3
+#define METHOD_RANGE       4
+#define METHOD_SUM         5
+#define METHOD_MEAN        6
+#define METHOD_STDDEV      7
+#define METHOD_VARIANCE    8
+#define METHOD_COEFF_VAR   9
+#define METHOD_MEDIAN     10
+#define METHOD_PERCENTILE 11
+#define METHOD_SKEWNESS   12
+#define METHOD_TRIMMEAN   13
+
+/* main.c */
+int scan_bounds(LASReaderH, int, double);
+
+/* support.c */
+int blank_array(void *, int, int, RASTER_MAP_TYPE, int);
+int update_n(void *, int, int, int);
+int update_min(void *, int, int, int, RASTER_MAP_TYPE, double);
+int update_max(void *, int, int, int, RASTER_MAP_TYPE, double);
+int update_sum(void *, int, int, int, RASTER_MAP_TYPE, double);
+int update_sumsq(void *, int, int, int, RASTER_MAP_TYPE, double);
+
+
+#endif /* __LOCAL_PROTO_H__ */

文件差异内容过多而无法显示
+ 1192 - 0
raster/r.in.lidar/main.c


+ 275 - 0
raster/r.in.lidar/r.in.lidar.html

@@ -0,0 +1,275 @@
+<h2>DESCRIPTION</h2>
+
+The <em>r.in.lidar</em> module will load and bin LAS LiDAR point clouds
+into a new raster map. The user may choose from a variety of statistical
+methods in creating the new raster. Gridded data provided as a stream of
+x,y,z points may also be imported.
+<p>
+
+<em>r.in.lidar</em> is designed for processing massive point cloud datasets,
+for example raw LIDAR or sidescan sonar swath data. It has been tested with
+datasets as large as tens of billion of points (705GB in a single file).
+ <!-- Doug Newcomb, US Fish & Wildlife Service -->
+<p>
+
+Available statistics for populating the raster are:<br>
+<ul>
+<li>
+<table>
+<tr><td><em>n</em></td>        <td>number of points in cell</td></tr>
+<tr><td><em>min</em></td>      <td>minimum value of points in cell</td></tr>
+<tr><td><em>max</em></td>      <td>maximum value of points in cell</td></tr>
+<tr><td><em>range</em></td>    <td>range of points in cell</td></tr>
+<tr><td><em>sum</em></td>      <td>sum of points in cell</td></tr>
+<tr><td><em>mean</em></td>     <td>average value of points in cell</td></tr>
+<tr><td><em>stddev</em></td>   <td>standard deviation of points in cell</td></tr>
+<tr><td><em>variance</em></td> <td>variance of points in cell</td></tr>
+<tr><td><em>coeff_var</em></td><td>coefficient of variance of points in cell</td></tr>
+<tr><td><em>median</em></td>   <td>median value of points in cell</td></tr>
+<tr valign="baseline"><td><em>percentile</em>&nbsp;</td>
+   <td>p<sup><i>th</i></sup> percentile of points in cell</td></tr>
+<tr><td><em>skewness</em></td> <td>skewness of points in cell</td></tr>
+<tr><td><em>trimmean</em></td> <td>trimmed mean of points in cell</td></tr>
+</table><br>
+
+<li><em>Variance</em> and derivatives use the biased estimator (n). [subject to change]
+<li><em>Coefficient of variance</em> is given in percentage and defined as
+<tt>(stddev/mean)*100</tt>.
+</ul>
+<br>
+
+<h2>NOTES</h2>
+
+<h3>Gridded data</h3>
+
+If data is known to be on a regular grid <em>r.in.lidar</em> can reconstruct
+the map perfectly as long as some care is taken to set up the region
+correctly and that the data's native map projection is used. A typical
+method would involve determining the grid resolution either by examining
+the data's associated documentation or by studying the text file. Next scan
+the data with <em>r.in.lidar</em>'s <b>-s</b> (or <b>-g</b>) flag to find the
+input data's bounds. GRASS uses the cell-center raster convention where
+data points fall within the center of a cell, as opposed to the grid-node
+convention. Therefore you will need to grow the region out by half a cell
+in all directions beyond what the scan found in the file. After the region
+bounds and resolution are set correctly with <em>g.region</em>, run
+<em>r.in.lidar</em> using the <i>n</i> method and verify that n=1 at all places.
+<em>r.univar</em> can help. Once you are confident that the region exactly
+matches the data proceed to run <em>r.in.lidar</em> using one of the <i>mean,
+min, max</i>, or <i>median</i> methods. With n=1 throughout, the result
+should be identical regardless of which of those methods are used.
+
+
+<h3>Memory use</h3>
+
+While the <b>input</b> file can be arbitrarily large, <em>r.in.lidar</em>
+will use a large amount of system memory for large raster regions (10000x10000).
+If the module refuses to start complaining that there isn't enough memory,
+use the <b>percent</b> parameter to run the module in several passes.
+In addition using a less precise map format (<tt>CELL</tt> [integer] or
+<tt>FCELL</tt> [floating point]) will use less memory than a <tt>DCELL</tt>
+[double precision floating point] <b>output</b> map. Methods such as <em>n,
+min, max, sum</em> will also use less memory, while <em>stddev, variance,
+and coeff_var</em> will use more.
+
+The aggregate functions <em>median, percentile, skewness</em> and
+<em>trimmed mean</em> will use even more memory and may not be appropriate
+for use with arbitrarily large input files<!-- without a small value for percent= -->.
+<!-- explained: memory use for regular stats will be based solely on region size,
+ but for the aggregate fns it will also depend on the number of data points. (?) -->
+
+<p>
+
+The default map <b>type</b>=<tt>FCELL</tt> is intended as compromise between
+preserving data precision and limiting system resource consumption.
+If reading data from a <tt>stdin</tt> stream, the program can only run using
+a single pass.
+
+<h3>Setting region bounds and resolution</h3>
+
+You can use the <b>-s</b> scan flag to find the extent of the input data
+(and thus point density) before performing the full import. Use
+<em>g.region</em> to adjust the region bounds to match. The <b>-g</b> shell
+style flag prints the extent suitable as parameters for <em>g.region</em>.
+A suitable resolution can be found by dividing the number of input points
+by the area covered. e.g.
+
+<div class="code"><pre>
+  wc -l inputfile.txt
+  g.region -p
+  # points_per_cell = n_points / (rows * cols)
+
+  g.region -e
+  # UTM location:
+  # points_per_sq_m = n_points / (ns_extent * ew_extent)
+
+  # Lat/Lon location:
+  # points_per_sq_m = n_points / (ns_extent * ew_extent*cos(lat) * (1852*60)^2)
+</pre></div>
+
+<p>
+
+If you only intend to interpolate the data with <em>r.to.vect</em> and
+<em>v.surf.rst</em>, then there is little point to setting the region
+resolution so fine that you only catch one data point per cell -- you might
+as well use "<tt>v.in.ascii&nbsp;-zbt</tt>" directly.
+
+
+<h3>Filtering</h3>
+
+Points falling outside the current region will be skipped. This includes
+points falling <em>exactly</em> on the southern region bound.
+(to capture those adjust the region with "<tt>g.region s=s-0.000001</tt>";
+see <em>g.region</em>)
+<p>
+Blank lines and comment lines starting with the hash symbol (<tt>#</tt>)
+will be skipped.
+
+<p>
+
+The <b>zrange</b> parameter may be used for filtering the input data by
+vertical extent. Example uses might include preparing multiple raster
+sections to be combined into a 3D raster array with <em>r.to.rast3</em>, or
+for filtering outliers on relatively flat terrain.
+
+<p>
+
+In varied terrain the user may find that <em>min</em> maps make for a good
+noise filter as most LIDAR noise is from premature hits. The <em>min</em> map
+may also be useful to find the underlying topography in a forested or urban
+environment if the cells are over sampled.
+
+<p>
+
+The user can use a combination of <em>r.in.lidar</em> <b>output</b> maps to create
+custom filters. e.g. use <em>r.mapcalc</em> to create a <tt>mean-(2*stddev)</tt>
+map. [In this example the user may want to include a lower bound filter in
+<em>r.mapcalc</em> to remove highly variable points (small <em>n</em>) or run
+<em>r.neighbors</em> to smooth the stddev map before further use.]
+
+
+<h3>Interpolation into a DEM</h3>
+
+The vector engine's topological abilities introduce a finite memory overhead
+per vector point which will limit the size a vector map relative to 
+available RAM. If you want more, use the <em>r.to.vect</em>
+<b>-b</b> flag to skip building topology. Without topology, however, all
+you'll be able to do with the vector map is display with <em>d.vect</em> and
+interpolate with <em>v.surf.rst</em>.
+Run <em>r.univar</em> on your raster map to check the number of non-NULL cells
+and adjust bounds and/or resolution as needed before proceeding.
+
+<p>
+
+Typical commands to create a DEM using a regularized spline fit:
+<div class="code"><pre>
+  r.univar lidar_min
+  r.to.vect -z feature=point in=lidar_min out=lidar_min_pt
+  v.surf.rst layer=0 in=lidar_min_pt elev=lidar_min.rst
+</pre></div>
+
+Typical commands to create a DEM using bsplines:
+<div class="code"><pre>
+  r.resamp.bspline in=lidar_min out=lidar_min.bspline
+</pre></div>
+<br>
+
+<h2>EXAMPLE</h2>
+
+Import the <a href="http://www.grassbook.org/data_menu2nd.phtml">Jockey's
+Ridge, NC, LIDAR dataset</a>, and process into a clean DEM:
+
+<div class="code"><pre>
+    # scan and set region bounds
+  r.in.xyz -s fs=, in=lidaratm2.txt out=test
+  g.region n=35.969493 s=35.949693 e=-75.620999 w=-75.639999
+  g.region res=0:00:00.075 -a
+    # create "n" map containing count of points per cell for checking density
+  r.in.xyz in=lidaratm2.txt out=lidar_n fs=, method=n zrange=-2,50
+    # check point density [rho = n_sum / (rows*cols)]
+  r.univar lidar_n | grep sum
+    # create "min" map (elevation filtered for premature hits)
+  r.in.xyz in=lidaratm2.txt out=lidar_min fs=, method=min zrange=-2,50
+    # zoom to area of interest
+  g.region n=35:57:56.25N s=35:57:13.575N w=75:38:23.7W e=75:37:15.675W
+    # check number of non-null cells (try and keep under a few million)
+  r.univar lidar_min | grep '^n:'
+    # convert to points 
+  r.to.vect -z feature=point in=lidar_min out=lidar_min_pt
+    # interpolate using a regularized spline fit
+  v.surf.rst layer=0 in=lidar_min_pt elev=lidar_min.rst
+    # set color scale to something interesting
+  r.colors lidar_min.rst rule=bcyr -n -e
+    # prepare a 1:1:1 scaled version for NVIZ visualization (for lat/lon input)
+  r.mapcalc "lidar_min.rst_scaled = lidar_min.rst / (1852*60)"
+  r.colors lidar_min.rst_scaled rule=bcyr -n -e
+</pre></div>
+<br>
+
+
+<h2>TODO</h2>
+
+<ul>
+<li> Support for multiple map output from a single run.<br>
+     <tt>method=string[,string,...] output=name[,name,...]</tt>
+<li> Merge with r.in.xyz.
+     <!-- Bob Covill has supplied patches for MBIO interface already -->
+</ul>
+
+<h2>BUGS</h2>
+
+<ul>
+<li> <em>n</em> map sum can be ever-so-slightly more than `<tt>wc -l</tt>`
+  with e.g. <tt>percent=10</tt> or less.
+  <br>Cause unknown.
+
+<li> <em>n</em> map <tt>percent=100</tt> and <tt>percent=xx</tt> maps
+  differ slightly (point will fall above/below the segmentation line)
+  <br>Investigate with "<tt>r.mapcalc diff = bin_n.100 - bin_n.33</tt>" etc.
+  <br>Cause unknown.
+
+<li> "<tt>nan</tt>" can leak into <em>coeff_var</em> maps.
+  <br>Cause unknown. Possible work-around: "<tt>r.null setnull=nan</tt>"
+<!-- Another method:  r.mapcalc 'No_nan = if(map == map, map, null() )' -->
+</ul>
+
+If you encounter any problems (or solutions!) please contact the GRASS
+Development Team.
+
+<h2>SEE ALSO</h2>
+
+<i>
+<a href="g.region.html">g.region</a><br>
+<a href="m.proj.html">m.proj</a><br>
+<a href="r.in.xyz.html">r.in.xyz</a><br>
+<a href="r.fillnulls.html">r.fillnulls</a><br>
+<a href="r.in.ascii.html">r.in.ascii</a><br>
+<a href="r.mapcalc.html">r.mapcalc</a><br>
+<a href="r.neighbors.html">r.neighbors</a><br>
+<a href="r.out.xyz.html">r.out.xyz</a><br>
+<a href="r.to.rast3.html">r.to.rast3</a><br>
+<a href="r.to.vect.html">r.to.vect</a><br>
+<a href="r.univar.html">r.univar</a><br>
+<a href="v.in.ascii.html">v.in.ascii</a><br>
+<a href="v.surf.rst.html">v.surf.rst</a><br>
+<br>
+<a href="v.lidar.correction.html">v.lidar.correction</a>,
+<a href="v.lidar.edgedetection.html">v.lidar.edgedetection</a>,
+<a href="v.lidar.growing.html">v.lidar.growing</a>,
+<a href="v.outlier.html">v.outlier</a>,
+<a href="v.surf.bspline.html">v.surf.bspline</a>
+</i>
+<p>
+<i><a href="http://www.ivarch.com/programs/pv.shtml">pv</a></i>
+ - The UNIX pipe viewer utility
+<BR><BR>
+
+
+<h2>AUTHORS</h2>
+
+Markus Metz<br>
+based on r.in.xyz by Hamish Bowman and Volker Wichmann<br>
+
+<br>
+<p>
+<i>Last changed: $Date$</i></p>

+ 130 - 0
raster/r.in.lidar/support.c

@@ -0,0 +1,130 @@
+/*
+ * r.in.lidar support fns, from r.in.xyz.
+ *   Copyright 2006 by M. Hamish Bowman, and The GRASS Development Team
+ *   Author: M. Hamish Bowman, University of Otago, Dunedin, New Zealand
+ *
+ *   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/raster.h>
+#include "local_proto.h"
+
+static void *get_cell_ptr(void *array, int cols, int row, int col,
+			  RASTER_MAP_TYPE map_type)
+{
+    return G_incr_void_ptr(array,
+			   ((row * (size_t) cols) +
+			    col) * Rast_cell_size(map_type));
+}
+
+int blank_array(void *array, int nrows, int ncols, RASTER_MAP_TYPE map_type,
+		int value)
+{
+    /* flood fill initialize the array to either 0 or NULL */
+    /*  "value" can be either 0 (for 0.0) or -1 (for NULL) */
+    int row, col;
+    void *ptr;
+
+    ptr = array;
+
+    switch (value) {
+    case 0:
+	/* fill with 0 */
+	/* simpler to use Rast_raster_cpy() or similar ?? */
+
+	for (row = 0; row < nrows; row++) {
+	    for (col = 0; col < ncols; col++) {
+		Rast_set_c_value(ptr, 0, map_type);
+		ptr = G_incr_void_ptr(ptr, Rast_cell_size(map_type));
+	    }
+	}
+	break;
+
+    case -1:
+	/* fill with NULL */
+	/* alloc for col+1, do we come up (nrows) short? no. */
+	Rast_set_null_value(array, nrows * ncols, map_type);
+	break;
+
+    default:
+	return -1;
+    }
+
+    return 0;
+}
+
+/* make all these fns return void? */
+int update_n(void *array, int cols, int row, int col)
+{
+    void *ptr = get_cell_ptr(array, cols, row, col, CELL_TYPE);
+    CELL old_n;
+
+    old_n = Rast_get_c_value(ptr, CELL_TYPE);
+    Rast_set_c_value(ptr, (1 + old_n), CELL_TYPE);
+
+    return 0;
+}
+
+
+int update_min(void *array, int cols, int row, int col,
+	       RASTER_MAP_TYPE map_type, double value)
+{
+    void *ptr = get_cell_ptr(array, cols, row, col, map_type);
+    DCELL old_val;
+
+    if (Rast_is_null_value(ptr, map_type))
+	Rast_set_d_value(ptr, (DCELL) value, map_type);
+    else {
+	old_val = Rast_get_d_value(ptr, map_type);
+	if (value < old_val)
+	    Rast_set_d_value(ptr, (DCELL) value, map_type);
+    }
+    return 0;
+}
+
+
+int update_max(void *array, int cols, int row, int col,
+	       RASTER_MAP_TYPE map_type, double value)
+{
+    void *ptr = get_cell_ptr(array, cols, row, col, map_type);
+    DCELL old_val;
+
+    if (Rast_is_null_value(ptr, map_type))
+	Rast_set_d_value(ptr, (DCELL) value, map_type);
+    else {
+	old_val = Rast_get_d_value(ptr, map_type);
+	if (value > old_val)
+	    Rast_set_d_value(ptr, (DCELL) value, map_type);
+    }
+
+    return 0;
+}
+
+
+int update_sum(void *array, int cols, int row, int col,
+	       RASTER_MAP_TYPE map_type, double value)
+{
+    void *ptr = get_cell_ptr(array, cols, row, col, map_type);
+    DCELL old_val;
+
+    old_val = Rast_get_d_value(ptr, map_type);
+    Rast_set_d_value(ptr, value + old_val, map_type);
+
+    return 0;
+}
+
+
+int update_sumsq(void *array, int cols, int row, int col,
+		 RASTER_MAP_TYPE map_type, double value)
+{
+    void *ptr = get_cell_ptr(array, cols, row, col, map_type);
+    DCELL old_val;
+
+    old_val = Rast_get_d_value(ptr, map_type);
+    Rast_set_d_value(ptr, (value * value) + old_val, map_type);
+
+    return 0;
+}