r.in.lidar: more detailed and extented docs

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

raster/r.in.lidar/main.c

@@ -106,6 +106,10 @@ int main(int argc, char *argv[])
     G_add_keyword(_("raster"));
     G_add_keyword(_("import"));
     G_add_keyword(_("LIDAR"));
+    G_add_keyword(_("statistics"));
+    G_add_keyword(_("conversion"));
+    G_add_keyword(_("aggregation"));
+    G_add_keyword(_("binning"));
     module->description =
 	_("Creates a raster map from LAS LiDAR points using univariate statistics.");
 
@@ -135,6 +139,33 @@ int main(int argc, char *argv[])
 	"n,min,max,range,sum,mean,stddev,variance,coeff_var,median,percentile,skewness,trimmean";
     method_opt->answer = "mean";
     method_opt->guisection = _("Statistic");
+    G_asprintf((char **)&(method_opt->descriptions),
+               "n;%s;"
+               "min;%s;"
+               "max;%s;"
+               "range;%s;"
+               "sum;%s;"
+               "mean;%s;"
+               "stddev;%s;"
+               "variance;%s;"
+               "coeff_var;%s;"
+               "median;%s;"
+               "percentile;%s;"
+               "skewness;%s;"
+               "trimmean;%s",
+               _("Number of points in cell"),
+               _("Minimum value of point values in cell"),
+               _("Maximum value of point values in cell"),
+               _("Range of point values in cell"),
+               _("Sum of point values in cell"),
+               _("Mean (average) value of point values in cell"),
+               _("Standard deviation of point values in cell"),
+               _("Variance of point values in cell"),
+               _("Coefficient of variance of point values in cell"),
+               _("Median value of point values in cell"),
+               _("pth (nth) percentile of point values in cell"),
+               _("Skewness of point values in cell"),
+               _("Trimmed mean of point values in cell"));
 
     type_opt = G_define_option();
     type_opt->key = "type";
@@ -147,8 +178,11 @@ int main(int argc, char *argv[])
     base_raster_opt = G_define_standard_option(G_OPT_R_INPUT);
     base_raster_opt->key = "base_raster";
     base_raster_opt->required = NO;
-    base_raster_opt->label = _("Subtract raster values from the z coordinates");
-    base_raster_opt->description = _("The scale for z is applied beforehand, the filter afterwards");
+    base_raster_opt->label =
+        _("Subtract raster values from the Z coordinates");
+    base_raster_opt->description =
+        _("The scale for Z is applied beforehand, the range filter for"
+          " Z afterwards");
     base_raster_opt->guisection = _("Transform");
 
     zrange_opt = G_define_option();
@@ -156,7 +190,7 @@ int main(int argc, char *argv[])
     zrange_opt->type = TYPE_DOUBLE;
     zrange_opt->required = NO;
     zrange_opt->key_desc = "min,max";
-    zrange_opt->description = _("Filter range for z data (min,max)");
+    zrange_opt->description = _("Filter range for Z data (min,max)");
     zrange_opt->guisection = _("Selection");
 
     zscale_opt = G_define_option();
@@ -164,7 +198,7 @@ int main(int argc, char *argv[])
     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->description = _("Scale to apply to Z data");
     zscale_opt->guisection = _("Transform");
 
     irange_opt = G_define_option();
@@ -246,8 +280,11 @@ int main(int argc, char *argv[])
 
     extents_flag = G_define_flag();
     extents_flag->key = 'e';
+    extents_flag->label =
+        _("Use the extent of the input for the raster extent");
     extents_flag->description =
-	_("Extend region extents based on new dataset");
+        _("Set internally computational region extents based on the"
+          " point cloud");
     extents_flag->guisection = _("Output");
 
     set_region_flag = G_define_flag();
@@ -278,18 +315,24 @@ int main(int argc, char *argv[])
 
     intens_flag = G_define_flag();
     intens_flag->key = 'i';
+    intens_flag->label =
+        _("Use intensity values rather than Z values");
     intens_flag->description =
-        _("Import intensity values rather than z values");
+        _("Uses intensity values everywhere as if they would be Z"
+          " coordinates");
 
     intens_import_flag = G_define_flag();
     intens_import_flag->key = 'j';
     intens_import_flag->description =
-        _("Use z values for filtering, but import intensity values");
+        _("Use Z values for filtering, but intensity values for statistics");
 
     base_rast_res_flag = G_define_flag();
     base_rast_res_flag->key = 'd';
+    base_rast_res_flag->label =
+        _("Use base raster resolution instead of computational region");
     base_rast_res_flag->description =
-        _("Use base raster actual resolution instead of computational region");
+        _("For getting values from base raster, use its actual"
+          " resolution instead of computational region resolution");
 
     G_option_exclusive(intens_flag, intens_import_flag, NULL);
 
@@ -710,7 +753,7 @@ int main(int argc, char *argv[])
     /* close raster file & write history */
     Rast_close(out_fd);
 
-    sprintf(title, "Raw x,y,z data binned into a raster grid by cell %s",
+    sprintf(title, "Raw X,Y,Z data binned into a raster grid by cell %s",
             method_opt->answer);
     Rast_put_cell_title(outmap, title);
 

raster/r.in.lidar/r.in.lidar.html

@@ -1,21 +1,21 @@
 <h2>DESCRIPTION</h2>
 
-The <em>r.in.lidar</em> module loads and bins LAS LiDAR point clouds
-into a new raster map. The user may choose from a variety of statistical
-methods in creating the new raster.
+The <em>r.in.lidar</em> module loads LAS LiDAR point clouds into a new
+raster map using binning. The user may choose from a variety of
+statistical methods which will be used for binning when creating
+the new raster.
 
 <p>
-Since the creation of raster maps depends on the computational 
-region settings (extent and resolution), as default the current 
-region extents and resolution are used for the import. When using 
-the <em>-e</em> flag along with the <em>resolution=value</em> 
-parameter, the region extents will be based on new dataset. It is therefore
-recommended to first use the <em>-s</em> flag to get the extents of the
-LiDAR point cloud to be imported, then adjust the current region extent
-and resolution accordingly, and only then proceed with the actual import.
-Another option is to automatically set the region extents based on the
-LAS dataset itself along with the desired raster resolution. See below
-for details.
+Since a new raster map is created during the binning, the binning of
+points depends on the current computational region settings
+(extent and resolution) by default (see more about binning below).
+When using the <b>-e</b> flag, the binning will be done in the extent
+of the point cloud, so the resulting raster will have extent based on
+the input point cloud.
+When the <em>resolution=value</em> parameter is used,
+the binning is done using the provided resolution and the resulting
+raster will have that resolution (see more below for more information
+about extent and resolution management).
 
 <p>
 <em>r.in.lidar</em> is designed for processing massive point cloud 
@@ -23,51 +23,311 @@ datasets, for example raw LiDAR or sidescan sonar swath data. It has
 been tested with large datasets (see below for memory management 
 notes).
 
+<h3>Binning</h3>
+
+The main different of <em>r.in.lidar</em> in comparison to
+<em><a href="r.in.lidar.html">r.in.lidar</a></em> is that
+<em>r.in.lidar</em> creates a raster instead of just importing the
+points into GRASS GIS. However, <em>r.in.lidar</em> does not merely
+rasterizes the points from the point cloud. <em>r.in.lidar</em>
+uses binning to derive values for individual raster cells,
+so the value of a cell is typically an aggregation of values
+of individual points falling into one cell.
+
+In general binning is the conversion of points into a regular grid. The bin-
+ning of points with X and Y coordinates starts with the overlay of
+a grid of bins over the points.
 <p>
-Available statistics for populating the output raster map 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>
+In the basic case, binning is a method which counts the number of
+points which fall into one raster cell, i.e. bin. The number of points
+per cell (bin) indicates the density of points in the point cloud.
+The cell (bin) is always square or rectangular in case of
+<em>r.in.lidar</em> because the result is GRASS GIS 2D raster.
+The result of binning where the number of point per cell is counted
+is sometimes called 2D (two dimensional) histogram because
+a histogram is used in univariate statistics (in one dimension)
+to count the number samples falling into a given bin.
 
-<h2>NOTES</h2>
+<center>
+<img src="r_in_lidar_binning_count.png">
+<img src="r_in_lidar_binning_mean.png">
+<p><em>
+    Figure: The binning on left was used to count number of points per
+    (sometimes also called 2D histogram). The numbers in cells are
+    examples of counts, the rest is represented by the color.
+    The binning on right was used with mean to create a surface
+    based on the values associated with the points. The numbers
+    show examples of cell values. Note also the cells without any points
+    which were assigned the NULL value.
+</em>
+</center>
 
-<h3>LAS file import preparations</h3>
+The basic concept of binning is extended when the points have another
+value associated with them. For LiDAR data this value can be the Z
+coordinate or intensity. The value for a given cell (bin) is computed
+using univariate statistics from the values of all points in the cell.
+For example, computing the mean value of Z coordinates can yield
+a raster representing the digital elevation model. Another example is
+the range of Z coordinates which can be used as a rough estimate of
+vegetation height.
+
+<h3>Statistics</h3>
+
+Available statistics for populating the output raster map are:
+
+<dl class="option_descriptions">
+<dt>n</dt>
+<dd>This computes the number (count) of points per cell. The result
+is a indicator of spatially variable density of points in the given
+area.</dd>
+<dt>min</dt>
+<dd>This finds the minimum of point values in each cell.
+It can be useful when finding topography in a forested or urban
+environment and there is a lot of points per one cells (terrain is
+oversampled considering the desired resolution).
+It can also create surfaces independent on the noise from premature
+hits as it will always select the lowest point.
+</dd>
+<dt>max</dt>
+<dd>This finds the maximum of point values in each cell.
+In connection with <b>base_raster</b> it can yield maximum vegetation
+of feature height per cell.
+For this purpose, it is usually much more appropriate than <em>mean</em>
+which would yield heights mostly influenced by the vertical
+distribution of points.
+</dd>
+<dt>range</dt>
+<dd>This computes the range of point values in each cell.
+The range of Z coordinates per cell can be used as a rough estimate of
+vegetation height when the cells are small enough, slopes low
+and the area is mostly vegetated.
+However, for more profound analysis, the base raster together with
+different statistics is recommended.</dd>
+<dt>sum</dt>
+<dd>This computes the sum of point values per cell.
+This is useful especially when intensity is used as a value
+(flags <b>-i</b> and <b>-j</b>).</dd>
+<dt>mean</dt>
+<dd>This is a mean (average) value of point values in cell.
+When used with Z coordinates (the default) and points from the ground
+class, the resulting raster is a digital elevation model.
+When intensity is used as a point value, the resulting raster contains
+mean intensity per cell.
+Note that <em>mean</em> gives heights influenced by the vertical
+distribution of points</dd>
+<dt>stddev</dt>
+<dd>This computes the standard deviation of point values for each
+cell.</dd>
+<dt>variance</dt>
+<dd>This computes the variance of point values for each cell.
+Variance and derivatives use the biased estimator (n)
+[note that this might be subject to change].</dd>
+<dt>coeff_var</dt>
+<dd>This computes the coefficient of variance of point values for each
+cell. Coefficient of variance is given in percentage and defined as
+<tt>(stddev/mean)*100</tt>.</dd>
+<dt>median</dt>
+<dd>This computes the median of point values for each cell</dd>
+<dt>percentile</dt>
+<dd>p<sup><i>th</i></sup> (nth) percentile of points in cell</dd>
+<dt>skewness</dt>
+<dd>This is a skewness of point values in cell</dd>
+<dt>trimmean</dt>
+<dd>This is a trimmed mean of point values in cell.
+Trimmed mean also know as truncated mean is a mean
+computed after discarding values at the low end and at the high end.
+How many values to discard is given by the <b>trim</b> option
+in percent. In statistics the usual percentage of trimmed values ranges
+from 5 to 25 percent.</dd>
+</dl>
+
+Note that different statistics have different memory requirements
+(see below for details).
+
+<h3>Filtering</h3>
+
+Points falling outside the current computational region will be skipped.
+This includes points falling <em>exactly</em> on the southern region
+bound. To capture those adjust the region with:
+
+<div class="code"><pre>
+g.region s=s-0.000001
+</pre></div>
+
+See <em><a href="g.region.html">g.region</a></em> for details about
+computation region handling in GRASS GIS.
+
+<p>
+The <b>zrange</b> parameter may be used for filtering the input data by
+vertical extent. Example uses include
+filtering out extreme outliers and outliers on relatively flat terrain.
+This parameter can be also used for cutting the point cloud into
+vertical sections preparing it for further processing
+by separate sections, together as if it would be an imagery group
+(see <em><a href="i.group.html">i.group</a></em>), or combined into
+a 3D raster using <em><a href="r.to.rast3.html">r.to.rast3</a></em>.
+In for these last examples, it might actually be more advantageous
+to use <em><a href="r3.in.lidar.html">r3.in.lidar</a></em> module.
+The <b>zrange</b> parameter is especially powerful when used
+together with the <b>base_raster</b> parameter. The <b>zrange</b>
+is applied to Z values after the <b>base_raster</b> reduction.
+
+<center>
+<img src="r_in_lidar_zrange.png">
+<p><em>
+    Figure: This is the principle of zrange filter. Points with the
+    Z coordinate value below the lower value in the range (here 180)
+    are filtered out (blue points) and same applies for points above
+    higher value in the range (here 250). All other points are preserved
+    (green points).
+</em>
+</center>
+
+<p>
+A LiDAR pulse can have multiple returns. The first return values can be 
+used to obtain a digital surface model (DSM) where e.g. canopy cover is 
+represented. The last return values can be used to obtain a digital 
+terrain model (DTM) where e.g. the forest floor instead of canopy 
+cover is represented. The <b>return_filter</b> option allows selecting
+one of first, mid, or last returns. Return number and number of returns
+in the pulse associated with each point are compared to determine
+if the point is first, mid, or last return.
+
+<p>
+LiDAR points often come as already classified into standardized classes.
+For example, class number 2 represents ground. For other classes see
+LAS format specification in references. The <b>class_filter</b> option
+allows selecting one or more classes using numbers (integers) separated
+by comma.
+
+<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 oversampled.
+
+<p>
+The user can use a combination of <em>r.in.lidar</em> <b>output</b> maps
+to create custom raster-based filters, for examplee, use
+<em><a href="r.mapcalc.html">r.mapcalc</a></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.)
+
+<p>
+Note that proper filtering of the input points in not only critical for
+the analysis itself but it can also speed up the processing.
+
+<h3>Reduction to a base raster</h3>
+
+For analysis of features on the terrain surface, especially vegetation
+it is advantageous to remove the influence of the terrain on heights
+because the height above the terrain is important (e.g. height of
+a tree) rather than height of the top of the tree above the see level.
+In this case, the base raster would be digital elevation model
+which can be one derived from the point cloud, or obtained in
+some other way. LiDAR data often come with precomputed DEMs
+(quality should be checked in this case) and there is often a DEM
+available for a given area (fit with the point cloud, especially
+vertical, and resolution should be checked).
+
+<center>
+<img src="r_in_lidar_base_raster.png">
+<p><em>
+    Figure: This is a profile of base raster (in orange) representing
+    digital elevation model and selected points, e.g. first return,
+    from point cloud (green dots). By default the points would create
+    a digital surface model (thin brown line) but after reducing the
+    Z coordinates using the base raster, the created surface is a
+    derived from the height of points relative to the base raster.
+</em>
+</center>
+
+The usage of base raster is not limited to digital elevation model.
+The base raster can be any surface which has some relation to the
+point values, for example digital surface model representing
+top of the canopy.
+
+<h3>Setting extent and resolution</h3>
+
+<p>
+Since the creation of raster maps depends on the computational 
+region settings (extent and resolution), as default the current 
+region extents and resolution are used for the import. When using 
+the <em>-e</em> flag along with the <em>resolution=value</em> 
+parameter, the region used for the new raster will be based
+the point cloud extent and the provided resolution. It is therefore
+recommended to first use the <em>-s</em> flag to get the extents of the
+LiDAR point cloud to be imported, then adjust the current region extent
+and resolution accordingly, and only then proceed with the actual import.
+Another option is to automatically set the region extents based on the
+LAS dataset itself (<em>-e</em> flag) along with the desired raster
+resolution. The best option is to know the point cloud extent ahead,
+e.g. from tiling scheme, and use it. See below for details.
 
+<p>
 Since the <em>r.in.lidar</em> generates a raster map through binning 
 from the original LiDAR points, the target computational region 
 extent and resolution have to be determined. A typical workflow 
 would involve the examination of the LAS data's associated 
 documentation or the scan of the LAS data file with
 <em>r.in.lidar</em>'s <b>-s</b> (or <b>-g</b>) flag to find the input
-data's bounds.<br>
+data's bounds.
+
+<p>
 Another option is to automatically set the region extents based on the
 LAS dataset extent (<b>-e</b> flag) along with the desired raster
 resolution using the <em>resolution</em> parameter.
 
+<p>
+Using the <b>-s</b> scan flag, the extent of the input data (and thus
+point density) is printed. To check this is recommended before performing
+the full import. The <b>-g</b> shell style flag prints the extent suitable
+as command line parameters for <em>g.region</em>.
+
+<p>
+A simpler option is to automatically set the region extents based on the
+LAS dataset (<b>-e</b> flag) along with the target raster resolution using
+the <em>resolution</em> parameter. Also here it is recommended to verify
+and optimize the resulting region settings with <em>g.region</em> prior
+to importing the dataset.
+
+
+<h2>NOTES</h2>
+
+<h3>Format and projection support</h3>
+
+The typical file extensions for the LAS format are .las and .laz
+(compressed). The compressed LAS (.laz) format can be imported only if
+libLAS has been compiled with LASzip support. It is also recommended to
+compile libLAS with GDAL which is used to test if the LAS projection
+matches that of the GRASS location.
+
+<h3>LAS file import preparations</h3>
+
+Note that the scanning (<b>-s</b> or <b>-g</b> flags) needs to iterate
+over the whole point cloud. This will take a long time for large
+datasets, so if the user knows the approximate extent of the dataset,
+for example because it dataset for one county or tiling scheme is
+available as vector polygons, it is much more advantageous to provide
+the extent information instead of retrieving it from the dataset.
+The same applies to the <b>-e</b> flag which also needs to perform
+scanning before the binning begins.
+
+<p>
+Also note that the scanning does not apply any filters, so the
+extent determined by scanning can be theoretically bigger than
+the extent actively used during the binning.
+This behavior ensures that the newly created raster has always
+the same extent regardless the filters used.
+However, for most cases (considering the point cloud and the resolution
+used) there is no difference between the extent without filters applied
+and the extent if the filters would be applied.
 
-<h3>Memory use</h3>
+<h3>Memory consumption</h3>
 
+<p>
 While the <b>input</b> file can be arbitrarily large, <em>r.in.lidar</em>
 will use a large amount of system memory (RAM) for large raster regions
 (&gt; 10000x10000 pixels).
@@ -75,37 +335,52 @@ 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. (?) -->
+[double precision floating point] <b>output</b> map.
+For <b>method</b>=<em>n</em>, the <tt>CELL</tt> format is used
+automatically.
+
 <p>
-A LiDAR pulse can have multiple returns. The first return values can be 
-used to obtain a digital surface model (DSM) where e.g. canopy cover is 
-represented. The last return values can be used to obtain a digital 
-terrain model (DTM) where e.g. the forest floor instead of canopy 
-cover is represented. The <b>return_filter</b> option allows selecting one of 
-first, mid, or last returns.
+The <em>mean</em> and <em>range</em> methods will use average amount
+of memory (comparing to other methods).
+Methods such as <em>n, min, max</em>, and <em>sum</em> will use
+less memory,
+while <em>stddev, variance</em>, and <em>coeff_var</em> will use more.
+
 <p>
-LiDAR points can be already classified into standardized classes. For example,
-class number 2 represents ground (for other classes see LAS format specification
-in references). The <b>class_filter</b> option allows selecting one or more 
-classes, as numbers (integers) separated by comma.
+The memory usage for regular statistics mentioned above is based solely
+on region (raster) size.
+However, the aggregate functions <em>median, percentile, skewness</em>
+and <em>trimmean</em> will use more memory and may not be
+appropriate for use with arbitrarily large input files without
+a small value for the <b>percent</b> option because unlike
+the other statistics memory use for these also depends 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.
 
-<h3>Setting region bounds and resolution</h3>
+<h2>EXAMPLES</h2>
+
+Simple example of binning of point from a LAS file into a newly created
+raster map in an existing location/mapset (using metric units):
+
+<div class="code"><pre>
+# set the computational region automatically, resol. for binning is 5m
+r.in.lidar -e -o input=points.las resolution=5 output=lidar_dem_mean
+g.region raster=lidar_dem_mean -p
+r.univar lidar_dem_mean
+</pre></div>
+
+<h3>Finding suitable extent and resolution</h3>
 
+<!-- TODO: there is duplication with the text in the description -->
 Using the <b>-s</b> scan flag, the extent of the input data (and thus
 point density) is printed. To check this is recommended before performing
 the full import. The <b>-g</b> shell style flag prints the extent suitable
-as command line parameters for <em>g.region</em>.<br>
+as command line parameters for <em>g.region</em>.
+
+<p>
 A simpler option is to automatically set the region extents based on the
 LAS dataset (<b>-e</b> flag) along with the target raster resolution using
 the <em>resolution</em> parameter. Also here it is recommended to verify
@@ -171,52 +446,6 @@ g.region -e
   # points_per_sq_m = n_points / (ns_extent * ew_extent*cos(lat) * (1852*60)^2)
 </pre></div>
 
-<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 oversampled.
-
-<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.]
-
-<h2>NOTES</h2>
-
-The typical file extensions for the LAS format are .las and .laz (compressed). 
-The compressed LAS (.laz) format can be imported only if libLAS has been compiled 
-with laszip support. It is also recommended to compile libLAS with GDAL which is
-used to test if the LAS projection matches that of the GRASS location.
-
-<h2>EXAMPLES</h2>
-
-Import of a LAS file into an existing location/mapset (metric units):
-
-<div class="code"><pre>
-# set the computational region automatically, resol. for binning is 5m
-r.in.lidar -e -o input=points.las resolution=5 output=lidar_dem_mean
-g.region raster=lidar_dem_mean -p
-r.univar lidar_dem_mean
-</pre></div>
-
 <h3>Serpent Mound dataset</h3>
 
 This example is analogous to the example used in the GRASS wiki page for
@@ -265,6 +494,7 @@ r.in.lidar input=points.las output=mean_height_above_ground base_raster=elevatio
 In this type of computation, it might be advantageous to change the resolution
 to match the precision of the points rather than deriving it from the base raster.
 <!-- TODO: say how -->
+
 <h3>Multiple file input</h3>
 
 The file option requres a file that contains a list of file names with the full 
@@ -297,29 +527,45 @@ g.region raster=elevation -p
 r.in.lidar file=/home/user/data/filelist.txt output=mean_height_above_ground base_raster=elevation method=mean
 </pre></div>
 
+In Python, the list of files can be created using the <em>glob</em>
+Python module:
+
+<div class="code"><pre>
+import glob
+import gscript
+    
+file_list_name = '/home/user/data/filelist.txt'
+with open(, mode='w') as file_list:
+    for path in glob.iglob('/home/user/data/lidar/*.las'):
+        file_list.write(path + "\n")
+
+gscript.run_command('r.in.lidar', file=file_list_name,
+                    output='mean_height_above_ground',
+                    base_raster='elevation' method='mean')
+</pre></div>
+
 
-<h2>TODO</h2>
+<h2>KNOWN ISSUES</h2>
 
 <ul>
-<li> Support for multiple map output from a single run.<br>
-     <tt>method=string[,string,...] output=name[,name,...]</tt>
+<li>The "<tt>nan</tt>" value (as defined in C language) can leak into
+    <em>coeff_var</em> raster maps. Cause is unknown. Possible
+    work-around is: <tt>r.null setnull=nan</tt> or
+    <tt>r.mapcalc 'no_nan = if(map == map, map, null())'</tt>.
+<li>Only one method can be applied for a single run and multiple map
+    output from a single run
+    (e.g. <tt>method=string[,string,...] output=name[,name,...]</tt>
+    or <tt>n=string mean=string</tt>) is no supported.
 <!-- not really:
 <li> Merge with r.in.xyz.
 -->
      <!-- Bob Covill has supplied patches for MBIO interface already -->
 </ul>
 
-<h2>KNOWN ISSUES</h2>
-
-<ul>
-<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>
 
 <em>
@@ -335,18 +581,37 @@ Development Team.
 <a href="v.outlier.html">v.outlier</a>,
 <a href="v.surf.bspline.html">v.surf.bspline</a>
 </em>
+<br>
+<a href="https://en.wikipedia.org/wiki/Truncated_mean">Trimmed mean</a>
+(Truncated mean, Wikipedia article),
+<a href="http://opentopography.org/">OpenTopography</a>
+(LiDAR point cloud repository)
+
 
 <h2>REFERENCES</h2>
 
+<ul>
+<li>
+<!-- TODO: link the PDF -->
+V. Petras, A. Petrasova, J. Jeziorska, H. Mitasova (2016):
+<em>Processing UAV and lidar point clouds in GRASS GIS</em>.
+XXIII ISPRS Congress 2016 [in press]
+<li>
 <a href="http://www.asprs.org/Committee-General/LASer-LAS-File-Format-Exchange-Activities.html">
-ASPRS LAS format</a><br>
-<a href="http://www.liblas.org/">LAS library</a> <br>
+ASPRS LAS format</a>
+<li>
+<a href="http://www.liblas.org/">LAS library</a>
+<li>
 <a href="http://test.liblas.org/doxygen/liblas_8h.htm">LAS library C API</a> documentation
+</ul>
 
 <h2>AUTHORS</h2>
 
 Markus Metz<br>
+Vaclav Petras,
+<a href="http://geospatial.ncsu.edu/osgeorel/">NCSU OSGeoREL</a>
+(base_raster option, documentation)
+<br>
 based on r.in.xyz by Hamish Bowman and Volker Wichmann<br>
 
-<br>
 <p><i>Last changed: $Date$</i>

raster/r.in.lidar/r_in_lidar_zrange.svg

@@ -0,0 +1,298 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<svg
+   xmlns:dc="http://purl.org/dc/elements/1.1/"
+   xmlns:cc="http://creativecommons.org/ns#"
+   xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
+   xmlns:svg="http://www.w3.org/2000/svg"
+   xmlns="http://www.w3.org/2000/svg"
+   version="1.1"
+   id="svg2"
+   viewBox="0 0 504.05695 354.90146"
+   height="100.16108mm"
+   width="142.25607mm">
+  <defs
+     id="defs4" />
+  <metadata
+     id="metadata7">
+    <rdf:RDF>
+      <cc:Work
+         rdf:about="">
+        <dc:format>image/svg+xml</dc:format>
+        <dc:type
+           rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
+        <dc:title></dc:title>
+      </cc:Work>
+    </rdf:RDF>
+  </metadata>
+  <g
+     transform="translate(-137.97152,-281.58319)"
+     id="layer1">
+    <circle
+       r="4.1224623"
+       cy="423.7908"
+       cx="227.14287"
+       id="path4136"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="479.50507"
+       cx="255.71428"
+       id="path4136-0"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="520.21936"
+       cx="175"
+       id="path4136-8"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="415.93362"
+       cx="250.71428"
+       id="path4136-3"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="435.21936"
+       cx="282.85715"
+       id="path4136-1"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="420.21936"
+       cx="316.42856"
+       id="path4136-4"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="403.07648"
+       cx="321.42856"
+       id="path4136-6"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="385.93362"
+       cx="332.85715"
+       id="path4136-10"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="375.21936"
+       cx="362.14285"
+       id="path4136-41"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="379.50507"
+       cx="370.71429"
+       id="path4136-56"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="391.64792"
+       cx="392.14285"
+       id="path4136-14"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="403.07648"
+       cx="400.71429"
+       id="path4136-7"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="421.64792"
+       cx="436.42856"
+       id="path4136-49"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="406.64792"
+       cx="445"
+       id="path4136-55"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="405.21936"
+       cx="472.85715"
+       id="path4136-2"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="426.64792"
+       cx="475"
+       id="path4136-00"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="458.79077"
+       cx="480.71429"
+       id="path4136-54"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="466.64792"
+       cx="485.71429"
+       id="path4136-16"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="477.36221"
+       cx="503.57144"
+       id="path4136-33"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="478.79077"
+       cx="522.14288"
+       id="path4136-89"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="485.93362"
+       cx="541.42859"
+       id="path4136-47"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="493.79077"
+       cx="556.42859"
+       id="path4136-58"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="507.36221"
+       cx="576.42859"
+       id="path4136-81"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="509.50507"
+       cx="590.71429"
+       id="path4136-42"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="510.93362"
+       cx="603.57141"
+       id="path4136-11"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="632.36218"
+       cx="309.28571"
+       id="path4136-59"
+       style="fill:#0095d7;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="620.36218"
+       cx="366.42856"
+       id="path4136-9"
+       style="fill:#0095d7;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="586.50507"
+       cx="333.71429"
+       id="path4136-75"
+       style="fill:#0095d7;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <path
+       id="path4329"
+       d="m 140,342.3622 500,0"
+       style="fill:none;fill-rule:evenodd;stroke:#f62929;stroke-width:4;stroke-linecap:round;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:32, 4;stroke-dashoffset:0;stroke-opacity:0.75117373" />
+    <text
+       id="text4331"
+       y="332.36221"
+       x="140"
+       style="font-style:normal;font-weight:normal;font-size:20px;line-height:125%;font-family:Sans;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+       xml:space="preserve"><tspan
+         y="332.36221"
+         x="140"
+         id="tspan4333">zrange=180,<tspan
+   id="tspan4335"
+   style="font-weight:bold">250</tspan></tspan></text>
+    <path
+       id="path4329-4"
+       d="m 140,560.3622 500,0"
+       style="fill:none;fill-rule:evenodd;stroke:#f62929;stroke-width:4;stroke-linecap:round;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:32, 4;stroke-dashoffset:0;stroke-opacity:0.75117373" />
+    <text
+       id="text4331-8"
+       y="550.36218"
+       x="140"
+       style="font-style:normal;font-weight:normal;font-size:20px;line-height:125%;font-family:Sans;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+       xml:space="preserve"><tspan
+         y="550.36218"
+         x="140"
+         id="tspan4333-5">zrange=<tspan
+   id="tspan4405"
+   style="font-weight:bold">180</tspan>,250</tspan></text>
+    <circle
+       r="4.1224623"
+       cy="318.07648"
+       cx="465"
+       id="path4136-5"
+       style="fill:#0095d7;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="298.07648"
+       cx="545"
+       id="path4136-5-1"
+       style="fill:#0095d7;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="498.07648"
+       cx="215"
+       id="path4136-5-8"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <circle
+       r="4.1224623"
+       cy="518.07648"
+       cx="195"
+       id="path4136-5-2"
+       style="fill:#00d710;fill-opacity:1;stroke:none;stroke-width:0.47379306;stroke-linecap:round;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" />
+    <text
+       id="text4459"
+       y="582.36218"
+       x="340"
+       style="font-style:normal;font-weight:normal;font-size:12.5px;line-height:125%;font-family:Sans;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+       xml:space="preserve"><tspan
+         y="582.36218"
+         x="340"
+         id="tspan4461">175</tspan></text>
+    <text
+       id="text4463"
+       y="626.36218"
+       x="312"
+       style="font-style:normal;font-weight:normal;font-size:12.5px;line-height:125%;font-family:Sans;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+       xml:space="preserve"><tspan
+         y="626.36218"
+         x="312"
+         id="tspan4465">160</tspan></text>
+    <text
+       id="text4467"
+       y="612.36218"
+       x="370"
+       style="font-style:normal;font-weight:normal;font-size:12.5px;line-height:125%;font-family:Sans;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+       xml:space="preserve"><tspan
+         y="612.36218"
+         x="370"
+         id="tspan4469">162</tspan></text>
+    <text
+       id="text4459-1"
+       y="310.36218"
+       x="469"
+       style="font-style:normal;font-weight:normal;font-size:12.5px;line-height:125%;font-family:Sans;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+       xml:space="preserve"><tspan
+         y="310.36218"
+         x="469"
+         id="tspan4461-5">257</tspan></text>
+    <text
+       id="text4459-9"
+       y="290.86053"
+       x="550.5116"
+       style="font-style:normal;font-weight:normal;font-size:12.5px;line-height:125%;font-family:Sans;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+       xml:space="preserve"><tspan
+         y="290.86053"
+         x="550.5116"
+         id="tspan4461-8">235</tspan></text>
+  </g>
+</svg>