r.in.pdal – a PDAL based replacement of r.in.lidar (#1200)

* r.in.pdal: point cloud binning with PDAL

Reads points using PDAL, uses several PDAL filters including reprojection,
filters using custom filter, writes into memory structures.

Reuses r.in.lidar functionality with moving code into PDAL Filter
and Writer classes. Preserves most of r.in.lidar functionality,
but misses several features most notably scanning and auto extent.
Furthermore, it now assumes that the whole output raster fits
into memory. The multi-pass code was removed and it should be
replaced by use of the improved Segment with all-in-memory mode.

In comparison to r.in.lidar, here are more binning options, including eigenvalue statistic option according to suggestion of "Mallet et al. 2011. Relevance assessment of full-waveform lidar data for urban area classification"; Summation is performed with Neumaiers improved Kahan–Babuska algorithm; Variance, stddev and coeff_var with Welfords algorithm.

Testing infrastructure is also updated to better handle minimal version required of PDAL library.

Co-authored-by: Vaclav Petras <wenzeslaus@gmail.com>
Co-authored-by: Māris Nartišs <maris.nartiss@lu.lv>

.github/workflows/build_ubuntu18.sh

@@ -0,0 +1,45 @@
+#!/usr/bin/env bash
+
+# The make step requires something like:
+# export LD_LIBRARY_PATH="$LD_LIBRARY_PATH:$PREFIX/lib"
+# further steps additionally require:
+# export PATH="$PATH:$PREFIX/bin"
+
+# fail on non-zero return code from a subprocess
+set -e
+
+# print commands
+set -x
+
+if [ -z "$1" ]; then
+    echo "Usage: $0 PREFIX"
+    exit 1
+fi
+
+# non-existent variables as an errors
+set -u
+
+export INSTALL_PREFIX=$1
+
+./configure \
+    --prefix="$INSTALL_PREFIX/" \
+    --enable-largefile \
+    --with-cxx \
+    --with-zstd \
+    --with-bzlib \
+    --with-blas \
+    --with-lapack \
+    --with-readline \
+    --with-openmp \
+    --with-pthread \
+    --with-tiff \
+    --with-freetype \
+    --with-freetype-includes="/usr/include/freetype2/" \
+    --with-proj-share=/usr/share/proj \
+    --with-geos \
+    --with-sqlite \
+    --with-fftw \
+    --with-netcdf
+
+make
+make install

.github/workflows/ci.yml

@@ -1,4 +1,4 @@
-name: CI
+name: Ubuntu18
 
 on:
   - push
@@ -13,7 +13,6 @@ jobs:
       matrix:
         os:
           - ubuntu-18.04
-          - ubuntu-20.04
       fail-fast: false
 
     steps:
@@ -34,7 +33,7 @@ jobs:
         run: |
           echo "LD_LIBRARY_PATH=$HOME/install/lib" >> $GITHUB_ENV
       - name: Build
-        run: .github/workflows/build.sh $HOME/install
+        run: .github/workflows/build_ubuntu18.sh $HOME/install
 
   test:
     name: ${{ matrix.os }} tests
@@ -44,7 +43,6 @@ jobs:
       matrix:
         os:
           - ubuntu-18.04
-          - ubuntu-20.04
       fail-fast: false
 
     steps:
@@ -65,7 +63,7 @@ jobs:
         run: |
           echo "LD_LIBRARY_PATH=$HOME/install/lib" >> $GITHUB_ENV
       - name: Build
-        run: .github/workflows/build.sh $HOME/install
+        run: .github/workflows/build_ubuntu18.sh $HOME/install
       - name: Add the bin directory to PATH
         run: |
           echo "$HOME/install/bin" >> $GITHUB_PATH

.github/workflows/ci_ubuntu20.yml

@@ -0,0 +1,81 @@
+name: Ubuntu20
+
+on:
+  - push
+  - pull_request
+
+jobs:
+  build:
+    name: ${{ matrix.os }} build
+
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os:
+          - ubuntu-20.04
+      fail-fast: false
+
+    steps:
+      - uses: actions/checkout@v2
+      - name: Get dependencies
+        run: |
+          sudo apt-get update -y
+          sudo apt-get install -y wget git gawk findutils
+          xargs -a <(awk '! /^ *(#|$)/' ".github/workflows/apt.txt") -r -- \
+              sudo apt-get install -y --no-install-recommends --no-install-suggests
+      - name: Create installation directory
+        run: |
+          mkdir $HOME/install
+      - name: Ensure one core for compilation
+        run: |
+          echo "MAKEFLAGS=-j1" >> $GITHUB_ENV
+      - name: Set LD_LIBRARY_PATH for compilation
+        run: |
+          echo "LD_LIBRARY_PATH=$HOME/install/lib" >> $GITHUB_ENV
+      - name: Build
+        run: .github/workflows/build_ubuntu20.sh $HOME/install
+
+  test:
+    name: ${{ matrix.os }} tests
+
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os:
+          - ubuntu-20.04
+      fail-fast: false
+
+    steps:
+      - uses: actions/checkout@v2
+      - name: Get dependencies
+        run: |
+          sudo apt-get update -y
+          sudo apt-get install -y wget git gawk findutils
+          xargs -a <(awk '! /^ *(#|$)/' ".github/workflows/apt.txt") -r -- \
+              sudo apt-get install -y --no-install-recommends --no-install-suggests
+      - name: Create installation directory
+        run: |
+          mkdir $HOME/install
+      - name: Set number of cores for compilation
+        run: |
+          echo "MAKEFLAGS=-j$(nproc)" >> $GITHUB_ENV
+      - name: Set LD_LIBRARY_PATH for compilation
+        run: |
+          echo "LD_LIBRARY_PATH=$HOME/install/lib" >> $GITHUB_ENV
+      - name: Build
+        run: .github/workflows/build_ubuntu20.sh $HOME/install
+      - name: Add the bin directory to PATH
+        run: |
+          echo "$HOME/install/bin" >> $GITHUB_PATH
+      - name: Test executing of the grass command
+        run: .github/workflows/test_simple.sh
+      - name: Run tests
+        run: .github/workflows/test_thorough.sh
+
+      - name: Make HTML test report available
+        if: ${{ always() }}
+        uses: actions/upload-artifact@v2
+        with:
+          name: testreport-${{ matrix.os }}
+          path: testreport
+          retention-days: 3

.github/workflows/codeql-analysis.yml

@@ -50,7 +50,7 @@ jobs:
         env:
           CFLAGS: "-std=gnu11"
           CXXFLAGS: "-std=c++11"
-        run: .github/workflows/build.sh $HOME/install
+        run: .github/workflows/build_ubuntu20.sh $HOME/install
 
       - name: Perform CodeQL Analysis
         uses: github/codeql-action/analyze@v1

.github/workflows/gcc.yml

@@ -45,4 +45,4 @@ jobs:
           CFLAGS: "-std=${{ matrix.c }}"
           # TODO: -pedantic-errors here won't compile
           CXXFLAGS: "-std=${{ matrix.cpp }}"
-        run: .github/workflows/build.sh $HOME/install
+        run: .github/workflows/build_ubuntu20.sh $HOME/install

configure.in

@@ -1011,9 +1011,13 @@ else
   LIBS="$LIBS $PDAL_LIBS"
   CFLAGS="$CFLAGS $PDAL_CFLAGS"
   CPPFLAGS="$CPPFLAGS $PDAL_CPPFLAGS $PDAL_INC"
-  AC_TRY_LINK([#include <pdal/PointTable.hpp>],[pdal::PointTable table;],,[
-  AC_TRY_LINK([#include <pdal/PointTable.hpp>],[pdal::PointTable table;],PDAL_LIBS="$PDAL_LIBS",[
-  AC_MSG_ERROR([*** Unable to locate PDAL library.])
+  AC_TRY_LINK([#include <pdal/PointTable.hpp>
+  #include <pdal/Streamable.hpp>
+  class St:public pdal::Streamable {};],[pdal::PointTable table;],,[
+  AC_TRY_LINK([#include <pdal/PointTable.hpp>
+  #include <pdal/Streamable.hpp>
+  class St:public pdal::Streamable {};],[pdal::PointTable table;],PDAL_LIBS="$PDAL_LIBS",[
+  AC_MSG_ERROR([*** Unable to locate suitable (>=1.7.1) PDAL library.])
   ])
   ])
   LIBS=${ac_save_libs}

gui/wxpython/xml/toolboxes.xml

@@ -302,6 +302,9 @@
       <module-item name="r.in.lidar">
         <label>LAS LiDAR points import</label>
       </module-item>
+      <module-item name="r.in.pdal">
+        <label>Point cloud (LAS LiDAR) import</label>
+      </module-item>
       <separator/>
       <module-item name="r.unpack">
         <label>Unpack raster map</label>

raster/Makefile

@@ -35,6 +35,7 @@ SUBDIRS = \
 	r.in.gridatb \
 	r.in.lidar \
 	r.in.mat \
+	r.in.pdal \
 	r.in.png \
 	r.in.poly \
 	r.in.xyz \

raster/r.in.pdal/Makefile

@@ -0,0 +1,23 @@
+MODULE_TOPDIR = ../..
+
+PGM=r.in.pdal
+
+LIBES = $(RASTERLIB) $(SEGMENTLIB) $(GPROJLIB) $(VECTORLIB) $(DBMILIB) $(GISLIB) $(MATHLIB) $(PDALLIBS) $(GMATHLIB)
+DEPENDENCIES = $(GPROJDEP) $(VECTORDEP) $(DBMIDEP) $(RASTERDEP) $(GISDEP) $(SEGMENTDEP) $(GMATHDEP)
+
+EXTRA_INC = $(VECT_INC) $(PROJINC) $(PDALINC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+$(OBJDIR)/main.o : EXTRA_CFLAGS += $(PDALCPPFLAGS)
+
+$(OBJDIR)/grasslidarfilter.o : EXTRA_CFLAGS += $(PDALCPPFLAGS)
+
+LINK = $(CXX)
+
+ifneq ($(strip $(CXX)),)
+ifneq ($(USE_PDAL),)
+default: cmd
+endif
+endif

raster/r.in.pdal/bin_update.c

@@ -0,0 +1,381 @@
+/*
+ * r.in.pdal Functions performing value updates on each incoming point
+ *            during point binning process
+ *   Copyright 2006 by M. Hamish Bowman, and The GRASS Development Team
+ *   Author: M. Hamish Bowman, University of Otago, Dunedin, New Zealand
+ *     Maris Nartiss code refactoring for r.in.pdal
+ *
+ *   This program is free software licensed under the GPL (>=v2).
+ *   Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+#include <stddef.h>
+#include <stdlib.h>
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+
+#include "point_binning.h"
+#include "bin_update.h"
+
+
+static int new_node(struct BinIndex *bin_index, size_t size)
+{
+    int n = bin_index->num_nodes++;
+
+    if (bin_index->num_nodes >= bin_index->max_nodes) {
+        bin_index->max_nodes += SIZE_INCREMENT;
+        bin_index->nodes = G_realloc(bin_index->nodes,
+                                     (size_t)bin_index->max_nodes * size);
+    }
+
+    return n;
+}
+
+
+void update_val(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);
+
+    Rast_set_d_value(ptr, value, map_type);
+    return;
+}
+
+
+void 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;
+}
+
+
+void 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;
+}
+
+
+void 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;
+}
+
+/* Implements improved Kahan–Babuska algorithm by Neumaier, A. 1974 */
+void update_sum(void *sum_array, void *c_array, int cols, int row, int col,
+                RASTER_MAP_TYPE map_type, double value)
+{
+    DCELL old_sum, old_c, tmp;
+    void *s_ptr = get_cell_ptr(sum_array, cols, row, col, map_type);
+    void *c_ptr = get_cell_ptr(c_array, cols, row, col, map_type);
+
+    old_sum = Rast_get_d_value(s_ptr, map_type);
+    old_c = Rast_get_d_value(c_ptr, map_type);
+    tmp = old_sum + value;
+
+    if (abs(old_sum) >= abs(value))
+        Rast_set_d_value(c_ptr, old_c + (old_sum - tmp) + value, map_type);
+    else
+        Rast_set_d_value(c_ptr, old_c + (value - tmp) + old_sum, map_type);
+
+    Rast_set_d_value(s_ptr, tmp, map_type);
+
+    return;
+}
+
+
+/* Implements Welford algorithm */
+void update_m2(void *n_array, void *mean_array, void *m2_array, int cols,
+               int row, int col, RASTER_MAP_TYPE map_type, double value)
+{
+    int n;
+    double m2, mean, d1, d2;
+    void *n_ptr = get_cell_ptr(n_array, cols, row, col, CELL_TYPE);
+    void *mean_ptr = get_cell_ptr(mean_array, cols, row, col, map_type);
+    void *m2_ptr = get_cell_ptr(m2_array, cols, row, col, map_type);
+
+    n = Rast_get_c_value(n_ptr, CELL_TYPE);
+    mean = Rast_get_d_value(mean_ptr, map_type);
+    m2 = Rast_get_d_value(m2_ptr, map_type);
+
+    n++;
+    Rast_set_c_value(n_ptr, n, CELL_TYPE);
+
+    if (n == 1) {
+        Rast_set_d_value(mean_ptr, value, map_type);
+        return;
+    }
+
+    d1 = value - mean;
+    mean = mean + (d1 / n);
+    d2 = value - mean;
+    m2 = m2 + (d1 * d2);
+
+    Rast_set_d_value(mean_ptr, mean, map_type);
+    Rast_set_d_value(m2_ptr, m2, map_type);
+
+    return;
+}
+
+
+void update_moving_mean(void *array, int cols, int row, int col,
+                        RASTER_MAP_TYPE rtype, double value, int n)
+{
+    /* for xy we do this check twice */
+    if (n != 0) {
+        double m_v;
+
+        row_array_get_value_row_col(array, row, col, cols, rtype, &m_v);
+        value = m_v + (value - m_v) / n;
+    }
+    /* else we just write the initial value */
+    update_val(array, cols, row, col, rtype, value);
+    return;
+}
+
+
+/* add node to sorted, single linked list
+ * returns id if head has to be saved to index array, otherwise -1 */
+int add_z_node(struct BinIndex *bin_index, int head, double z)
+{
+    int node_id, last_id, newnode_id, head_id;
+
+    head_id = head;
+    node_id = head_id;
+    last_id = head_id;
+
+    while (node_id != -1 &&
+           ((struct z_node *)bin_index->nodes)[node_id].z < z) {
+        last_id = node_id;
+        node_id = ((struct z_node *)bin_index->nodes)[node_id].next;
+    }
+
+    /* end of list, simply append */
+    if (node_id == -1) {
+        newnode_id = new_node(bin_index, sizeof(struct z_node));
+        ((struct z_node *)bin_index->nodes)[newnode_id].next = -1;
+        ((struct z_node *)bin_index->nodes)[newnode_id].z = z;
+        ((struct z_node *)bin_index->nodes)[last_id].next = newnode_id;
+        return -1;
+    }
+    else if (node_id == head_id) {      /* pole position, insert as head */
+        newnode_id = new_node(bin_index, sizeof(struct z_node));
+        ((struct z_node *)bin_index->nodes)[newnode_id].next = head_id;
+        head_id = newnode_id;
+        ((struct z_node *)bin_index->nodes)[newnode_id].z = z;
+        return (head_id);
+    }
+    else {                      /* somewhere in the middle, insert */
+        newnode_id = new_node(bin_index, sizeof(struct z_node));
+        ((struct z_node *)bin_index->nodes)[newnode_id].z = z;
+        ((struct z_node *)bin_index->nodes)[newnode_id].next = node_id;
+        ((struct z_node *)bin_index->nodes)[last_id].next = newnode_id;
+        return -1;
+    }
+}
+
+
+void add_cnt_node(struct BinIndex *bin_index, int head, int value)
+{
+    int node_id, newnode_id, head_id, next;
+
+    head_id = head;
+    node_id = head_id;
+    next = node_id;
+
+    while (next != -1) {
+        node_id = next;
+        if (((struct cnt_node *)bin_index->nodes)[node_id].value == value) {
+            ((struct cnt_node *)bin_index->nodes)[node_id].count++;
+            return;
+        }
+        next = ((struct cnt_node *)bin_index->nodes)[node_id].next;
+    }
+
+    newnode_id = new_node(bin_index, sizeof(struct cnt_node));
+    ((struct cnt_node *)bin_index->nodes)[newnode_id].next = -1;
+    ((struct cnt_node *)bin_index->nodes)[newnode_id].value = value;
+    ((struct cnt_node *)bin_index->nodes)[newnode_id].count = 1;
+    ((struct cnt_node *)bin_index->nodes)[node_id].next = newnode_id;
+    return;
+}
+
+
+/* Unlike the other functions, this one is not using map_type (RASTER_MAP_TYPE)
+ * because the values (z) are always doubles and the index is integer. */
+void update_bin_z_index(struct BinIndex *bin_index, void *index_array,
+                        int cols, int row, int col, double value)
+{
+    int head_id;
+    void *ptr = index_array;
+
+    ptr =
+        G_incr_void_ptr(ptr,
+                        (((size_t)row * cols) +
+                         col) * Rast_cell_size(CELL_TYPE));
+
+    /* first node */
+    if (Rast_is_null_value(ptr, CELL_TYPE)) {
+        head_id = new_node(bin_index, sizeof(struct z_node));
+        ((struct z_node *)bin_index->nodes)[head_id].next = -1;
+        ((struct z_node *)bin_index->nodes)[head_id].z = value;
+        /* store index to head */
+        Rast_set_c_value(ptr, head_id, CELL_TYPE);
+    }
+    /* head is already there */
+    else {
+
+        /* get index to head */
+        head_id = Rast_get_c_value(ptr, CELL_TYPE);
+        head_id = add_z_node(bin_index, head_id, value);
+        /* if id valid, store index to head */
+        if (head_id != -1)
+            Rast_set_c_value(ptr, head_id, CELL_TYPE);
+    }
+
+    return;
+}
+
+
+void update_bin_cnt_index(struct BinIndex *bin_index, void *index_array,
+                          int cols, int row, int col, int value)
+{
+    int head_id;
+    void *ptr = index_array;
+
+    ptr =
+        G_incr_void_ptr(ptr,
+                        (((size_t)row * cols) +
+                         col) * Rast_cell_size(CELL_TYPE));
+
+    /* first node */
+    if (Rast_is_null_value(ptr, CELL_TYPE)) {
+        head_id = new_node(bin_index, sizeof(struct cnt_node));
+        ((struct cnt_node *)bin_index->nodes)[head_id].next = -1;
+        ((struct cnt_node *)bin_index->nodes)[head_id].value = value;
+        ((struct cnt_node *)bin_index->nodes)[head_id].count = 1;
+        /* store index to head */
+        Rast_set_c_value(ptr, head_id, CELL_TYPE);
+    }
+    /* head is already there */
+    else {
+        head_id = Rast_get_c_value(ptr, CELL_TYPE);
+        add_cnt_node(bin_index, head_id, value);
+    }
+
+    return;
+}
+
+
+/* Co-moment value update */
+void update_com_node(struct com_node *cn, int item, double x, double y)
+{
+    double dx;
+
+    dx = x - cn->meanx[item];
+    cn->meanx[item] = cn->meanx[item] + (dx / cn->n);
+    cn->meany[item] = cn->meany[item] + (y - cn->meany[item]) / cn->n;
+    cn->comoment[item] = cn->comoment[item] + dx * (y - cn->meany[item]);
+}
+
+
+void update_bin_com_index(struct BinIndex *bin_index, void *index_array,
+                          int cols, int row, int col,
+                          double x, double y, double z)
+{
+    int node_id;
+    void *ptr = index_array;
+
+    ptr =
+        G_incr_void_ptr(ptr,
+                        (((size_t)row * cols) +
+                         col) * Rast_cell_size(CELL_TYPE));
+
+    if (Rast_is_null_value(ptr, CELL_TYPE)) {
+        node_id = new_node(bin_index, sizeof(struct com_node));
+
+        ((struct com_node *)bin_index->nodes)[node_id].n = 0;
+        ((struct com_node *)bin_index->nodes)[node_id].meanx =
+            (double *)G_calloc(6, sizeof(double));
+        ((struct com_node *)bin_index->nodes)[node_id].meany =
+            (double *)G_calloc(6, sizeof(double));
+        ((struct com_node *)bin_index->nodes)[node_id].comoment =
+            (double *)G_calloc(6, sizeof(double));
+
+        /* store index to node */
+        Rast_set_c_value(ptr, node_id, CELL_TYPE);
+    }
+    /* node is already there */
+    else {
+        node_id = Rast_get_c_value(ptr, CELL_TYPE);
+    }
+
+    /* update values */
+    ((struct com_node *)bin_index->nodes)[node_id].n++;
+    update_com_node(&(((struct com_node *)bin_index->nodes)[node_id]),
+                    0, x, x);
+    update_com_node(&(((struct com_node *)bin_index->nodes)[node_id]),
+                    1, x, y);
+    update_com_node(&(((struct com_node *)bin_index->nodes)[node_id]),
+                    2, x, z);
+    update_com_node(&(((struct com_node *)bin_index->nodes)[node_id]),
+                    3, y, y);
+    update_com_node(&(((struct com_node *)bin_index->nodes)[node_id]),
+                    4, x, z);
+    update_com_node(&(((struct com_node *)bin_index->nodes)[node_id]),
+                    5, z, z);
+
+    return;
+}
+
+
+/* 0 on NULL, 1 on success */
+int row_array_get_value_row_col(void *array, int arr_row, int arr_col,
+                                int cols, RASTER_MAP_TYPE rtype,
+                                double *value)
+{
+    void *ptr = array;
+
+    ptr =
+        G_incr_void_ptr(ptr,
+                        (((size_t)arr_row * cols) +
+                         arr_col) * Rast_cell_size(rtype));
+    if (Rast_is_null_value(ptr, rtype))
+        return 0;
+    if (rtype == DCELL_TYPE)
+        *value = (double)*(DCELL *) ptr;
+    else if (rtype == FCELL_TYPE)
+        *value = (double)*(FCELL *) ptr;
+    else
+        *value = (double)*(CELL *) ptr;
+    return 1;
+}

raster/r.in.pdal/bin_update.h

@@ -0,0 +1,54 @@
+
+ /****************************************************************************
+ *
+ * MODULE:    r.in.pdal
+ *
+ * AUTHOR(S): Vaclav Petras
+ *            Based on r.in.xyz and r.in.lidar by Markus Metz,
+ *            Hamish Bowman, Volker Wichmann
+ *            Maris Nartiss code refactoring
+ *
+ * PURPOSE:   Functions performing value updates on each incoming point
+ *            during point binning process
+ *
+ * COPYRIGHT: (C) 2011-2021 by Vaclav Petras and 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 __BIN_UPDATE_H__
+#define __BIN_UPDATE_H__
+
+#include <grass/raster.h>
+
+#define SIZE_INCREMENT 16
+
+/* forward declarations for point_binning.h */
+struct BinIndex;
+struct com_node;
+
+void update_n(void *, int, int, int);
+void update_min(void *, int, int, int, RASTER_MAP_TYPE, double);
+void update_max(void *, int, int, int, RASTER_MAP_TYPE, double);
+void update_sum(void *, void *, int, int, int, RASTER_MAP_TYPE, double);
+void update_m2(void *, void *, void *, int, int, int, RASTER_MAP_TYPE,
+               double);
+void update_moving_mean(void *, int, int, int, RASTER_MAP_TYPE, double, int);
+
+int add_z_node(struct BinIndex *, int, double);
+void add_cnt_node(struct BinIndex *, int, int);
+
+void update_bin_z_index(struct BinIndex *, void *, int, int, int, double);
+void update_bin_cnt_index(struct BinIndex *, void *, int, int, int, int);
+void update_com_node(struct com_node *, int, double, double);
+void update_bin_com_index(struct BinIndex *, void *,
+                          int, int, int, double, double, double);
+
+int row_array_get_value_row_col(void *, int, int,
+                                int, RASTER_MAP_TYPE, double *);
+
+
+#endif /* __BIN_UPDATE_H__ */

raster/r.in.pdal/bin_write.c

@@ -0,0 +1,467 @@
+/*
+ * r.in.pdal point binning output value calculation
+ *
+ * Copyright 2011-2015, 2020 by Markus Metz, and The GRASS Development Team
+ * Authors:
+ *  Markus Metz (r.in.lidar)
+ *  Vaclav Petras (move code to separate functions)
+ *  Maris Nartiss (refactoring for r.in.pdal)
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+#include <stddef.h>
+#include <math.h>
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/gmath.h>
+
+#include "point_binning.h"
+#include "bin_write.h"
+
+
+/* Get error corrected sum */
+double get_sum(void *sum_array, void *c_array,
+               int row, int cols, int col, RASTER_MAP_TYPE rtype)
+{
+    size_t offset = ((size_t)row * cols + col) * Rast_cell_size(rtype);
+    double sum = Rast_get_d_value(((char *)sum_array) + offset, rtype);
+    double c = Rast_get_d_value(((char *)c_array) + offset, rtype);
+
+    return sum + c;
+}
+
+
+void write_sum(void *raster_row, void *sum_array, void *c_array,
+               int row, int cols, RASTER_MAP_TYPE rtype)
+{
+    int col;
+    void *ptr = raster_row;
+
+    for (col = 0; col < cols; col++) {
+        double sum = get_sum(sum_array, c_array, row, cols, col, rtype);
+
+        Rast_set_d_value(ptr, sum, rtype);
+        ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
+    }
+}
+
+
+void write_variance(void *raster_row, void *n_array, void *mean_array,
+                    void *m2_array, int row, int cols,
+                    RASTER_MAP_TYPE rtype, int method)
+{
+    double variance;
+    int col;
+    void *ptr = raster_row;
+
+    for (col = 0; col < cols; col++) {
+        size_t offset = ((size_t)row * cols + col) * Rast_cell_size(rtype);
+        size_t n_offset =
+            ((size_t)row * cols + col) * Rast_cell_size(CELL_TYPE);
+        int n = Rast_get_c_value(((char *)n_array) + n_offset, CELL_TYPE);
+        double mean = Rast_get_d_value(((char *)mean_array) + offset, rtype);
+        double m2 = Rast_get_d_value(((char *)m2_array) + offset, rtype);
+
+        if (n == 0)
+            Rast_set_null_value(ptr, 1, rtype);
+        else if (n == 1)
+            Rast_set_d_value(ptr, 0.0, rtype);
+        else {
+            variance = m2 / n;
+
+            if (variance < GRASS_EPSILON)
+                variance = 0.0;
+
+            if (method == METHOD_STDDEV)
+                variance = sqrt(variance);
+
+            else if (method == METHOD_COEFF_VAR)
+                variance = 100 * sqrt(variance) / mean;
+
+            Rast_set_d_value(ptr, variance, rtype);
+        }
+        ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
+    }
+}
+
+
+void write_median(struct BinIndex *bin_index, void *raster_row,
+                  void *index_array, int row, int cols, RASTER_MAP_TYPE rtype)
+{
+    int n;
+    int j;
+    double z;
+    int col;
+    int node_id, head_id;
+    void *ptr = raster_row;
+
+    for (col = 0; col < cols; col++) {
+        size_t n_offset =
+            ((size_t)row * cols + col) * Rast_cell_size(CELL_TYPE);
+        if (Rast_is_null_value(((char *)index_array) + n_offset, CELL_TYPE))    /* no points in cell */
+            Rast_set_null_value(ptr, 1, rtype);
+        else {                  /* one or more points in cell */
+
+            head_id =
+                Rast_get_c_value(((char *)index_array) + n_offset, CELL_TYPE);
+            node_id = head_id;
+
+            n = 0;
+
+            while (node_id != -1) {     /* count number of points in cell */
+                n++;
+                node_id = ((struct z_node *)bin_index->nodes)[node_id].next;
+            }
+
+            if (n == 1)         /* only one point, use that */
+                Rast_set_d_value(ptr,
+                                 ((struct z_node *)bin_index->
+                                  nodes)[head_id].z, rtype);
+            else if (n % 2 != 0) {      /* odd number of points: median_i = (n + 1) / 2 */
+                n = (n + 1) / 2;
+                node_id = head_id;
+                for (j = 1; j < n; j++) /* get "median element" */
+                    node_id =
+                        ((struct z_node *)bin_index->nodes)[node_id].next;
+
+                Rast_set_d_value(ptr,
+                                 ((struct z_node *)bin_index->
+                                  nodes)[node_id].z, rtype);
+            }
+            else {              /* even number of points: median = (val_below + val_above) / 2 */
+
+                z = (n + 1) / 2.0;
+                n = floor(z);
+                node_id = head_id;
+                for (j = 1; j < n; j++) /* get element "below" */
+                    node_id =
+                        ((struct z_node *)bin_index->nodes)[node_id].next;
+
+                z = (((struct z_node *)bin_index->nodes)[node_id].z +
+                     ((struct z_node *)
+                      bin_index->nodes)[((struct z_node *)bin_index->
+                                         nodes)[node_id].next].z) / 2;
+                Rast_set_d_value(ptr, z, rtype);
+            }
+        }
+        ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
+    }
+}
+
+
+void write_mode(struct BinIndex *bin_index, void *raster_row,
+                void *index_array, int row, int cols)
+{
+    int col;
+    int node_id;
+    void *ptr = raster_row;
+
+    for (col = 0; col < cols; col++) {
+        size_t n_offset =
+            ((size_t)row * cols + col) * Rast_cell_size(CELL_TYPE);
+        if (Rast_is_null_value(((char *)index_array) + n_offset, CELL_TYPE))    /* no points in cell */
+            Rast_set_null_value(ptr, 1, CELL_TYPE);
+        else {
+            int mode_node = -1;
+
+            node_id =
+                Rast_get_c_value(((char *)index_array) + n_offset, CELL_TYPE);
+
+            while (node_id != -1) {
+                if (mode_node == -1)
+                    mode_node = node_id;
+                else if (((struct cnt_node *)bin_index->
+                          nodes)[node_id].count >
+                         ((struct cnt_node *)bin_index->
+                          nodes)[mode_node].count)
+                    mode_node = node_id;
+                node_id = ((struct cnt_node *)bin_index->nodes)[node_id].next;
+            }
+            Rast_set_c_value(ptr,
+                             ((struct cnt_node *)bin_index->
+                              nodes)[mode_node].value, CELL_TYPE);
+        }
+        ptr = G_incr_void_ptr(ptr, Rast_cell_size(CELL_TYPE));
+    }
+}
+
+
+void write_percentile(struct BinIndex *bin_index, void *raster_row,
+                      void *index_array, int row, int cols,
+                      RASTER_MAP_TYPE rtype, int pth)
+{
+    int n;
+    int j;
+    double z;
+    int col;
+    int node_id, head_id;
+    int r_low, r_up;
+    void *ptr = raster_row;
+
+    for (col = 0; col < cols; col++) {
+        size_t n_offset =
+            ((size_t)row * cols + col) * Rast_cell_size(CELL_TYPE);
+        if (Rast_is_null_value(((char *)index_array) + n_offset, CELL_TYPE))    /* no points in cell */
+            Rast_set_null_value(ptr, 1, rtype);
+        else {
+            head_id =
+                Rast_get_c_value(((char *)index_array) + n_offset, CELL_TYPE);
+            node_id = head_id;
+            n = 0;
+
+            while (node_id != -1) {     /* count number of points in cell */
+                n++;
+                node_id = ((struct z_node *)bin_index->nodes)[node_id].next;
+            }
+
+            z = (pth * (n + 1)) / 100.0;
+            r_low = floor(z);   /* lower rank */
+            if (r_low < 1)
+                r_low = 1;
+            else if (r_low > n)
+                r_low = n;
+
+            r_up = ceil(z);     /* upper rank */
+            if (r_up > n)
+                r_up = n;
+
+            node_id = head_id;
+            for (j = 1; j < r_low; j++) /* search lower value */
+                node_id = ((struct z_node *)bin_index->nodes)[node_id].next;
+
+            z = ((struct z_node *)bin_index->nodes)[node_id].z; /* save lower value */
+            node_id = head_id;
+            for (j = 1; j < r_up; j++)  /* search upper value */
+                node_id = ((struct z_node *)bin_index->nodes)[node_id].next;
+
+            z = (z + ((struct z_node *)bin_index->nodes)[node_id].z) / 2;
+            Rast_set_d_value(ptr, z, rtype);
+        }
+        ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
+    }
+}
+
+
+void write_skewness(struct BinIndex *bin_index, void *raster_row,
+                    void *index_array, int row, int cols,
+                    RASTER_MAP_TYPE rtype)
+{
+    int n;
+    double z;
+    int col;
+    int node_id, head_id;
+    double variance, mean, skew, sumdev;
+    double sum;
+    double sumsq;
+    void *ptr = raster_row;
+
+    for (col = 0; col < cols; col++) {
+        size_t n_offset =
+            ((size_t)row * cols + col) * Rast_cell_size(CELL_TYPE);
+        if (Rast_is_null_value(((char *)index_array) + n_offset, CELL_TYPE))    /* no points in cell */
+            Rast_set_null_value(ptr, 1, rtype);
+        else {
+            head_id =
+                Rast_get_c_value(((char *)index_array) + n_offset, CELL_TYPE);
+            node_id = head_id;
+
+            n = 0;              /* count */
+            sum = 0.0;          /* sum */
+            sumsq = 0.0;        /* sum of squares */
+            sumdev = 0.0;       /* sum of (xi - mean)^3 */
+            skew = 0.0;         /* skewness */
+
+            while (node_id != -1) {
+                z = ((struct z_node *)bin_index->nodes)[node_id].z;
+                n++;
+                sum += z;
+                sumsq += (z * z);
+                node_id = ((struct z_node *)bin_index->nodes)[node_id].next;
+            }
+
+            if (n > 1) {        /* if n == 1, skew is "0.0" */
+                mean = sum / n;
+                node_id = head_id;
+                while (node_id != -1) {
+                    z = ((struct z_node *)bin_index->nodes)[node_id].z;
+                    sumdev += pow((z - mean), 3);
+                    node_id =
+                        ((struct z_node *)bin_index->nodes)[node_id].next;
+                }
+
+                variance = (sumsq - sum * sum / n) / n;
+                if (variance < GRASS_EPSILON)
+                    skew = 0.0;
+                else
+                    skew = sumdev / ((n - 1) * pow(sqrt(variance), 3));
+            }
+            Rast_set_d_value(ptr, skew, rtype);
+        }
+        ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
+    }
+}
+
+
+void write_trimmean(struct BinIndex *bin_index, void *raster_row,
+                    void *index_array, int row, int cols,
+                    RASTER_MAP_TYPE rtype, double trim)
+{
+    int n;
+    int j, k;
+    int col;
+    int node_id, head_id;
+    double mean;
+    double sum;
+    void *ptr = raster_row;
+
+    for (col = 0; col < cols; col++) {
+        size_t n_offset =
+            ((size_t)row * cols + col) * Rast_cell_size(CELL_TYPE);
+        if (Rast_is_null_value(((char *)index_array) + n_offset, CELL_TYPE))    /* no points in cell */
+            Rast_set_null_value(ptr, 1, rtype);
+        else {
+            head_id =
+                Rast_get_c_value(((char *)index_array) + n_offset, CELL_TYPE);
+
+            node_id = head_id;
+            n = 0;
+            while (node_id != -1) {     /* count number of points in cell */
+                n++;
+                node_id = ((struct z_node *)bin_index->nodes)[node_id].next;
+            }
+
+            if (1 == n)
+                mean = ((struct z_node *)bin_index->nodes)[head_id].z;
+            else {
+                k = floor(trim * n + 0.5);      /* number of ranks to discard on each tail */
+
+                if (k > 0 && (n - 2 * k) > 0) { /* enough elements to discard */
+                    node_id = head_id;
+                    for (j = 0; j < k; j++)     /* move to first rank to consider */
+                        node_id =
+                            ((struct z_node *)bin_index->nodes)[node_id].next;
+
+                    j = k + 1;
+                    k = n - k;
+                    n = 0;
+                    sum = 0.0;
+
+                    while (j <= k) {    /* get values in interval */
+                        n++;
+                        sum += ((struct z_node *)bin_index->nodes)[node_id].z;
+                        node_id =
+                            ((struct z_node *)bin_index->nodes)[node_id].next;
+                        j++;
+                    }
+                }
+                else {
+                    node_id = head_id;
+                    n = 0;
+                    sum = 0.0;
+                    while (node_id != -1) {
+                        n++;
+                        sum += ((struct z_node *)bin_index->nodes)[node_id].z;
+                        node_id =
+                            ((struct z_node *)bin_index->nodes)[node_id].next;
+                    }
+                }
+                mean = sum / n;
+            }
+            Rast_set_d_value(ptr, mean, rtype);
+        }
+        ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
+    }
+}
+
+
+void write_sidn(struct BinIndex *bin_index, void *raster_row,
+                void *index_array, int row, int cols, int min)
+{
+    int col;
+    int node_id;
+    void *ptr = raster_row;
+    int count;
+
+    for (col = 0; col < cols; col++) {
+        size_t n_offset =
+            ((size_t)row * cols + col) * Rast_cell_size(CELL_TYPE);
+        if (Rast_is_null_value(((char *)index_array) + n_offset, CELL_TYPE))    /* no points in cell */
+            Rast_set_c_value(ptr, 0, CELL_TYPE);
+        else {
+
+            node_id =
+                Rast_get_c_value(((char *)index_array) + n_offset, CELL_TYPE);
+
+            count = ((struct cnt_node *)bin_index->nodes)[node_id].count;
+            node_id = ((struct cnt_node *)bin_index->nodes)[node_id].next;
+            while (node_id != -1) {
+                if (min &&
+                    ((struct cnt_node *)bin_index->nodes)[node_id].count <
+                    count)
+                    count =
+                        ((struct cnt_node *)bin_index->nodes)[node_id].count;
+                else if (!min &&
+                         ((struct cnt_node *)bin_index->
+                          nodes)[node_id].count > count)
+                    count =
+                        ((struct cnt_node *)bin_index->nodes)[node_id].count;
+                node_id = ((struct cnt_node *)bin_index->nodes)[node_id].next;
+            }
+            Rast_set_c_value(ptr, count, CELL_TYPE);
+        }
+        ptr = G_incr_void_ptr(ptr, Rast_cell_size(CELL_TYPE));
+    }
+}
+
+
+void write_ev(struct BinIndex *bin_index, void *raster_row, void *index_array,
+              int row, int cols, RASTER_MAP_TYPE rtype, int method)
+{
+    void *ptr = raster_row;
+    double **cov_matrix = G_alloc_matrix(3, 3);
+    double *ev = (double *)G_malloc(3 * sizeof(double));
+
+    for (int col = 0; col < cols; col++) {
+        size_t n_offset =
+            ((size_t)row * cols + col) * Rast_cell_size(CELL_TYPE);
+        if (Rast_is_null_value(((char *)index_array) + n_offset, CELL_TYPE))    /* no points in cell */
+            Rast_set_null_value(ptr, 1, rtype);
+        else {
+            int node_id;
+            struct com_node cn;
+
+            node_id =
+                Rast_get_c_value(((char *)index_array) + n_offset, CELL_TYPE);
+            cn = ((struct com_node *)bin_index->nodes)[node_id];
+
+            cov_matrix[0][0] = cn.comoment[0] / cn.n;
+            cov_matrix[0][1] = cov_matrix[1][0] = cn.comoment[1] / cn.n;
+            cov_matrix[0][2] = cov_matrix[2][0] = cn.comoment[2] / cn.n;
+            cov_matrix[1][1] = cn.comoment[3] / cn.n;
+            cov_matrix[1][2] = cov_matrix[2][1] = cn.comoment[4] / cn.n;
+            cov_matrix[2][2] = cn.comoment[5] / cn.n;
+
+            G_math_eigval(cov_matrix, ev, 3);
+
+            switch (method) {
+            case METHOD_EV1:
+                Rast_set_d_value(ptr, ev[0], rtype);
+                break;
+            case METHOD_EV2:
+                Rast_set_d_value(ptr, ev[1], rtype);
+                break;
+            case METHOD_EV3:
+                Rast_set_d_value(ptr, ev[1], rtype);
+                break;
+            }
+        }
+        ptr = G_incr_void_ptr(ptr, Rast_cell_size(CELL_TYPE));
+    }
+
+    G_free_matrix(cov_matrix);
+    G_free(ev);
+}

raster/r.in.pdal/bin_write.h

@@ -0,0 +1,41 @@
+/*
+ * r.in.pdal point binning output value calculation
+ *
+ * Copyright 2011-2015, 2020 by Markus Metz, and The GRASS Development Team
+ * Authors:
+ *  Markus Metz (r.in.lidar)
+ *  Vaclav Petras (move code to separate functions)
+ *  Maris Nartiss (refactoring for r.in.pdal)
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+#ifndef __BIN_WRITE_H__
+#define __BIN_WRITE_H__
+
+#include <grass/raster.h>
+
+/* forward declaration for point_binning.h */
+struct BinIndex;
+
+double get_sum(void *, void *, int, int, int, RASTER_MAP_TYPE);
+void write_sum(void *, void *, void *, int, int, RASTER_MAP_TYPE);
+void write_variance(void *, void *, void *,
+                    void *, int, int, RASTER_MAP_TYPE, int);
+void write_median(struct BinIndex *, void *,
+                  void *, int, int, RASTER_MAP_TYPE);
+void write_mode(struct BinIndex *, void *, void *, int, int);
+void write_percentile(struct BinIndex *, void *,
+                      void *, int, int, RASTER_MAP_TYPE, int);
+void write_skewness(struct BinIndex *, void *,
+                    void *, int, int, RASTER_MAP_TYPE);
+void write_trimmean(struct BinIndex *, void *,
+                    void *, int, int, RASTER_MAP_TYPE, double);
+void write_sidn(struct BinIndex *, void *, void *, int, int, int);
+void write_ev(struct BinIndex *, void *, void *,
+              int, int, RASTER_MAP_TYPE, int);
+
+
+#endif /* __BIN_WRITE_H__ */

raster/r.in.pdal/filters.c

@@ -0,0 +1,157 @@
+/*
+ * r.in.pdal filtering functions
+ * adapted from v.in.lidar
+ *
+ * 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)
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+
+#include "filters.h"
+
+#include "lidar.h"
+
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+
+int spatial_filter_from_option(struct Option *option, double *xmin,
+                               double *ymin, double *xmax, double *ymax)
+{
+    if (option->answer)
+        return FALSE;
+    int arg_s_num = 0;
+    int i = 0;
+
+    while (option->answers[i]) {
+        if (i == 0)
+            *xmin = atof(option->answers[i]);
+        if (i == 1)
+            *ymin = atof(option->answers[i]);
+        if (i == 2)
+            *xmax = atof(option->answers[i]);
+        if (i == 3)
+            *ymax = atof(option->answers[i]);
+        arg_s_num++;
+        i++;
+    }
+    if (arg_s_num != 4)
+        G_fatal_error(_("4 values required for '%s' option"), option->key);
+    return TRUE;
+}
+
+int spatial_filter_from_current_region(double *xmin, double *ymin,
+                                       double *xmax, double *ymax)
+{
+    struct Cell_head region;
+
+    G_get_window(&region);
+    *xmin = region.west;
+    *xmax = region.east;
+    *ymin = region.south;
+    *ymax = region.north;
+    return TRUE;
+}
+
+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)
+{
+    return_filter->filter = LAS_ALL;
+    if (name) {
+        if (strcmp(name, "first") == 0)
+            return_filter->filter = LAS_FIRST;
+        else if (strcmp(name, "last") == 0)
+            return_filter->filter = LAS_LAST;
+        else if (strcmp(name, "mid") == 0)
+            return_filter->filter = LAS_MID;
+        else
+            G_fatal_error(_("Unknown return filter value <%s>"), name);
+    }
+    if (return_filter->filter == LAS_ALL)
+        return FALSE;
+    else
+        return TRUE;
+}
+
+int return_filter_is_out(struct ReturnFilter *return_filter, int return_n,
+                         int n_returns)
+{
+    if (return_filter->filter == LAS_ALL)
+        return FALSE;
+    int skipme = 1;
+
+    switch (return_filter->filter) {
+    case LAS_FIRST:
+        if (return_n == 1)
+            skipme = 0;
+        break;
+    case LAS_MID:
+        if (return_n > 1 && return_n < n_returns)
+            skipme = 0;
+        break;
+    case LAS_LAST:
+        if (n_returns > 1 && return_n == n_returns)
+            skipme = 0;
+        break;
+    }
+    if (skipme)
+        return TRUE;
+    return FALSE;
+}
+
+int class_filter_create_from_strings(struct ClassFilter *class_filter,
+                                     char **classes)
+{
+    class_filter->str_classes = classes;
+    if (classes)
+        return TRUE;
+    else
+        return FALSE;
+}
+
+int class_filter_is_out(struct ClassFilter *class_filter, int class_n)
+{
+    if (!class_filter->str_classes)
+        return FALSE;
+    int i = 0;
+    int skipme = TRUE;
+
+    while (class_filter->str_classes[i]) {
+        if (class_n == atoi(class_filter->str_classes[i])) {
+            skipme = FALSE;
+            break;
+        }
+        i++;
+    }
+    if (skipme)
+        return TRUE;
+    return FALSE;
+}

raster/r.in.pdal/filters.h

@@ -0,0 +1,47 @@
+/*
+ * r.in.pdal filtering functions
+ * adapded from v.in.lidar
+ *
+ * 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)
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+#ifndef __FILTERS_H__
+#define __FILTERS_H__
+
+struct ReturnFilter
+{
+    int filter;
+};
+
+struct ClassFilter
+{
+
+    /** NULL terminated list of class numbers represented as string */
+    char **str_classes;
+};
+
+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 return_filter_create_from_string(struct ReturnFilter *return_filter,
+                                     const char *name);
+int return_filter_is_out(struct ReturnFilter *return_filter, int return_n,
+                         int n_returns);
+int class_filter_create_from_strings(struct ClassFilter *class_filter,
+                                     char **classes);
+int class_filter_is_out(struct ClassFilter *class_filter, int class_n);
+
+#endif /* __FILTERS_H__ */

raster/r.in.pdal/grasslidarfilter.cpp

@@ -0,0 +1,97 @@
+
+ /****************************************************************************
+ *
+ * MODULE:    r.in.pdal
+ *
+ * AUTHOR(S): Vaclav Petras
+ *
+ * PURPOSE:   Imports LAS LiDAR point clouds to a raster map using
+ *            aggregate statistics.
+ *
+ * COPYRIGHT: (C) 2011-2019 by Vaclav Petras and 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 "grasslidarfilter.h"
+
+
+bool GrassLidarFilter::processOne(pdal::PointRef & point)
+{
+    using pdal::Dimension::Id;
+
+    double x = point.getFieldAs < double >(Id::X);
+    double y = point.getFieldAs < double >(Id::Y);
+    double z = point.getFieldAs < double >(Id::Z);
+
+    n_processed_++;
+
+    z *= zscale_;
+
+    if (use_spatial_filter_) {
+        if (x < xmin_ || x > xmax_ || y < ymin_ || y > ymax_) {
+            n_outside_++;
+            return false;
+        }
+    }
+    if (use_irange_) {
+        double intensity = point.getFieldAs < double >(Id::Intensity);
+
+        intensity *= iscale_;
+        if (intensity < imin_ || intensity > imax_) {
+            irange_filtered_++;
+            return false;
+        }
+    }
+    if (use_drange_) {
+        double value = point.getFieldAs < double >(dim_to_import_);
+
+        value *= dscale_;
+        if (value < dmin_ || value > dmax_) {
+            drange_filtered_++;
+            return false;
+        }
+    }
+    if (base_segment_) {
+        double base_z;
+
+        if (rast_segment_get_value_xy(base_segment_, input_region_,
+                                      base_raster_data_type_, x, y, &base_z))
+            z -= base_z;
+        else
+            return false;       // skip points outside of base raster
+    }
+    if (use_zrange_) {
+        if (z < zmin_ || z > zmax_) {
+            zrange_filtered_++;
+            return false;
+        }
+    }
+    if (use_return_filter_) {
+        int return_n = point.getFieldAs < int >(Id::ReturnNumber);
+        int n_returns = point.getFieldAs < int >(Id::NumberOfReturns);
+
+        if (return_filter_is_out(&return_filter_, return_n, n_returns)) {
+            return_filtered_++;
+            return false;
+        }
+    }
+    if (use_class_filter_) {
+        int point_class = point.getFieldAs < int >(Id::Classification);
+
+        if (class_filter_is_out(&class_filter_, point_class)) {
+            n_class_filtered_++;
+            return false;
+        }
+    }
+    n_passed_++;
+
+    // PDAL filter streaming function should return false when
+    // point is filtered out (and should not be used by next stage)
+    // and true otherwise.
+    return true;
+}

raster/r.in.pdal/grasslidarfilter.h

@@ -0,0 +1,233 @@
+
+ /****************************************************************************
+ *
+ * MODULE:    r.in.pdal
+ *
+ * AUTHOR(S): Vaclav Petras
+ *
+ * PURPOSE:   Imports LAS LiDAR point clouds to a raster map using
+ *            aggregate statistics.
+ *
+ * COPYRIGHT: (C) 2011-2019 by Vaclav Petras and 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 GRASSLIDARFILTER_H
+#define GRASSLIDARFILTER_H
+
+extern "C"
+{
+#include "filters.h"
+#include "lidar.h"
+#include "rast_segment.h"
+}
+
+#include <pdal/Filter.hpp>
+#include <pdal/Streamable.hpp>
+#include <pdal/Dimension.hpp>
+
+
+/* All GRASS GIS filters which are similar across multiple modules
+ * put together as one PDAL Stage class.
+ */
+class GrassLidarFilter:public pdal::Filter, public pdal::Streamable
+{
+  public:
+    GrassLidarFilter():
+        dim_to_import_(pdal::Dimension::Id::Z),
+        use_spatial_filter_(false),
+        use_zrange_(false),
+        use_irange_(false),
+        use_drange_(false),
+        xmin_(0),
+        xmax_(0),
+        ymin_(0),
+        ymax_(0),
+        zmin_(0),
+        zmax_(0),
+        imin_(0),
+        imax_(0),
+        dmin_(0),
+        dmax_(0),
+        n_processed_(0),
+        n_passed_(0),
+        n_outside_(0),
+        zrange_filtered_(0),
+        irange_filtered_(0),
+        drange_filtered_(0),
+        return_filtered_(0),
+        n_class_filtered_(0),
+        zscale_(1),
+        iscale_(1),
+        dscale_(1),
+        use_class_filter_(false),
+        use_return_filter_(false), base_segment_(nullptr)
+    {
+    }
+    std::string getName() const
+    {
+        return "filters.grasslidar";
+    }
+
+    void dim_to_import(pdal::Dimension::Id dim_to_import)
+    {
+        dim_to_import_ = dim_to_import;
+    }
+
+    void set_spatial_filter(double xmin, double xmax,
+                            double ymin, double ymax)
+    {
+        use_spatial_filter_ = true;
+        xmin_ = xmin;
+        xmax_ = xmax;
+        ymin_ = ymin;
+        ymax_ = ymax;
+        n_outside_ = 0;
+    }
+    void set_zrange_filter(double min, double max)
+    {
+        use_zrange_ = true;
+        zmin_ = min;
+        zmax_ = max;
+        zrange_filtered_ = 0;
+    }
+    void set_irange_filter(double min, double max)
+    {
+        use_irange_ = true;
+        imin_ = min;
+        imax_ = max;
+        irange_filtered_ = 0;
+    }
+    void set_drange_filter(double min, double max)
+    {
+        use_drange_ = true;
+        dmin_ = min;
+        dmax_ = max;
+        drange_filtered_ = 0;
+    }
+    void set_return_filter(ReturnFilter return_filter)
+    {
+        use_return_filter_ = true;
+        return_filter_ = return_filter;
+        return_filtered_ = 0;
+    }
+    void set_class_filter(ClassFilter class_filter)
+    {
+        use_class_filter_ = true;
+        class_filter_ = class_filter;
+        n_class_filtered_ = 0;
+    }
+    void set_base_raster(SEGMENT * base_segment,
+                         struct Cell_head *region, RASTER_MAP_TYPE rtype)
+    {
+        base_segment_ = base_segment;
+        input_region_ = region;
+        base_raster_data_type_ = rtype;
+    }
+    void set_z_scale(double scale)
+    {
+        zscale_ = scale;
+    }
+    void set_intensity_scale(double scale)
+    {
+        iscale_ = scale;
+    }
+    void set_d_scale(double scale)
+    {
+        dscale_ = scale;
+    }
+
+    gpoint_count num_processed()
+    {
+        return n_processed_;
+    }
+    gpoint_count num_passed()
+    {
+        return n_passed_;
+    }
+    gpoint_count num_return_filtered()
+    {
+        return return_filtered_;
+    }
+    gpoint_count num_class_filtered()
+    {
+        return n_class_filtered_;
+    }
+    gpoint_count num_zrange_filtered()
+    {
+        return zrange_filtered_;
+    }
+    gpoint_count num_irange_filtered()
+    {
+        return irange_filtered_;
+    }
+    gpoint_count num_drange_filtered()
+    {
+        return drange_filtered_;
+    }
+    gpoint_count num_spatially_filtered()
+    {
+        return n_outside_;
+    }
+
+  private:
+    virtual void filter(pdal::PointView & view)
+    {
+        pdal::PointRef p(view, 0);
+        for (pdal::PointId idx = 0; idx < view.size(); ++idx) {
+            p.setPointId(idx);
+            processOne(p);
+        }
+    }
+    virtual bool processOne(pdal::PointRef & point);
+
+    pdal::Dimension::Id dim_to_import_;
+
+    bool use_spatial_filter_;
+    bool use_zrange_;
+    bool use_irange_;
+    bool use_drange_;
+    double xmin_;
+    double xmax_;
+    double ymin_;
+    double ymax_;
+    double zmin_;
+    double zmax_;
+    double imin_;
+    double imax_;
+    double dmin_;
+    double dmax_;
+    gpoint_count n_processed_;
+    gpoint_count n_passed_;
+    gpoint_count n_outside_;
+    gpoint_count zrange_filtered_;
+    gpoint_count irange_filtered_;
+    gpoint_count drange_filtered_;
+    gpoint_count return_filtered_;
+    gpoint_count n_class_filtered_;
+
+    double zscale_;
+    double iscale_;
+    double dscale_;
+
+    bool use_class_filter_;
+    ClassFilter class_filter_;
+    bool use_return_filter_;
+    ReturnFilter return_filter_;
+
+    SEGMENT *base_segment_;
+    struct Cell_head *input_region_;
+    RASTER_MAP_TYPE base_raster_data_type_;
+
+    // not implemented
+    GrassLidarFilter & operator=(const GrassLidarFilter &);
+    GrassLidarFilter(const GrassLidarFilter &);
+};
+
+
+#endif // GRASSLIDARFILTER_H

raster/r.in.pdal/grassrasterwriter.h

@@ -0,0 +1,143 @@
+
+/***************************************************************************
+ *
+ * MODULE:    r.in.pdal
+ *
+ * AUTHOR(S): Vaclav Petras
+ *
+ * PURPOSE:   Imports LAS LiDAR point clouds to a raster map using
+ *            aggregate statistics.
+ *
+ * COPYRIGHT: (C) 2019 by Vaclav Petras and 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 GRASSRASTERWRITER_H
+#define GRASSRASTERWRITER_H
+
+extern "C"
+{
+#include "lidar.h"
+#include "point_binning.h"
+}
+
+extern "C"
+{
+#include <grass/gis.h>
+#include <grass/raster.h>
+}
+
+#include <pdal/Streamable.hpp>
+#include <pdal/Writer.hpp>
+
+
+/* Binning code wrapped as a PDAL Writer class */
+class GrassRasterWriter:public pdal::Writer, public pdal::Streamable
+{
+  public:
+    GrassRasterWriter():n_processed(0)
+    {
+    }
+
+    std::string getName() const
+    {
+        return "writers.grassbinning";
+    }
+
+    void set_binning(struct Cell_head *region,
+                     struct PointBinning *point_binning,
+                     struct BinIndex *bin_index_nodes,
+                     RASTER_MAP_TYPE rtype, int cols)
+    {
+        region_ = region;
+        point_binning_ = point_binning;
+        bin_index_nodes_ = bin_index_nodes;
+        rtype_ = rtype;
+        cols_ = cols;
+    }
+
+    void dim_to_import(pdal::Dimension::Id dim_to_import)
+    {
+        dim_to_import_ = dim_to_import;
+    }
+
+    void set_output_scale(double scale)
+    {
+        scale_ = scale;
+    }
+
+    void set_base_raster(SEGMENT * base_segment,
+                         struct Cell_head *region, RASTER_MAP_TYPE rtype)
+    {
+        base_segment_ = base_segment;
+        input_region_ = region;
+        base_raster_data_type_ = rtype;
+    }
+
+    virtual void write(const pdal::PointViewPtr view)
+    {
+        pdal::PointRef p(*view, 0);
+        for (pdal::PointId idx = 0; idx < view->size(); ++idx) {
+            p.setPointId(idx);
+            processOne(p);
+        }
+    }
+
+    virtual bool processOne(pdal::PointRef & point)
+    {
+        using namespace pdal::Dimension;
+
+        double x = point.getFieldAs < double >(Id::X);
+        double y = point.getFieldAs < double >(Id::Y);
+        double z = point.getFieldAs < double >(dim_to_import_);
+
+        z *= scale_;
+        if (base_segment_) {
+            double base_z;
+
+            rast_segment_get_value_xy(base_segment_, input_region_,
+                                      base_raster_data_type_, x, y, &base_z);
+            z -= base_z;
+        }
+
+        // TODO: check the bounds and report discrepancies in
+        // number of filtered out vs processed to the user
+        // (alternativelly, change the spatial bounds test to
+        // give same results as this, but it might be actually helpful
+        // to tell user that they have points right on the border)
+        int arr_row = (int)((region_->north - y) / region_->ns_res);
+        int arr_col = (int)((x - region_->west) / region_->ew_res);
+
+        if (arr_row >= region_->rows || arr_col >= region_->cols) {
+            G_message(_("A point on the edge of computational region detected. Ignoring."));
+            return false;
+        }
+
+        update_value(point_binning_, bin_index_nodes_, cols_,
+                     arr_row, arr_col, rtype_, x, y, z);
+        n_processed++;
+        return true;
+    }
+
+    gpoint_count n_processed;
+
+  private:
+    struct Cell_head *region_;
+    struct PointBinning *point_binning_;
+    struct BinIndex *bin_index_nodes_;
+    RASTER_MAP_TYPE rtype_;
+    int cols_;
+    double scale_;
+
+    pdal::Dimension::Id dim_to_import_;
+
+    SEGMENT *base_segment_;
+    struct Cell_head *input_region_;
+    RASTER_MAP_TYPE base_raster_data_type_;
+};
+
+#endif // GRASSRASTERWRITER_H

raster/r.in.pdal/info.cpp

@@ -0,0 +1,120 @@
+/*
+ * r.in.pdal Functions printing out various information on input LAS files
+ *  
+ *   Copyright 2021 by Maris Nartiss, and The GRASS Development Team
+ *   Author: Maris Nartiss
+ *
+ *   This program is free software licensed under the GPL (>=v2).
+ *   Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+#include "info.h"
+
+
+void get_extent(struct StringList *infiles, double *min_x, double *max_x,
+                double *min_y, double *max_y, double *min_z, double *max_z)
+{
+    pdal::StageFactory factory;
+    bool first = 1;
+
+    *min_x = *max_x = *min_y = *max_y = *min_z = *max_z = 0.0 / 0.0;
+
+    for (int i = 0; i < infiles->num_items; i++) {
+        const char *infile = infiles->items[i];
+
+        std::string pdal_read_driver = factory.inferReaderDriver(infile);
+        if (pdal_read_driver.empty())
+            G_fatal_error("Cannot determine input file type of <%s>", infile);
+
+        pdal::PointTable table;
+        pdal::Options las_opts;
+        pdal::Option las_opt("filename", infile);
+        las_opts.add(las_opt);
+        pdal::LasReader las_reader;
+        las_reader.setOptions(las_opts);
+        las_reader.prepare(table);
+        pdal::LasHeader las_header = las_reader.header();
+        if (first) {
+            *min_x = las_header.minX();
+            *min_y = las_header.minY();
+            *min_z = las_header.minZ();
+            *max_x = las_header.maxX();
+            *max_y = las_header.maxY();
+            *max_z = las_header.maxZ();
+
+            first = 0;
+        }
+        else {
+            if (*min_x > las_header.minX())
+                *min_x = las_header.minX();
+            if (*min_y > las_header.minY())
+                *min_y = las_header.minY();
+            if (*min_z > las_header.minZ())
+                *min_z = las_header.minZ();
+            if (*max_x < las_header.maxX())
+                *max_x = las_header.maxX();
+            if (*max_y < las_header.maxY())
+                *max_y = las_header.maxY();
+            if (*max_z < las_header.maxZ())
+                *max_z = las_header.maxZ();
+        }
+    }
+}
+
+
+void print_extent(struct StringList *infiles)
+{
+    double min_x, max_x, min_y, max_y, min_z, max_z;
+
+    get_extent(infiles, &min_x, &max_x, &min_y, &max_y, &min_z, &max_z);
+    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, max_z);
+}
+
+
+void print_lasinfo(struct StringList *infiles)
+{
+    pdal::StageFactory factory;
+    pdal::MetadataNode meta_node;
+
+    std::cout << std::endl << "Using PDAL library version '" <<
+        pdal::Config::fullVersionString() << "'" << std::endl << std::endl;
+
+    for (int i = 0; i < infiles->num_items; i++) {
+        const char *infile = infiles->items[i];
+
+        std::string pdal_read_driver = factory.inferReaderDriver(infile);
+        if (pdal_read_driver.empty())
+            G_fatal_error("Cannot determine input file type of <%s>", infile);
+
+        pdal::PointTable table;
+        pdal::Options las_opts;
+        pdal::Option las_opt("filename", infile);
+        las_opts.add(las_opt);
+        pdal::LasReader las_reader;
+        las_reader.setOptions(las_opts);
+        las_reader.prepare(table);
+        pdal::LasHeader las_header = las_reader.header();
+        pdal::PointLayoutPtr point_layout = table.layout();
+        const pdal::Dimension::IdList & dims = point_layout->dims();
+
+        std::cout << "File: " << infile << std::endl;
+        std::cout << las_header;
+
+        bool first = 1;
+
+        for (auto di = dims.begin(); di != dims.end(); ++di) {
+            pdal::Dimension::Id d = *di;
+
+            if (first) {
+                std::cout << "Dimensions: " << point_layout->dimName(d);
+                first = 0;
+            }
+            else {
+                std::cout << ", " << point_layout->dimName(d);
+            }
+        }
+        std::cout << std::endl << std::endl;
+    }
+}

raster/r.in.pdal/info.h

@@ -0,0 +1,34 @@
+/*
+ * r.in.pdal Functions printing out various information on input LAS files
+ *  
+ *   Copyright 2021 by Maris Nartiss, and The GRASS Development Team
+ *   Author: Maris Nartiss
+ *
+ *   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 <pdal/pdal_config.hpp>
+#include <pdal/PointTable.hpp>
+#include <pdal/StageFactory.hpp>
+#include <pdal/Options.hpp>
+#include <pdal/io/LasReader.hpp>
+#include <pdal/io/LasHeader.hpp>
+
+extern "C"
+{
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "string_list.h"
+}
+
+void get_extent(struct StringList *, double *, double *,
+                double *, double *, double *, double *);
+void print_extent(struct StringList *);
+void print_lasinfo(struct StringList *);
+
+#endif // INFO_H

raster/r.in.pdal/lidar.c

@@ -0,0 +1,60 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.in.pdal
+ *               adapted from v.in.lidar
+ * AUTHOR(S):    Vaclav Petras
+ * PURPOSE:      common lidar-related definitions
+ * COPYRIGHT:    (C) 2015 by the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *               License (>=v2). Read the COPYING file that comes with GRASS
+ *               for details.
+ *
+ *****************************************************************************/
+
+
+#include "lidar.h"
+
+void GLidarLayers_set_no_layers(struct GLidarLayers *layers)
+{
+    layers->id_layer = 0;
+    layers->return_layer = 0;
+    layers->class_layer = 0;
+    layers->rgb_layer = 0;
+}
+
+struct class_table class_val[] = {
+    {0, "Created, never classified"},
+    {1, "Unclassified"},
+    {2, "Ground"},
+    {3, "Low Vegetation"},
+    {4, "Medium Vegetation"},
+    {5, "High Vegetation"},
+    {6, "Building"},
+    {7, "Low Point (noise)"},
+    {8, "Model Key-point (mass point)"},
+    {9, "Water"},
+    {10, "Reserved for ASPRS Definition"},
+    {11, "Reserved for ASPRS Definition"},
+    {12, "Overlap Points"},
+    {13 /* 13 - 31 */ , "Reserved for ASPRS Definition"},
+    {0, 0}
+};
+
+struct class_table class_type[] = {
+    {5, "Synthetic"},
+    {6, "Key-point"},
+    {7, "Withheld"},
+    {0, 0}
+};
+
+int return_to_cat(int return_n, int n_returns)
+{
+    if (return_n == 1)
+        return LAS_FIRST;
+    else if (n_returns > 1 && return_n == n_returns)
+        return LAS_LAST;
+    else
+        return LAS_MID;
+}

raster/r.in.pdal/lidar.h

@@ -0,0 +1,93 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.in.pdal
+ *               adapted from v.in.lidar
+ * AUTHOR(S):    Vaclav Petras
+ * PURPOSE:      common lidar-related definitions
+ * COPYRIGHT:    (C) 2015 by the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *               License (>=v2). Read the COPYING file that comes with GRASS
+ *               for details.
+ *
+ *****************************************************************************/
+
+
+#ifndef GRASS_LIDAR_H
+#define GRASS_LIDAR_H
+
+#define LAS_ALL 0
+#define LAS_FIRST 1
+#define LAS_MID 2
+#define LAS_LAST 3
+
+/* Type and format specifier for count of points */
+#ifdef HAVE_LONG_LONG_INT
+typedef unsigned long long gpoint_count;
+
+#define GPOINT_COUNT_FORMAT "%llu"
+#else
+typedef unsigned long gpoint_count;
+
+#define GPOINT_COUNT_FORMAT "%lu"
+#endif
+
+struct GLidarLayers
+{
+    int id_layer;
+    int return_layer;
+    int class_layer;
+    int rgb_layer;
+};
+
+void GLidarLayers_set_no_layers(struct GLidarLayers *layers);
+
+/*
+ * ASPRS Standard LIDAR Point Classes
+ * Classification Value (bits 0:4) : Meaning
+ *      0 : Created, never classified
+ *      1 : Unclassified
+ *      2 : Ground
+ *      3 : Low Vegetation
+ *      4 : Medium Vegetation
+ *      5 : High Vegetation
+ *      6 : Building
+ *      7 : Low Point (noise)
+ *      8 : Model Key-point (mass point)
+ *      9 : Water
+ *     10 : Reserved for ASPRS Definition
+ *     11 : Reserved for ASPRS Definition
+ *     12 : Overlap Points
+ *  13-31 : Reserved for ASPRS Definition
+ */
+
+/* Classification Bit Field Encoding
+ * Bits | Field Name     | Description
+ *  0-4 | Classification | Standard ASPRS classification as defined in the
+ *                         above classification table.
+ *    5 | Synthetic      | If set then this point was created by a technique
+ *                         other than LIDAR collection such as digitized from
+ *                         a photogrammetric stereo model or by traversing
+ *                         a waveform.
+ *    6 | Key-point      | If set, this point is considered to be a model
+ *                         key-point and thus generally should not be withheld
+ *                         in a thinning algorithm.
+ *    7 | Withheld       | If set, this point should not be included in
+ *                         processing (synonymous with Deleted).
+ */
+
+/* keep the comments above in sync with the .c file */
+
+struct class_table
+{
+    int code;
+    char *name;
+};
+
+extern struct class_table class_val[];
+extern struct class_table class_type[];
+
+int return_to_cat(int return_n, int n_returns);
+
+#endif /* GRASS_LIDAR_H */

raster/r.in.pdal/main.cpp

@@ -0,0 +1,937 @@
+
+ /****************************************************************************
+ *
+ * MODULE:    r.in.pdal
+ *
+ * AUTHOR(S): Vaclav Petras
+ *            Based on r.in.xyz and r.in.lidar by Markus Metz,
+ *            Hamish Bowman, Volker Wichmann
+ *
+ * PURPOSE:   Imports LAS LiDAR point clouds to a raster map using
+ *            aggregate statistics.
+ *
+ * COPYRIGHT: (C) 2019-2021 by Vaclav Petras and 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 <pdal/PointTable.hpp>
+#include <pdal/PointLayout.hpp>
+#include <pdal/StageFactory.hpp>
+#include <pdal/io/LasReader.hpp>
+#include <pdal/io/LasHeader.hpp>
+#include <pdal/Options.hpp>
+#include <pdal/filters/MergeFilter.hpp>
+#include <pdal/filters/ReprojectionFilter.hpp>
+#include <pdal/filters/StreamCallbackFilter.hpp>
+
+extern "C"
+{
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/gprojects.h>
+#include <grass/glocale.h>
+}
+
+#include "grasslidarfilter.h"
+#include "grassrasterwriter.h"
+#include "info.h"
+
+extern "C"
+{
+#include "lidar.h"
+#include "projection.h"
+#include "filters.h"
+#include "point_binning.h"
+#include "string_list.h"
+}
+
+#define BUFFSIZE GPATH_MAX
+
+
+int main(int argc, char *argv[])
+{
+    int out_fd;
+    char *outmap;
+
+    RASTER_MAP_TYPE rtype, base_raster_data_type;
+    struct History history;
+    char title[64];
+    SEGMENT base_segment;
+    struct PointBinning point_binning;
+    void *raster_row;
+    struct Cell_head region;
+    struct Cell_head input_region;
+    int rows, cols;             /* scan box size */
+
+    char buff[BUFFSIZE];
+
+    struct BinIndex bin_index_nodes;
+
+    bin_index_nodes.num_nodes = 0;
+    bin_index_nodes.max_nodes = 0;
+    bin_index_nodes.nodes = NULL;
+
+    struct Cell_head loc_wind;
+
+    G_gisinit(argv[0]);
+
+    GModule *module = G_define_module();
+
+    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.");
+
+    Option *input_opt = G_define_standard_option(G_OPT_F_BIN_INPUT);
+
+    input_opt->required = NO;
+    input_opt->label = _("LAS input file");
+    input_opt->description =
+        _("LiDAR input files in LAS format (*.las or *.laz)");
+    input_opt->guisection = _("Input");
+
+    Option *output_opt = G_define_standard_option(G_OPT_R_OUTPUT);
+
+    output_opt->required = NO;
+    output_opt->guisection = _("Output");
+
+    Option *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");
+
+    Option *method_opt = G_define_option();
+
+    method_opt->key = "method";
+    method_opt->type = TYPE_STRING;
+    method_opt->required = NO;
+    method_opt->description = _("Statistic to use for raster values");
+    method_opt->options =
+        "n,min,max,range,sum,mean,stddev,variance,coeff_var,median,mode,"
+        "percentile,skewness,trimmean,sidnmax,sidnmin,ev1,ev2,ev3";
+    method_opt->answer = const_cast < char *>("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;"
+               "mode;%s;"
+               "percentile;%s;"
+               "skewness;%s;"
+               "trimmean;%s;"
+               "sidnmax;%s;"
+               "sidnmin;%s;"
+               "ev1;%s;"
+               "ev2;%s;"
+               "ev3;%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"),
+               _("Mode 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"),
+               _("Maximum number of points in cell per source ID"),
+               _("Minimum number of points in cell per source ID"),
+               _("First eigenvalue of point x, y, z coordinates"),
+               _("Second eigenvalue of point x, y, z coordinates"),
+               _("Third eigenvalue of point x, y, z coordinates"));
+
+    Option *type_opt = G_define_standard_option(G_OPT_R_TYPE);
+
+    type_opt->required = NO;
+    type_opt->answer = const_cast < char *>("FCELL");
+
+    Option *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 range filter for"
+          " Z afterwards");
+    base_raster_opt->guisection = _("Transform");
+
+    Option *zrange_opt = G_define_option();
+
+    zrange_opt->key = "zrange";
+    zrange_opt->type = TYPE_DOUBLE;
+    zrange_opt->required = NO;
+    zrange_opt->key_desc = "min,max";
+    zrange_opt->label = _("Filter range for Z data (min,max)");
+    zrange_opt->description =
+        _("Applied after base_raster transformation step");
+    zrange_opt->guisection = _("Selection");
+
+    Option *zscale_opt = G_define_option();
+
+    zscale_opt->key = "zscale";
+    zscale_opt->type = TYPE_DOUBLE;
+    zscale_opt->required = NO;
+    zscale_opt->answer = const_cast < char *>("1.0");
+
+    zscale_opt->description = _("Scale to apply to Z data");
+    zscale_opt->guisection = _("Transform");
+
+    Option *irange_opt = G_define_option();
+
+    irange_opt->key = "irange";
+    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");
+
+    Option *iscale_opt = G_define_option();
+
+    iscale_opt->key = "iscale";
+    iscale_opt->type = TYPE_DOUBLE;
+    iscale_opt->required = NO;
+    iscale_opt->answer = const_cast < char *>("1.0");
+
+    iscale_opt->description = _("Scale to apply to intensity values");
+    iscale_opt->guisection = _("Transform");
+
+    Option *drange_opt = G_define_option();
+
+    drange_opt->key = "drange";
+    drange_opt->type = TYPE_DOUBLE;
+    drange_opt->required = NO;
+    drange_opt->key_desc = "min,max";
+    drange_opt->description =
+        _("Filter range for output dimension values (min,max)");
+    drange_opt->guisection = _("Selection");
+
+    Option *dscale_opt = G_define_option();
+
+    dscale_opt->key = "dscale";
+    dscale_opt->type = TYPE_DOUBLE;
+    dscale_opt->required = NO;
+    dscale_opt->label = _("Scale to apply to output dimension values");
+    dscale_opt->description =
+        _("Use if output dimension is not Z or intensity");
+    dscale_opt->guisection = _("Transform");
+
+    Flag *reproject_flag = G_define_flag();
+
+    reproject_flag->key = 'w';
+    reproject_flag->label =
+        _("Reproject to location's coordinate system if needed");
+    reproject_flag->description =
+        _("Reprojects input dataset to the coordinate system of"
+          " the GRASS location (by default only datasets with the"
+          " matching cordinate system can be imported");
+    reproject_flag->guisection = _("Projection");
+
+    // TODO: from the API it seems that also prj file path and proj string will work
+    Option *input_srs_opt = G_define_option();
+
+    input_srs_opt->key = "input_srs";
+    input_srs_opt->type = TYPE_STRING;
+    input_srs_opt->required = NO;
+    input_srs_opt->label =
+        _("Input dataset projection (WKT or EPSG, e.g. EPSG:4326)");
+    input_srs_opt->description =
+        _("Override input dataset coordinate system using EPSG code"
+          " or WKT definition");
+    input_srs_opt->guisection = _("Projection");
+
+    /* I would prefer to call the following "percentile", but that has too
+     * much namespace overlap with the "percent" option above */
+    Option *pth_opt = G_define_option();
+
+    pth_opt->key = "pth";
+    pth_opt->type = TYPE_INTEGER;
+    pth_opt->required = NO;
+    pth_opt->options = "1-100";
+    pth_opt->description = _("pth percentile of the values");
+    pth_opt->guisection = _("Statistic");
+
+    Option *trim_opt = G_define_option();
+
+    trim_opt->key = "trim";
+    trim_opt->type = TYPE_DOUBLE;
+    trim_opt->required = NO;
+    trim_opt->options = "0-50";
+    trim_opt->label =
+        _("Discard given percentage of the smallest and largest values");
+    trim_opt->description =
+        _("Discard <trim> percent of the smallest and <trim> percent of the largest observations");
+    trim_opt->guisection = _("Statistic");
+
+    Option *res_opt = G_define_option();
+
+    res_opt->key = "resolution";
+    res_opt->type = TYPE_DOUBLE;
+    res_opt->required = NO;
+    res_opt->description = _("Output raster resolution");
+    res_opt->guisection = _("Output");
+
+    Option *return_filter_opt = G_define_option();
+
+    return_filter_opt->key = "return_filter";
+    return_filter_opt->type = TYPE_STRING;
+    return_filter_opt->required = NO;
+    return_filter_opt->label =
+        _("Only import points of selected return type");
+    return_filter_opt->description =
+        _("If not specified, all points are imported");
+    return_filter_opt->options = "first,last,mid";
+    return_filter_opt->guisection = _("Selection");
+
+    Option *class_filter_opt = G_define_option();
+
+    class_filter_opt->key = "class_filter";
+    class_filter_opt->type = TYPE_INTEGER;
+    class_filter_opt->multiple = YES;
+    class_filter_opt->required = NO;
+    class_filter_opt->label = _("Only import points of selected class(es)");
+    class_filter_opt->description = _("Input is comma separated integers. "
+                                      "If not specified, all points are imported.");
+    class_filter_opt->guisection = _("Selection");
+
+    Option *dimension_opt = G_define_option();
+
+    dimension_opt->key = "dimension";
+    dimension_opt->type = TYPE_STRING;
+    dimension_opt->required = NO;
+    dimension_opt->label = _("Dimension (variable) to use for raster values");
+    dimension_opt->options =
+        "z,intensity,number,returns,direction,angle,class,source";
+    dimension_opt->answer = const_cast < char *>("z");
+
+    dimension_opt->guisection = _("Selection");
+    G_asprintf((char **)&(dimension_opt->descriptions),
+               "z;%s;"
+               "intensity;%s;"
+               "number;%s;"
+               "returns;%s;"
+               "direction;%s;"
+               "angle;%s;" "class;%s;" "source;%s", _("Z coordinate"),
+               /* GTC: LAS LiDAR point property */
+               _("Intensity"),
+               /* GTC: LAS LiDAR point property */
+               _("Return number"),
+               /* GTC: LAS LiDAR point property */
+               _("Number of returns"),
+               /* GTC: LAS LiDAR point property */
+               _("Scan direction"),
+               /* GTC: LAS LiDAR point property */
+               _("Scan angle"),
+               /* GTC: LAS LiDAR point property */
+               _("Point class value"),
+               /* GTC: LAS LiDAR point property */
+               _("Source ID"));
+
+    Option *user_dimension_opt = G_define_option();
+
+    user_dimension_opt->key = "user_dimension";
+    user_dimension_opt->type = TYPE_STRING;
+    user_dimension_opt->required = NO;
+    user_dimension_opt->label =
+        _("Custom dimension (variable) to use for raster values");
+    user_dimension_opt->description = _("PDAL dimension name");
+    user_dimension_opt->guisection = _("Selection");
+
+    Flag *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 =
+        _("Set internally computational region extents based on the"
+          " point cloud");
+    extents_flag->guisection = _("Output");
+
+    Flag *set_region_flag = G_define_flag();
+
+    set_region_flag->key = 'n';
+    set_region_flag->label =
+        _("Set computation region to match the new raster map");
+    set_region_flag->description =
+        _("Set computation region to match the 2D extent and resolution"
+          " of the newly created new raster map");
+    set_region_flag->guisection = _("Output");
+
+    Flag *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");
+    over_flag->guisection = _("Projection");
+
+    Flag *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 =
+        _("For getting values from base raster, use its actual"
+          " resolution instead of computational region resolution");
+    base_rast_res_flag->guisection = _("Transform");
+
+    Flag *print_info_flag = G_define_flag();
+
+    print_info_flag->key = 'p';
+    print_info_flag->description = _("Print LAS file info and exit");
+
+    Flag *print_extent_flag = G_define_flag();
+
+    print_extent_flag->key = 'g';
+    print_extent_flag->description =
+        _("Print data file extent in shell script style and then exit");
+
+
+    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);
+    G_option_exclusive(base_rast_res_flag, res_opt, NULL);
+    G_option_exclusive(reproject_flag, over_flag, NULL);
+    G_option_required(output_opt, print_extent_flag, print_info_flag, NULL);
+
+
+    if (G_parser(argc, argv))
+        return EXIT_FAILURE;
+
+    /* Get input file list. Needs to be done before printing extent. */
+    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);
+    }
+
+    /* If we print extent, there is no need to validate rest of the input */
+    if (print_extent_flag->answer) {
+        print_extent(&infiles);
+        exit(EXIT_SUCCESS);
+    }
+
+    if (print_info_flag->answer) {
+        print_lasinfo(&infiles);
+        exit(EXIT_SUCCESS);
+    }
+
+    /* we could use rules but this gives more info and allows continuing */
+    if (set_region_flag->answer &&
+        !(extents_flag->answer || res_opt->answer ||
+          base_rast_res_flag->answer)) {
+        G_warning(_("Flag %c makes sense only with %s option or -%c flag or -%c flag"),
+                  set_region_flag->key, res_opt->key, extents_flag->key,
+                  base_rast_res_flag->key);
+        /* avoid the call later on */
+        set_region_flag->answer = '\0';
+    }
+
+    /* Trim option is used only for trimmean method */
+    if (trim_opt->answer != NULL &&
+        strcmp(method_opt->answer, "trimmean") != 0) {
+        G_fatal_error(_("Trim option can be used only with trimmean method"));
+    }
+
+    /* Point density counting does not require any dimension information */
+    if ((strcmp(method_opt->answer, "sidnmax") == 0 ||
+         strcmp(method_opt->answer, "sidnmin") == 0 ||
+         strcmp(method_opt->answer, "n") == 0 ||
+         strcmp(method_opt->answer, "ev1") == 0 ||
+         strcmp(method_opt->answer, "ev2") == 0 ||
+         strcmp(method_opt->answer, "ev3") == 0) &&
+        (user_dimension_opt->answer ||
+         !(strcmp(dimension_opt->answer, "z") == 0)))
+        G_warning(_("Binning methods 'n', 'sidnmax', 'sidnmin' and "
+                    "eigenvalues are ignoring specified dimension"));
+
+    /* parse input values */
+    outmap = output_opt->answer;
+    if (input_opt->answer && access(input_opt->answer, F_OK) != 0) {
+        G_fatal_error(_("Input file <%s> does not exist"), input_opt->answer);
+    }
+
+    /* Set up input extent for point spatial filter */
+    double xmin = 0;
+    double ymin = 0;
+    double xmax = 0;
+    double ymax = 0;
+    bool use_spatial_filter = false;
+
+    Rast_get_window(&region);
+    /* G_get_window seems to be unreliable if the location has been changed */
+    G_get_set_window(&loc_wind);        /* TODO: v.in.lidar uses G_get_default_window() */
+
+    /* Region is set based on whole point cloud that could be larger than imported part */
+    if (extents_flag->answer) {
+        double min_x, max_x, min_y, max_y, min_z, max_z;
+
+        get_extent(&infiles, &min_x, &max_x, &min_y, &max_y, &min_z, &max_z);
+
+        region.east = xmax = max_x;
+        region.west = xmin = min_x;
+        region.north = ymax = max_y;
+        region.south = ymin = min_y;
+
+        use_spatial_filter = true;
+    }
+
+    /* Set up filtering options */
+    if (!extents_flag->answer) {
+        use_spatial_filter = spatial_filter_from_current_region(&xmin,
+                                                                &ymin,
+                                                                &xmax, &ymax);
+    }
+
+    double zrange_min, zrange_max;
+    bool use_zrange = range_filter_from_option(zrange_opt, &zrange_min,
+                                               &zrange_max);
+    double irange_min, irange_max;
+    bool use_irange = range_filter_from_option(irange_opt, &irange_min,
+                                               &irange_max);
+    double drange_min, drange_max;
+    bool use_drange = range_filter_from_option(drange_opt, &drange_min,
+                                               &drange_max);
+    struct ReturnFilter return_filter_struct;
+    bool use_return_filter =
+        return_filter_create_from_string(&return_filter_struct,
+                                         return_filter_opt->answer);
+    struct ClassFilter class_filter;
+    bool use_class_filter = class_filter_create_from_strings(&class_filter,
+                                                             class_filter_opt->
+                                                             answers);
+
+    point_binning_set(&point_binning, method_opt->answer, pth_opt->answer,
+                      trim_opt->answer);
+
+    /* Set up output map type */
+    if (strcmp("CELL", type_opt->answer) == 0)
+        rtype = CELL_TYPE;
+    else if (strcmp("DCELL", type_opt->answer) == 0)
+        rtype = DCELL_TYPE;
+    else
+        rtype = FCELL_TYPE;
+
+    if (point_binning.method == METHOD_N ||
+        point_binning.method == METHOD_MODE ||
+        point_binning.method == METHOD_SIDNMAX ||
+        point_binning.method == METHOD_SIDNMIN) {
+        if (rtype != CELL_TYPE)
+            G_warning(_("Output map type set to CELL"));
+        rtype = CELL_TYPE;
+    }
+
+    /* Set up output dimension */
+    // we use full qualification because the dim ns contains too general names
+    pdal::Dimension::Id dim_to_import = pdal::Dimension::Id::Z;
+
+    if (!user_dimension_opt->answer &&
+        !(strcmp(dimension_opt->answer, "z") == 0)) {
+        /* Should we enfocte the CELL type? */
+        if (rtype != CELL_TYPE)
+            G_warning(_("Output map type set to CELL"));
+        rtype = CELL_TYPE;
+
+        if (strcmp(dimension_opt->answer, "intensity") == 0) {
+            dim_to_import = pdal::Dimension::Id::Intensity;
+        }
+        else if (strcmp(dimension_opt->answer, "number") == 0) {
+            dim_to_import = pdal::Dimension::Id::ReturnNumber;
+        }
+        else if (strcmp(dimension_opt->answer, "returns") == 0) {
+            dim_to_import = pdal::Dimension::Id::NumberOfReturns;
+        }
+        else if (strcmp(dimension_opt->answer, "direction") == 0) {
+            dim_to_import = pdal::Dimension::Id::ScanDirectionFlag;
+        }
+        else if (strcmp(dimension_opt->answer, "angle") == 0) {
+            dim_to_import = pdal::Dimension::Id::ScanAngleRank;
+        }
+        else if (strcmp(dimension_opt->answer, "class") == 0) {
+            dim_to_import = pdal::Dimension::Id::Classification;
+        }
+        else if (strcmp(dimension_opt->answer, "source") == 0) {
+            dim_to_import = pdal::Dimension::Id::PointSourceId;
+        }
+    }
+
+    if (point_binning.method == METHOD_SIDNMAX ||
+        point_binning.method == METHOD_SIDNMIN)
+        dim_to_import = pdal::Dimension::Id::PointSourceId;
+
+    if (dim_to_import != pdal::Dimension::Id::Z &&
+        (strcmp(method_opt->answer, "ev1") == 0 ||
+         strcmp(method_opt->answer, "ev2") == 0 ||
+         strcmp(method_opt->answer, "ev3") == 0))
+        dim_to_import = pdal::Dimension::Id::Z;
+
+    /* Set up axis and output value scaling */
+    double zscale = 1.0;
+    double iscale = 1.0;
+    double dscale = 1.0;
+    double output_scale = 1.0;
+
+    if (zscale_opt->answer)
+        zscale = atof(zscale_opt->answer);
+    if (iscale_opt->answer)
+        iscale = atof(iscale_opt->answer);
+    if (dscale_opt->answer)
+        dscale = atof(dscale_opt->answer);
+
+    if (zscale_opt->answer && dim_to_import == pdal::Dimension::Id::Z)
+        output_scale = zscale;
+    if (iscale_opt->answer && dim_to_import == pdal::Dimension::Id::Intensity)
+        output_scale = iscale;
+    if (dscale_opt->answer)
+        output_scale = dscale;
+
+    double res = 0.0;
+
+    if (res_opt->answer) {
+        /* align to resolution */
+        res = atof(res_opt->answer);
+
+        if (!G_scan_resolution(res_opt->answer, &res, region.proj))
+            G_fatal_error(_("Invalid input <%s=%s>"), res_opt->key,
+                          res_opt->answer);
+
+        if (res <= 0)
+            G_fatal_error(_("Option '%s' must be > 0.0"), res_opt->key);
+
+        region.ns_res = region.ew_res = res;
+
+        region.north = ceil(region.north / res) * res;
+        region.south = floor(region.south / res) * res;
+        region.east = ceil(region.east / res) * res;
+        region.west = floor(region.west / res) * res;
+
+        G_adjust_Cell_head(&region, 0, 0);
+    }
+    else if (extents_flag->answer) {
+        /* align to current region */
+        Rast_align_window(&region, &loc_wind);
+    }
+    if (base_rast_res_flag->answer) {
+        Rast_get_cellhd(base_raster_opt->answer, "", &input_region);
+        region.ns_res = input_region.ns_res;
+        region.ew_res = input_region.ew_res;
+        G_adjust_Cell_head(&region, 0, 0);
+    }
+
+    Rast_set_output_window(&region);
+    rows = region.rows;
+    cols = region.cols;
+
+    G_debug(2, "region.n=%f  region.s=%f  region.ns_res=%f", region.north,
+            region.south, region.ns_res);
+    G_debug(2, "region.rows=%d  [box_rows=%d]  region.cols=%d", region.rows,
+            rows, region.cols);
+
+    /* using segment library for the base raster */
+    // TODO: use segment library also for the binning removing the
+    // current memory limitations
+    // TODO: remove hardcoded memory requirements, let user supply it
+    int use_base_raster_res = 0;
+
+    /* TODO: see if the input region extent is smaller than the raster
+     * if yes, the we need to load the whole base raster if the -e
+     * flag was defined (alternatively clip the regions) */
+    if (base_rast_res_flag->answer)
+        use_base_raster_res = 1;
+    if (base_raster_opt->answer) {
+        if (use_base_raster_res) {
+            /* read raster actual extent and resolution */
+            Rast_get_cellhd(base_raster_opt->answer, "", &input_region);
+            /* TODO: make it only as small as the output is or points are */
+            Rast_set_input_window(&input_region);       /* we have split window */
+        }
+        else {
+            Rast_get_input_window(&input_region);
+        }
+        rast_segment_open(&base_segment, base_raster_opt->answer,
+                          &base_raster_data_type);
+    }
+
+    // TODO: use memory requirements supplied by user
+    // TODO: use segment library for binning
+    point_binning_allocate(&point_binning, rows, cols, rtype);
+
+    /* open output map */
+    out_fd = Rast_open_new(outmap, rtype);
+
+    /* allocate memory for a single row of output data */
+    raster_row = Rast_allocate_output_buf(rtype);
+
+    G_message(_("Reading data..."));
+
+    std::vector < pdal::Stage * >readers;
+    pdal::StageFactory factory;
+    pdal::MergeFilter merge_filter;
+    /* loop of input files */
+    for (int i = 0; i < infiles.num_items; i++) {
+        const char *infile = infiles.items[i];
+
+        std::string pdal_read_driver = factory.inferReaderDriver(infile);
+        if (pdal_read_driver.empty())
+            G_fatal_error("Cannot determine input file type of <%s>", infile);
+
+        pdal::Options las_opts;
+        pdal::Option las_opt("filename", infile);
+        las_opts.add(las_opt);
+        // stages created by factory are destroyed with the factory
+        pdal::Stage * reader = factory.createStage(pdal_read_driver);
+        if (!reader)
+            G_fatal_error
+                ("PDAL reader creation failed, a wrong format of <%s>",
+                 infile);
+        reader->setOptions(las_opts);
+        readers.push_back(reader);
+        merge_filter.setInput(*reader);
+    }
+
+    // we need to keep pointer to the last stage
+    // merge filter puts the n readers into one stage,
+    // so we don't have to worry about the list of stages later
+    pdal::Stage * last_stage = &merge_filter;
+    pdal::ReprojectionFilter reprojection_filter;
+
+    // we reproject when requested regardless of the input projection
+    if (reproject_flag->answer) {
+        G_message(_("Reprojecting the input to the location projection"));
+        char *proj_wkt = location_projection_as_wkt(false);
+
+        pdal::Options o4;
+        // TODO: try catch for user input error
+        if (input_srs_opt->answer)
+            o4.add < std::string > ("in_srs", input_srs_opt->answer);
+        o4.add < std::string > ("out_srs", proj_wkt);
+        reprojection_filter.setOptions(o4);
+        reprojection_filter.setInput(*last_stage);
+        last_stage = &reprojection_filter;
+    }
+
+    /* Enable all filters */
+    GrassLidarFilter grass_filter;
+
+    if (base_raster_opt->answer)
+        grass_filter.set_base_raster(&base_segment, &input_region,
+                                     base_raster_data_type);
+    if (use_spatial_filter)
+        grass_filter.set_spatial_filter(xmin, xmax, ymin, ymax);
+    if (use_zrange)
+        grass_filter.set_zrange_filter(zrange_min, zrange_max);
+    if (use_irange)
+        grass_filter.set_irange_filter(irange_min, irange_max);
+    if (use_drange)
+        grass_filter.set_drange_filter(drange_min, drange_max);
+    if (use_return_filter)
+        grass_filter.set_return_filter(return_filter_struct);
+    if (use_class_filter)
+        grass_filter.set_class_filter(class_filter);
+    grass_filter.set_z_scale(zscale);   // Default is 1 == no scale
+    grass_filter.set_intensity_scale(iscale);
+    grass_filter.set_d_scale(dscale);
+    grass_filter.setInput(*last_stage);
+
+    GrassRasterWriter binning_writer;
+
+    binning_writer.set_output_scale(output_scale);
+    binning_writer.setInput(grass_filter);
+    //stream_filter.setInput(*last_stage);
+    // there is no difference between 1 and 10k points in memory
+    // consumption, so using 10k in case it is faster for some cases
+    pdal::point_count_t point_table_capacity = 10000;
+    pdal::FixedPointTable point_table(point_table_capacity);
+    binning_writer.prepare(point_table);
+
+    // getting projection is possible only after prepare
+    if (over_flag->answer) {
+        G_important_message(_("Overriding projection check and assuming"
+                              " that the projection of input matches"
+                              " the location projection"));
+    }
+    else if (!reproject_flag->answer) {
+        pdal::SpatialReference spatial_reference =
+            merge_filter.getSpatialReference();
+        if (spatial_reference.empty())
+            G_fatal_error(_("The input dataset has undefined projection"));
+        std::string dataset_wkt = spatial_reference.getWKT();
+        bool proj_match = is_wkt_projection_same_as_loc(dataset_wkt.c_str());
+
+        if (!proj_match)
+            wkt_projection_mismatch_report(dataset_wkt.c_str());
+    }
+
+    G_important_message(_("Running PDAL algorithms..."));
+
+    // get the layout to see the dimensions
+    pdal::PointLayoutPtr point_layout = point_table.layout();
+
+    // TODO: test also z
+    // TODO: the falses for filters should be perhaps fatal error
+    // (bad input) or warning if filter was requested by the user
+
+    // update layers we are writing based on what is in the data
+    // update usage of our filters as well
+    if (point_layout->hasDim(pdal::Dimension::Id::ReturnNumber) &&
+        point_layout->hasDim(pdal::Dimension::Id::NumberOfReturns)) {
+        use_return_filter = true;
+    }
+    else {
+        use_return_filter = false;
+    }
+
+    if (point_layout->hasDim(pdal::Dimension::Id::Classification))
+        use_class_filter = true;
+    else
+        use_class_filter = false;
+
+    G_message(_("Scanning points..."));
+
+    if (user_dimension_opt->answer) {
+        dim_to_import = point_layout->findDim(user_dimension_opt->answer);
+        if (dim_to_import == pdal::Dimension::Id::Unknown)
+            G_fatal_error(_("Cannot identify the requested dimension. "
+                            "Check dimension name spelling."));
+        if (!(strcmp(dimension_opt->answer, "z") == 0))
+            G_warning(_("Both dimension and user dimension parameters are specified. "
+                       "Using '%s' as the dimension to import."),
+                      user_dimension_opt->answer);
+    }
+
+    // this is just for sure, we tested the individual dimensions before
+    // TODO: should we test Z explicitly as well?
+    if (!point_layout->hasDim(dim_to_import))
+        G_fatal_error(_("Dataset doesn't have requested dimension '%s'"
+                        " (possibly a programming error)"),
+                      pdal::Dimension::name(dim_to_import).c_str());
+
+    // TODO: add percentage printing to one of the filters
+    binning_writer.set_binning(&region, &point_binning, &bin_index_nodes,
+                               rtype, cols);
+    binning_writer.dim_to_import(dim_to_import);
+    if (base_raster_opt->answer)
+        binning_writer.set_base_raster(&base_segment, &input_region,
+                                       base_raster_data_type);
+    grass_filter.dim_to_import(dim_to_import);
+
+    // run the actual processing
+    binning_writer.execute(point_table);
+
+    /* calc stats and output */
+    G_message(_("Writing output raster map..."));
+    for (int row = 0; row < rows; row++) {
+        /* assemble final values into a row */
+        write_values(&point_binning, &bin_index_nodes, raster_row, row,
+                     cols, rtype);
+        G_percent(row, rows, 10);
+
+        /* write out line of raster data */
+        Rast_put_row(out_fd, raster_row, rtype);
+    }
+    /* free memory */
+    point_binning_free(&point_binning, &bin_index_nodes);
+    if (base_raster_opt->answer)
+        Segment_close(&base_segment);
+
+    G_percent(1, 1, 1);         /* flush */
+    G_free(raster_row);
+
+    G_message(_(GPOINT_COUNT_FORMAT " points found in input file(s)"),
+              grass_filter.num_processed());
+
+    /* close raster file & write history */
+    Rast_close(out_fd);
+
+    sprintf(title, "Raw X,Y,Z data binned into a raster grid by cell %s",
+            method_opt->answer);
+    Rast_put_cell_title(outmap, title);
+
+    Rast_short_history(outmap, "raster", &history);
+    Rast_command_history(&history);
+
+    // Hist fields are limited to 4096
+    char file_list[4096];
+
+    if (file_list_opt->answer)
+        G_snprintf(file_list, sizeof(file_list), "%s", file_list_opt->answer);
+    else
+        G_snprintf(file_list, sizeof(file_list), "%s", input_opt->answer);
+
+    Rast_set_history(&history, HIST_DATSRC_1, file_list);
+    Rast_write_history(outmap, &history);
+
+    /* set computation region to the new raster map */
+    /* TODO: should be in the done message */
+    if (set_region_flag->answer)
+        G_put_window(&region);
+
+    if (infiles.num_items > 1) {
+        sprintf(buff, _("Raster map <%s> created."
+                        " " GPOINT_COUNT_FORMAT
+                        " points from %d files found in region."), outmap,
+                grass_filter.num_passed(), infiles.num_items);
+    }
+    else {
+        sprintf(buff, _("Raster map <%s> created."
+                        " " GPOINT_COUNT_FORMAT " points found in region."),
+                outmap, grass_filter.num_passed());
+    }
+
+    G_done_msg("%s", buff);
+    G_message("Filtered spatially " GPOINT_COUNT_FORMAT " points.",
+              grass_filter.num_spatially_filtered());
+    G_message("Filtered z range " GPOINT_COUNT_FORMAT " points.",
+              grass_filter.num_zrange_filtered());
+    G_message("Filtered i range " GPOINT_COUNT_FORMAT " points.",
+              grass_filter.num_irange_filtered());
+    G_message("Filtered d range " GPOINT_COUNT_FORMAT " points.",
+              grass_filter.num_drange_filtered());
+    G_message("Filtered class " GPOINT_COUNT_FORMAT " points.",
+              grass_filter.num_class_filtered());
+    G_message("Filtered return " GPOINT_COUNT_FORMAT " points.",
+              grass_filter.num_return_filtered());
+
+    G_message("Processed into raster " GPOINT_COUNT_FORMAT " points.",
+              binning_writer.n_processed);
+
+    G_debug(1, "Processed " GPOINT_COUNT_FORMAT " points.",
+            grass_filter.num_processed());
+
+    string_list_free(&infiles);
+
+    exit(EXIT_SUCCESS);
+}

raster/r.in.pdal/point_binning.c

@@ -0,0 +1,446 @@
+/*
+ * r.in.pdal point binning logic
+ *
+ * Copyright 2011-2015, 2020 by Markus Metz, and The GRASS Development Team
+ * Authors:
+ *  Markus Metz (r.in.lidar)
+ *  Vaclav Petras (move code to separate functions)
+ *  Maris Nartiss (refactoring for r.in.pdal)
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+#include <stddef.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/raster.h>
+
+#include "point_binning.h"
+#include "bin_update.h"
+#include "bin_write.h"
+
+
+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;
+}
+
+
+void point_binning_set(struct PointBinning *point_binning, char *method,
+                       char *percentile, char *trim)
+{
+    /* figure out what maps we need in memory */
+    /*  n               n
+       min              min
+       max              max
+       range            min max         max - min
+       sum              sum c
+       mean             sum n           sum/n
+       stddev           mean m2 n       sqrt((sumsq - sum*sum/n)/n)
+       variance         mean m2 n       (sumsq - sum*sum/n)/n
+       coeff_var        mean m2 n       sqrt((sumsq - sum*sum/n)/n) / (sum/n)
+       median           n               array index to linked list
+       mode             n               array index to linked list
+       percentile       n               array index to linked list
+       skewness         n               array index to linked list
+       trimmean         n               array index to linked list
+       sidnmax          n               array index to linked list
+       sidnmin          n               array index to linked list
+       ev1, ev2, ev3    n               array index to linked list
+     */
+    point_binning->method = METHOD_NONE;
+    point_binning->bin_n = FALSE;
+    point_binning->bin_min = FALSE;
+    point_binning->bin_max = FALSE;
+    point_binning->bin_sum = FALSE;
+    point_binning->bin_m2 = FALSE;
+    point_binning->bin_z_index = FALSE;
+    point_binning->bin_cnt_index = FALSE;
+    point_binning->bin_eigenvalues = FALSE;
+    point_binning->bin_coordinates = FALSE;
+
+    point_binning->n_array = NULL;
+    point_binning->min_array = NULL;
+    point_binning->max_array = NULL;
+    point_binning->mean_array = NULL;
+    point_binning->sum_array = NULL;
+    point_binning->c_array = NULL;
+    point_binning->m2_array = NULL;
+    point_binning->index_array = NULL;
+    point_binning->x_array = NULL;
+    point_binning->y_array = NULL;
+
+    if (strcmp(method, "n") == 0) {
+        point_binning->method = METHOD_N;
+        point_binning->bin_n = TRUE;
+    }
+    if (strcmp(method, "min") == 0) {
+        point_binning->method = METHOD_MIN;
+        point_binning->bin_min = TRUE;
+    }
+    if (strcmp(method, "max") == 0) {
+        point_binning->method = METHOD_MAX;
+        point_binning->bin_max = TRUE;
+    }
+    if (strcmp(method, "range") == 0) {
+        point_binning->method = METHOD_RANGE;
+        point_binning->bin_min = TRUE;
+        point_binning->bin_max = TRUE;
+    }
+    if (strcmp(method, "sum") == 0) {
+        point_binning->method = METHOD_SUM;
+        point_binning->bin_sum = TRUE;
+    }
+    if (strcmp(method, "mean") == 0) {
+        point_binning->method = METHOD_MEAN;
+        point_binning->bin_sum = TRUE;
+        point_binning->bin_n = TRUE;
+    }
+    if (strcmp(method, "stddev") == 0) {
+        point_binning->method = METHOD_STDDEV;
+        point_binning->bin_m2 = TRUE;
+    }
+    if (strcmp(method, "variance") == 0) {
+        point_binning->method = METHOD_VARIANCE;
+        point_binning->bin_m2 = TRUE;
+    }
+    if (strcmp(method, "coeff_var") == 0) {
+        point_binning->method = METHOD_COEFF_VAR;
+        point_binning->bin_m2 = TRUE;
+    }
+    if (strcmp(method, "median") == 0) {
+        point_binning->method = METHOD_MEDIAN;
+        point_binning->bin_z_index = TRUE;
+    }
+    if (strcmp(method, "mode") == 0) {
+        point_binning->method = METHOD_MODE;
+        point_binning->bin_cnt_index = TRUE;
+    }
+    if (strcmp(method, "percentile") == 0) {
+        if (percentile != NULL)
+            point_binning->pth = atoi(percentile);
+        else
+            G_fatal_error(_("Unable to calculate percentile without the pth option specified!"));
+        point_binning->method = METHOD_PERCENTILE;
+        point_binning->bin_z_index = TRUE;
+    }
+    if (strcmp(method, "skewness") == 0) {
+        point_binning->method = METHOD_SKEWNESS;
+        point_binning->bin_z_index = TRUE;
+    }
+    if (strcmp(method, "trimmean") == 0) {
+        if (trim != NULL)
+            point_binning->trim = atof(trim) / 100.0;
+        else
+            G_fatal_error(_("Unable to calculate trimmed mean without the trim option specified!"));
+        point_binning->method = METHOD_TRIMMEAN;
+        point_binning->bin_z_index = TRUE;
+    }
+    if (strcmp(method, "sidnmax") == 0) {
+        point_binning->method = METHOD_SIDNMAX;
+        point_binning->bin_cnt_index = TRUE;
+    }
+    if (strcmp(method, "sidnmin") == 0) {
+        point_binning->method = METHOD_SIDNMIN;
+        point_binning->bin_cnt_index = TRUE;
+    }
+    if (strcmp(method, "ev1") == 0) {
+        point_binning->method = METHOD_EV1;
+        point_binning->bin_eigenvalues = TRUE;
+    }
+    if (strcmp(method, "ev2") == 0) {
+        point_binning->method = METHOD_EV2;
+        point_binning->bin_eigenvalues = TRUE;
+    }
+    if (strcmp(method, "ev3") == 0) {
+        point_binning->method = METHOD_EV3;
+        point_binning->bin_eigenvalues = TRUE;
+    }
+}
+
+
+void point_binning_allocate(struct PointBinning *point_binning, int rows,
+                            int cols, RASTER_MAP_TYPE rtype)
+{
+    if (point_binning->bin_n) {
+        G_debug(2, "allocating n_array");
+        point_binning->n_array =
+            G_calloc((size_t)rows * (cols + 1), Rast_cell_size(CELL_TYPE));
+        blank_array(point_binning->n_array, rows, cols, CELL_TYPE, 0);
+    }
+    if (point_binning->bin_min) {
+        G_debug(2, "allocating min_array");
+        point_binning->min_array =
+            G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
+        blank_array(point_binning->min_array, rows, cols, rtype, -1);   /* fill with NULLs */
+    }
+    if (point_binning->bin_max) {
+        G_debug(2, "allocating max_array");
+        point_binning->max_array =
+            G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
+        blank_array(point_binning->max_array, rows, cols, rtype, -1);   /* fill with NULLs */
+    }
+    if (point_binning->bin_sum) {
+        G_debug(2, "allocating sum_array");
+        point_binning->sum_array =
+            G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
+        blank_array(point_binning->sum_array, rows, cols, rtype, 0);
+        point_binning->c_array =
+            G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
+        blank_array(point_binning->c_array, rows, cols, rtype, 0);
+    }
+    if (point_binning->bin_m2) {
+        G_debug(2, "allocating m2_array");
+        point_binning->m2_array =
+            G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
+        blank_array(point_binning->m2_array, rows, cols, rtype, 0);
+        point_binning->mean_array =
+            G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
+        blank_array(point_binning->mean_array, rows, cols, rtype, -1);
+        point_binning->n_array =
+            G_calloc((size_t)rows * (cols + 1), Rast_cell_size(CELL_TYPE));
+        blank_array(point_binning->n_array, rows, cols, CELL_TYPE, 0);
+    }
+    if (point_binning->bin_z_index ||
+        point_binning->bin_cnt_index || point_binning->bin_eigenvalues) {
+        G_debug(2, "allocating index_array");
+        point_binning->index_array =
+            G_calloc((size_t)rows * (cols + 1), Rast_cell_size(CELL_TYPE));
+        blank_array(point_binning->index_array, rows, cols, CELL_TYPE, -1);     /* fill with NULLs */
+    }
+}
+
+
+void point_binning_free(struct PointBinning *point_binning,
+                        struct BinIndex *bin_index_nodes)
+{
+    if (point_binning->bin_n)
+        G_free(point_binning->n_array);
+    if (point_binning->bin_min)
+        G_free(point_binning->min_array);
+    if (point_binning->bin_max)
+        G_free(point_binning->max_array);
+    if (point_binning->bin_sum) {
+        G_free(point_binning->sum_array);
+        G_free(point_binning->c_array);
+    }
+    if (point_binning->bin_m2) {
+        G_free(point_binning->m2_array);
+        G_free(point_binning->mean_array);
+        G_free(point_binning->n_array);
+    }
+    if (point_binning->bin_z_index ||
+        point_binning->bin_cnt_index || point_binning->bin_eigenvalues) {
+        G_free(point_binning->index_array);
+        G_free(bin_index_nodes->nodes);
+        bin_index_nodes->num_nodes = 0;
+        bin_index_nodes->max_nodes = 0;
+        bin_index_nodes->nodes = NULL;
+    }
+}
+
+
+void write_values(struct PointBinning *point_binning,
+                  struct BinIndex *bin_index_nodes, void *raster_row, int row,
+                  int cols, RASTER_MAP_TYPE rtype)
+{
+    void *ptr = NULL;
+    int col;
+
+    switch (point_binning->method) {
+    case METHOD_N:             /* n is a straight copy */
+        Rast_raster_cpy(raster_row,
+                        ((char *)point_binning->n_array) +
+                        ((size_t)row * cols * Rast_cell_size(CELL_TYPE)),
+                        cols, CELL_TYPE);
+        break;
+
+    case METHOD_MIN:
+        Rast_raster_cpy(raster_row,
+                        ((char *)point_binning->min_array) +
+                        ((size_t)row * cols * Rast_cell_size(rtype)), cols,
+                        rtype);
+        break;
+
+    case METHOD_MAX:
+        Rast_raster_cpy(raster_row,
+                        ((char *)point_binning->max_array) +
+                        ((size_t)row * cols * Rast_cell_size(rtype)), cols,
+                        rtype);
+        break;
+
+    case METHOD_SUM:
+        write_sum(raster_row, point_binning->sum_array,
+                  point_binning->c_array, row, cols, rtype);
+        break;
+
+    case METHOD_RANGE:         /* (max-min) */
+        ptr = raster_row;
+        for (col = 0; col < cols; col++) {
+            size_t offset =
+                ((size_t)row * cols + col) * Rast_cell_size(rtype);
+            double min =
+                Rast_get_d_value(((char *)point_binning->min_array) + offset,
+                                 rtype);
+            double max =
+                Rast_get_d_value(((char *)point_binning->max_array) + offset,
+                                 rtype);
+
+            Rast_set_d_value(ptr, max - min, rtype);
+            ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
+        }
+        break;
+
+    case METHOD_MEAN:          /* (sum / n) */
+        ptr = raster_row;
+        for (col = 0; col < cols; col++) {
+            size_t n_offset =
+                ((size_t)row * cols + col) * Rast_cell_size(CELL_TYPE);
+            int n =
+                Rast_get_c_value(((char *)point_binning->n_array) + n_offset,
+                                 CELL_TYPE);
+            double sum =
+                get_sum(point_binning->sum_array, point_binning->c_array,
+                        row, cols, col, rtype);
+
+            if (n == 0)
+                Rast_set_null_value(ptr, 1, rtype);
+            else
+                Rast_set_d_value(ptr, (sum / n), rtype);
+
+            ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
+        }
+        break;
+
+    case METHOD_STDDEV:        /*  sqrt(variance)        */
+    case METHOD_VARIANCE:      /*  (sumsq - sum*sum/n)/n */
+    case METHOD_COEFF_VAR:     /*  100 * stdev / mean    */
+        write_variance(raster_row, point_binning->n_array,
+                       point_binning->mean_array, point_binning->m2_array,
+                       row, cols, rtype, point_binning->method);
+        break;
+
+    case METHOD_MEDIAN:        /* median, if only one point in cell we will use that */
+        write_median(bin_index_nodes, raster_row, point_binning->index_array,
+                     row, cols, rtype);
+        break;
+
+    case METHOD_MODE:
+        write_mode(bin_index_nodes, raster_row, point_binning->index_array,
+                   row, cols);
+        break;
+
+    case METHOD_PERCENTILE:    /* rank = (pth*(n+1))/100; interpolate linearly */
+        write_percentile(bin_index_nodes, raster_row,
+                         point_binning->index_array, row, cols, rtype,
+                         point_binning->pth);
+        break;
+
+    case METHOD_SKEWNESS:      /* skewness = sum(xi-mean)^3/(N-1)*s^3 */
+        write_skewness(bin_index_nodes, raster_row,
+                       point_binning->index_array, row, cols, rtype);
+        break;
+
+    case METHOD_TRIMMEAN:
+        write_trimmean(bin_index_nodes, raster_row,
+                       point_binning->index_array, row, cols, rtype,
+                       point_binning->trim);
+        break;
+
+    case METHOD_SIDNMAX:
+        write_sidn(bin_index_nodes, raster_row, point_binning->index_array,
+                   row, cols, 0);
+        break;
+
+    case METHOD_SIDNMIN:
+        write_sidn(bin_index_nodes, raster_row, point_binning->index_array,
+                   row, cols, 1);
+        break;
+
+    case METHOD_EV1:
+    case METHOD_EV2:
+    case METHOD_EV3:
+        write_ev(bin_index_nodes, raster_row, point_binning->index_array,
+                 row, cols, rtype, point_binning->method);
+        break;
+
+    default:
+        G_debug(2, "No method selected");
+    }
+}
+
+
+void update_value(struct PointBinning *point_binning,
+                  struct BinIndex *bin_index_nodes, int cols, int arr_row,
+                  int arr_col, RASTER_MAP_TYPE rtype, double x, double y,
+                  double z)
+{
+    if (point_binning->bin_n)
+        update_n(point_binning->n_array, cols, arr_row, arr_col);
+    if (point_binning->bin_min)
+        update_min(point_binning->min_array, cols, arr_row, arr_col, rtype,
+                   z);
+    if (point_binning->bin_max)
+        update_max(point_binning->max_array, cols, arr_row, arr_col, rtype,
+                   z);
+    if (point_binning->bin_sum)
+        update_sum(point_binning->sum_array, point_binning->c_array, cols,
+                   arr_row, arr_col, rtype, z);
+    if (point_binning->bin_m2)
+        update_m2(point_binning->n_array, point_binning->mean_array,
+                  point_binning->m2_array, cols, arr_row, arr_col, rtype, z);
+    if (point_binning->bin_z_index)
+        update_bin_z_index(bin_index_nodes, point_binning->index_array, cols,
+                           arr_row, arr_col, z);
+    if (point_binning->bin_cnt_index)
+        update_bin_cnt_index(bin_index_nodes, point_binning->index_array,
+                             cols, arr_row, arr_col, z);
+    if (point_binning->bin_eigenvalues)
+        update_bin_com_index(bin_index_nodes, point_binning->index_array,
+                             cols, arr_row, arr_col, x, y, z);
+}

raster/r.in.pdal/point_binning.h

@@ -0,0 +1,116 @@
+/*
+ * r.in.pdal point binning logic
+ *
+ * Copyright 2011-2015, 2020 by Markus Metz, and The GRASS Development Team
+ * Authors:
+ *  Markus Metz (r.in.lidar)
+ *  Vaclav Petras (move code to separate functions)
+ *  Maris Nartiss (refactoring for r.in.pdal)
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+#ifndef __POINT_BINNING_H__
+#define __POINT_BINNING_H__
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+
+/* Point binning methods: */
+#define METHOD_NONE        0
+#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_MODE       11
+#define METHOD_PERCENTILE 12
+#define METHOD_SKEWNESS   13
+#define METHOD_TRIMMEAN   14
+#define METHOD_SIDNMAX    15
+#define METHOD_SIDNMIN    16
+#define METHOD_EV1        17
+#define METHOD_EV2        18
+#define METHOD_EV3        19
+
+
+struct z_node
+{
+    int next;
+    double z;
+};
+
+struct cnt_node
+{
+    int next;
+    CELL value;
+    int count;
+};
+
+struct com_node
+{
+    int n;
+    double *meanx;
+    double *meany;
+    double *comoment;
+};
+
+struct BinIndex
+{
+    int num_nodes;
+    int max_nodes;
+    void *nodes;
+};
+
+struct PointBinning
+{
+    int method;
+
+    int bin_n;
+    int bin_min;
+    int bin_max;
+    int bin_sum;
+    int bin_m2;
+    int bin_z_index;
+    int bin_cnt_index;
+    int bin_eigenvalues;
+    int bin_coordinates;
+
+    void *n_array;
+    void *min_array;
+    void *max_array;
+    void *sum_array;
+    void *c_array;
+    void *mean_array;
+    void *m2_array;
+    void *index_array;
+    void *x_array;
+    void *y_array;
+
+    int pth;
+    double trim;
+};
+
+void *get_cell_ptr(void *, int, int, int, RASTER_MAP_TYPE);
+int blank_array(void *, int, int, RASTER_MAP_TYPE, int);
+
+void point_binning_set(struct PointBinning *, char *, char *, char *);
+void point_binning_allocate(struct PointBinning *, int, int, RASTER_MAP_TYPE);
+void point_binning_free(struct PointBinning *, struct BinIndex *);
+
+
+void write_values(struct PointBinning *,
+                  struct BinIndex *, void *, int, int, RASTER_MAP_TYPE);
+void update_value(struct PointBinning *,
+                  struct BinIndex *, int, int,
+                  int, RASTER_MAP_TYPE, double, double, double);
+
+
+#endif /* __POINT_BINNING_H__ */

raster/r.in.pdal/projection.c

@@ -0,0 +1,214 @@
+/*
+ * projection checking
+ *
+ * Copyright 2011-2015 by Markus Metz, and The GRASS Development Team
+ * Authors:
+ *  Markus Metz (v.in.lidar)
+ *  Vaclav Petras (move code to standalone functions)
+ *
+ * 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>
+
+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(_("Overriding 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"));
+        }
+    }
+}
+
+
+/* Does the projection of the current location match the dataset? */
+int is_wkt_projection_same_as_loc(const char *wkt)
+{
+    struct Key_Value *loc_proj_info = NULL, *loc_proj_units = NULL;
+    struct Key_Value *proj_info = NULL, *proj_units = NULL;
+    struct Cell_head cellhd;
+    struct Cell_head loc_wind;
+
+    G_get_default_window(&loc_wind);
+
+    /* Projection only required for checking so convert non-interactively */
+    if (GPJ_wkt_to_grass(&cellhd, &proj_info, &proj_units, wkt, 0) < 0)
+        G_warning(_("Unable to convert input map projection information to "
+                    "GRASS format for checking"));
+
+    /* fetch LOCATION PROJ info */
+    if (loc_wind.proj != PROJECTION_XY) {
+        loc_proj_info = G_get_projinfo();
+        loc_proj_units = G_get_projunits();
+    }
+
+    if (loc_wind.proj != cellhd.proj) {
+        return FALSE;
+    }
+    else if (G_compare_projections(loc_proj_info, loc_proj_units,
+                                   proj_info, proj_units) != 1) {
+        return FALSE;
+    }
+    else {
+        return TRUE;
+    }
+}
+
+/* Does the projection of the current location match the dataset? */
+void wkt_projection_mismatch_report(const char *wkt)
+{
+    struct Key_Value *loc_proj_info = NULL, *loc_proj_units = NULL;
+    struct Key_Value *proj_info = NULL, *proj_units = NULL;
+    struct Cell_head cellhd;
+    struct Cell_head loc_wind;
+
+    G_get_default_window(&loc_wind);
+
+    /* Projection only required for checking so convert non-interactively */
+    if (GPJ_wkt_to_grass(&cellhd, &proj_info, &proj_units, wkt, 0) < 0)
+        G_warning(_("Unable to convert input map projection information to "
+                    "GRASS format for checking"));
+
+    /* fetch LOCATION PROJ info */
+    if (loc_wind.proj != PROJECTION_XY) {
+        loc_proj_info = G_get_projinfo();
+        loc_proj_units = G_get_projunits();
+    }
+    int err = G_compare_projections(loc_proj_info, loc_proj_units,
+                                    proj_info, proj_units);
+
+    projection_mismatch_report(cellhd, loc_wind, loc_proj_info,
+                               loc_proj_units, proj_info, proj_units, err);
+}
+
+/* caller should free the returned string */
+char *location_projection_as_wkt(int prettify)
+{
+    struct Key_Value *proj_info = G_get_projinfo();
+    struct Key_Value *proj_units = G_get_projunits();
+    char *proj_wkt = GPJ_grass_to_wkt(proj_info, proj_units, FALSE, prettify);
+
+    G_free_key_value(proj_info);
+    G_free_key_value(proj_units);
+    return proj_wkt;
+}

raster/r.in.pdal/projection.h

@@ -0,0 +1,36 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.in.pdal
+ *               adapted from v.in.lidar
+ * AUTHOR(S):    Vaclav Petras
+ * PURPOSE:      projection related functions
+ * COPYRIGHT:    (C) 2015 by the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *               License (>=v2). Read the COPYING file that comes with GRASS
+ *               for details.
+ *
+ *****************************************************************************/
+
+
+#ifndef PROJECTION_CHECKS_H
+#define PROJECTION_CHECKS_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 return_value, int verbose);
+
+int is_wkt_projection_same_as_loc(const char *wkt);
+void wkt_projection_mismatch_report(const char *wkt);
+char *location_projection_as_wkt(int prettify);
+
+#endif /* PROJECTION_CHECKS_H */

raster/r.in.pdal/r.in.pdal.html

@@ -0,0 +1,675 @@
+<h2>DESCRIPTION</h2>
+
+The <em>r.in.pdal</em> module loads PDAL library supported point clouds
+(with emphasis on LiDAR LAS files) 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 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.pdal</em> is designed for processing massive point cloud
+datasets, for example raw LiDAR or sidescan sonar swath data.
+
+<h3>Binning</h3>
+
+The main difference between <em>r.in.pdal</em> and
+<em><a href="v.in.pdal.html">v.in.pdal</a></em> is that
+<em>r.in.pdal</em> creates a raster instead of just importing the
+points into GRASS GIS. However, <em>r.in.pdal</em> does not merely
+rasterizes the points from the point cloud. <em>r.in.pdal</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 binning of points with X and Y coordinates starts with the overlay
+of a grid of bins over the points.
+<p>
+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.pdal</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.
+
+<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>
+
+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.</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).
+It is calculated by Welford algorithm.</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>100 * sqrt(variance) / mean</tt>.</dd>
+<dt>median</dt>
+<dd>This computes the median of point values for each cell</dd>
+<dt>mode</dt>
+<dd>This computes the mode (the most common value) of point values for each cell.
+If there are more than one mode for a cell, an arbitrary one will be selected.</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>
+<dt>sidnmax</dt>
+<dd>Maximum number of points in cell per source ID.
+Points are counted for each source separately and largest count
+is used as the final cell value. Corresponds to majority density (<i>Dmax</i>)
+as defined by Smeeckaert et al., 2013.</dd>
+<dt>sidnmin</dt>
+<dd>Minimum number of points in cell per source ID.
+Points are counted for each source separately and smallest count
+is used as the final cell value. Corresponds to minority density (<i>Dmin</i>)
+as defined by Smeeckaert et al., 2013.</dd>
+<dt>ev1, ev2, ev3</dt>
+<dd>1st, 2nd and 3rd eigenvalue of point x, y and z coordinates within
+a cell (Mallet et al. 2011).</dd>
+</dl>
+
+Note that different statistics have different memory requirements
+(see below for details).
+
+<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.
+
+<h3>Filtering and selection</h3>
+
+Points are always filtered by spatial extent (based on current computational
+region) unless <em>-e</em> flag is set. Afterwards import value scaling
+is performed, followed by filtering by having value in range.
+
+<p>
+Points falling outside the current computational region will be skipped.
+This includes points falling <em>exactly</em> on the southern or eastern
+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>
+Filtering can be performed also on intensity values with <em>irange</em>
+parameter. If imported dimension is neither Z nor intensity, filtering
+can be performed with <em>drange</em> parameter. Filtering by Z values
+and intensity can be performed independently on imported dimension.
+If imported dimension is not Z nor intensity, filtering can be performed
+with all three parameters simultaneously (<em>zrange</em>, <em>irange</em>
+and <em>drange</em>) or any of their combination.
+
+<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>
+The user can use a combination of <em>r.in.pdal</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.
+
+<p>
+LiDAR points can contain a lot of attributes besides their coordinates.
+By default <em>r.in.pdal</em> will use point Z coordinate as the output
+raster value source. The <em>dimension</em> option allows to use different
+point attributes as an output instead of Z value. In case if other than Z
+value is choosen, output map type will be set to CELL. Explanation
+of various attributes can be found in LAS Specification chapter "Point
+Data Records". Be ware - not all attributes might be present and also
+their value ranges and meaning might vary between LAS file producers.
+
+<p>
+PDAL supports creation of custom dimensions (point attributes) in LAS files.
+Use <em>user_dimension</em> option to provide name of any dimension for
+import. This option supersedes one specified by <em>dimension</em> option.
+
+<h3>Reduction to a base raster</h3>
+
+<p>
+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.pdal</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.pdal</em>'s <b>-g</b> flag to find the input
+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>
+The <b>-g</b> flag prints extent of input data set(s) in a shell style
+suitable as command line parameters for <em>g.region</em>. Note, extent
+is determined from data set(s) metadata and without scanning whole
+dataset and thus selectros and filters are not applied.
+
+<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.
+
+<h3>Transformation</h3>
+<p>
+Point Z axis values can be altered by scaling parameter set by <em>zscale</em>.
+Each point X axis value is multiplied by provided value before any
+filtering or further transformation takes place. In the same way act
+<em>iscale</em> and <em>dscale</em> parametres – values are multiplied
+before range filters <em>irange</em> or <em>drnage</em> are applied (if present).
+<p>
+Output value scaling can be performed independently of filtering by providing
+<em>dscale</em> parameter. This scaling parameter is applied only to
+output value independently of its dimension. Thus if <em>zrange</em> and
+<em>dscale</em> are set for Z output dimension, filtering by <em>zrange</em> 
+will be done with unscaled values but output will be scaled by <em>dscale</em>.
+<em>dscale</em> takes precedence over <em>zscale</em> and <em>iscale</em>
+if more than one parameter is provided.
+
+<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 consumption</h3>
+
+<p>
+While the <b>input</b> file can be arbitrarily large, <em>r.in.pdal</em>
+will use a large amount of system memory (RAM) for large raster regions
+(&gt; 10000x10000 pixels).
+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.
+For <b>methods</b>=<em>n, mode, sidnmin, sidnmax</em>, the <tt>CELL</tt>
+format is used automatically.
+
+<p>
+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>
+The memory usage for regular statistics mentioned above is based solely
+on region (raster) size.
+However, the aggregate functions <em>median, mode, 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>Trim option</h3>
+<p>
+Trim option value is used only when calculating trimmed mean values.
+Attempt to use it with other statistical methods will result in an error.
+
+<h3>Print flag</h3>
+<p>
+The <b>-p</b> flag will use PDAL library built-in LAS file metadata
+printing function. It will also include a list of supported dimension
+names. Any of listed dimensions can be passed to the <b>user_dimension</b>
+option. Dimensions are not checked for presence of data for any of points
+in the file and thus importing a dimension without actual data will
+result in an empty map.
+
+<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.pdal -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>
+
+<p>
+For the output raster map, a <b>suitable resolution</b> can be found by
+dividing the number of input points by the area covered (this requires
+an iterative approach as outlined here):
+
+<!-- points.las is from California -->
+<div class="code"><pre>
+# print LAS metadata (Number of Points)
+r.in.pdal -p input=points.las
+# Point count: 1287775
+
+# scan for LAS points cloud extent
+r.in.pdal -g input=points.las output=dummy -o
+# n=2193507.740000 s=2190053.450000 e=6070237.920000 w=6066629.860000 b=-3.600000 t=906.000000
+
+# set computation region to this extent
+g.region n=2193507.740000 s=2190053.450000 e=6070237.920000 w=6066629.860000 -p
+
+# print resulting extent
+g.region -p
+#  rows:       3454
+#  cols:       3608
+
+# points_per_cell = n_points / (rows * cols)
+# Here: 1287775 / (3454 * 3608) = 0.1033359 LiDAR points/raster cell
+# As this is too low, we need to select a lower raster resolution
+
+g.region res=5 -ap
+#  rows:       692
+#  cols:       723
+#  Now: 1287775 / (692 * 723) = 2.573923 LiDAR points/raster cell
+
+# import as mean
+r.in.pdal input=points.las output=lidar_dem_mean method=mean -o
+
+# import as max
+r.in.pdal input=points.las output=lidar_dem_max method=max -o
+
+# import as p'th percentile of the values
+r.in.pdal input=points.las output=lidar_dem_percentile_95 \
+           method=percentile pth=95 -o
+</pre></div>
+
+<center>
+<img src="r_in_lidar_dem_mean3D.jpg" alt="Mean value DEM in perspective view"><br>
+<i>Mean value DEM in perspective view, imported from LAS file</i>
+</center>
+
+<p>
+Further hints: how to calculate number of LiDAR points/square meter:
+<div class="code"><pre>
+g.region -e
+  # Metric 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>
+
+<h3>Serpent Mound dataset</h3>
+
+This example is analogous to the example used in the GRASS wiki page for
+<a href="http://grasswiki.osgeo.org/wiki/LIDAR#Import_LAS_as_raster_DEM">importing LAS as raster DEM</a>.
+<p>The sample LAS data are in the file "Serpent Mound Model LAS Data.las",
+available at
+<a href="http://www.appliedimagery.com/downloads/sampledata/Serpent%20Mound%20Model%20LAS%20Data.las">appliedimagery.com</a>:
+
+<div class="code"><pre>
+# print LAS file info
+r.in.pdal -p input="Serpent Mound Model LAS Data.las"
+
+# using v.in.lidar to create a new location
+# create location with projection information of the LAS data
+v.in.lidar -i input="Serpent Mound Model LAS Data.las" location=Serpent_Mound
+
+# quit and restart GRASS in the newly created location "Serpent_Mound"
+
+# scan the extents of the LAS data
+r.in.pdal -g input="Serpent Mound Model LAS Data.las"
+
+# set the region to the extents of the LAS data, align to resolution
+g.region n=4323641.57 s=4320942.61 w=289020.90 e=290106.02 res=1 -ap
+
+# import as raster DEM
+r.in.pdal input="Serpent Mound Model LAS Data.las" \
+           output=Serpent_Mound_Model_LAS_Data method=mean
+</pre></div>
+
+<center>
+<img src="r_in_lidar.png">
+<p><em>Figure: Elevation for the whole area of Serpent Mound dataset</em></p>
+</center>
+
+<h3>Height above ground</h3>
+
+The mean height above ground of the points can be computed for each
+raster cell (the ground elevation is given by the raster map
+<tt>elevation</tt>):
+
+<div class="code"><pre>
+g.region raster=elevation -p
+r.in.pdal input=points.las output=mean_height_above_ground base_raster=elevation method=mean
+</pre></div>
+
+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
+path. For example, a list of files in the directory /home/user/data:
+<div class="code"><pre>
+points1.laz
+points2.laz
+points3.laz
+</pre></div>
+
+would be lised in the file as:
+<div class="code"><pre>
+/home/user/data/points1.laz
+/home/user/data/points2.laz
+/home/user/data/points3.laz
+</pre></div>
+On Linux and OSX, this file can be automatically generated with the command:
+<div class="code"><pre>
+ls /home/user/data/*.laz > /home/user/data/filelist.txt
+</pre></div>
+On Windows:
+<div class="code"><pre>
+dir /b c:\users\user\data\*.laz > c:\users\user\data\filelist.txt
+</pre></div>
+
+The mean height above ground example above would then be:
+
+<div class="code"><pre>
+g.region raster=elevation -p
+r.in.pdal 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.pdal', file=file_list_name,
+                    output='mean_height_above_ground',
+                    base_raster='elevation' method='mean')
+</pre></div>
+
+
+<h2>KNOWN ISSUES</h2>
+
+<ul>
+<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.
+</ul>
+
+If you encounter any problems (or solutions!) please contact the GRASS
+Development Team.
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="g.region.html">g.region</a>,
+<a href="r.in.xyz.html">r.in.xyz</a>,
+<a href="r.mapcalc.html">r.mapcalc</a>,
+<a href="r.univar.html">r.univar</a>,
+<a href="v.in.lidar.html">v.in.lidar</a>,
+<a href="r3.in.lidar.html">r3.in.lidar</a>,
+<a href="v.vect.stats.html">v.vect.stats</a>
+<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>
+</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>
+V. Petras, A. Petrasova, J. Jeziorska, H. Mitasova (2016):
+<em>Processing UAV and lidar point clouds in GRASS GIS</em>.
+XXIII ISPRS Congress 2016 [<a href="http://www.int-arch-photogramm-remote-sens-spatial-inf-sci.net/XLI-B7/945/2016/">ISPRS Archives</a>, <a href="https://www.researchgate.net/publication/304340172_Processing_UAV_and_lidar_point_clouds_in_GRASS_GIS">ResearchGate</a>]
+<li>
+<a href="http://www.asprs.org/Committee-General/LASer-LAS-File-Format-Exchange-Activities.html">
+ASPRS LAS format</a>
+<li>
+<a href="https://pdal.io/">PDAL - Point Data Abstraction Library</a>
+</ul>
+
+<h2>AUTHORS</h2>
+
+Markus Metz<br>
+Vaclav Petras,
+<a href="http://geospatial.ncsu.edu/osgeorel/">NCSU GeoForAll Lab</a>
+(base_raster option, documentation),<br>
+Maris Nartiss, <a href="https://www.geo.lu.lv/">LU GZZF</a>
+(refactoring, additional filters, custom dimension support)
+<br>
+based on <em>r.in.xyz</em> by Hamish Bowman and Volker Wichmann<br>
+
+<!--
+<p>
+<i>Last changed: $Date$</i>
+-->

raster/r.in.pdal/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>

raster/r.in.pdal/rast_segment.c

@@ -0,0 +1,102 @@
+
+ /****************************************************************************
+ *
+ * MODULE:    r.in.pdal
+ *
+ * AUTHOR(S): Vaclav Petras
+ *
+ * PURPOSE:   This file is a wrapper for a subset of Segment functions
+ *
+ * COPYRIGHT: (C) 2015 by Vaclav Petras and 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 <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/raster.h>
+#include <grass/segment.h>
+
+#include "rast_segment.h"
+
+static void rast_segment_load(SEGMENT * segment, int rowio,
+                              RASTER_MAP_TYPE map_type)
+{
+    void *raster_row = Rast_allocate_input_buf(map_type);
+    int row;
+
+    for (row = 0; row < Rast_input_window_rows(); row++) {
+        /* TODO: free mem */
+        Rast_get_row(rowio, raster_row, row, map_type);
+        Segment_put_row(segment, raster_row, row);
+    }
+}
+
+/* TODO: close function */
+
+void rast_segment_open(SEGMENT * segment, const char *name,
+                       RASTER_MAP_TYPE * map_type)
+{
+    /* TODO: check if not passing the mapset is OK */
+    int rowio = Rast_open_old(name, "");
+
+    *map_type = Rast_get_map_type(rowio);
+    int segment_rows = 64;
+
+    /* we use long segments because this is how the values a binned */
+    int segment_cols = Rast_input_window_cols();
+    int segments_in_memory = 4;
+
+    if (Segment_open(segment, G_tempfile(), Rast_input_window_rows(),
+                     Rast_input_window_cols(), segment_rows, segment_cols,
+                     Rast_cell_size(*map_type), segments_in_memory) != 1)
+        G_fatal_error(_("Cannot create temporary file with segments of a raster map"));
+    rast_segment_load(segment, rowio, *map_type);
+    Rast_close(rowio);          /* we won't need the raster again */
+}
+
+
+/* 0 on out of region or NULL, 1 on success */
+int rast_segment_get_value_xy(SEGMENT * base_segment,
+                              struct Cell_head *input_region,
+                              RASTER_MAP_TYPE rtype, double x, double y,
+                              double *value)
+{
+    /* Rast gives double, Segment needs off_t */
+    off_t base_row = Rast_northing_to_row(y, input_region);
+    off_t base_col = Rast_easting_to_col(x, input_region);
+
+    /* skip points which are outside the base raster
+     * (null propagation) */
+    if (base_row < 0 || base_col < 0 ||
+        base_row >= input_region->rows || base_col >= input_region->cols)
+        return 0;
+    if (rtype == DCELL_TYPE) {
+        DCELL tmp;
+
+        Segment_get(base_segment, &tmp, base_row, base_col);
+        if (Rast_is_d_null_value(&tmp))
+            return 0;
+        *value = (double)tmp;
+    }
+    else if (rtype == FCELL_TYPE) {
+        FCELL tmp;
+
+        Segment_get(base_segment, &tmp, base_row, base_col);
+        if (Rast_is_f_null_value(&tmp))
+            return 0;
+        *value = (double)tmp;
+    }
+    else {
+        CELL tmp;
+
+        Segment_get(base_segment, &tmp, base_row, base_col);
+        if (Rast_is_c_null_value(&tmp))
+            return 0;
+        *value = (double)tmp;
+    }
+    return 1;
+}

raster/r.in.pdal/rast_segment.h

@@ -0,0 +1,31 @@
+
+ /****************************************************************************
+ *
+ * MODULE:    r.in.pdal
+ *
+ * AUTHOR(S): Vaclav Petras
+ *
+ * PURPOSE:   This file is a wrapper for a subset of Segment functions
+ *
+ * COPYRIGHT: (C) 2015 by Vaclav Petras and 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 GRASS_RAST_SEGMENT_H
+#define GRASS_RAST_SEGMENT_H
+
+#include <grass/raster.h>
+#include <grass/segment.h>
+
+void rast_segment_open(SEGMENT * segment, const char *name,
+                       RASTER_MAP_TYPE * map_type);
+int rast_segment_get_value_xy(SEGMENT * base_segment,
+                              struct Cell_head *input_region,
+                              RASTER_MAP_TYPE rtype, double x, double y,
+                              double *value);
+
+#endif /* GRASS_RAST_SEGMENT_H */

raster/r.in.pdal/string_list.c

@@ -0,0 +1,84 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.in.pdal
+ *
+ * AUTHOR(S):    Vaclav Petras
+ *
+ * PURPOSE:      Imports LAS LiDAR point clouds to a raster map using
+ *               aggregate statistics.
+ *
+ * COPYRIGHT:    (C) 2019 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 "string_list.h"
+
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include <stdio.h>
+#include <string.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;
+}

raster/r.in.pdal/string_list.h

@@ -0,0 +1,39 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.in.pdal
+ *
+ * AUTHOR(S):    Vaclav Petras
+ *
+ * PURPOSE:      Imports LAS LiDAR point clouds to a raster map using
+ *               aggregate statistics.
+ *
+ * COPYRIGHT:    (C) 2019 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.
+ *
+ *****************************************************************************/
+
+#ifndef __STRING_LIST_H__
+#define __STRING_LIST_H__
+
+
+/* A list which keeps multiple strings
+ *
+ * Intended for list of file names read from a file.
+ */
+
+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__ */

raster/r.in.pdal/testsuite/data/points.csv

@@ -0,0 +1,54 @@
+X,Y,Z,Intensity,Classification,PointSourceId,NumberOfReturns,ReturnNumber,CellID
+1,1,1,10,2,1,1,1,0
+2,1,1,12,2,1,2,2,0
+3,1,2,9,2,1,2,1,0
+4,1,1,11,3,1,1,1,0
+5,1,2,12,2,1,1,1,0
+2,2,2,10,2,1,2,1,0
+4,2,1,11,3,1,1,1,0
+5,2,2,10,2,1,3,3,0
+13,1,5,80,5,1,4,2,2
+13,5,2,10,2,2,1,1,2
+17,1,2,11,2,1,1,1,2
+2,7,3,15,2,1,1,1,3
+3,7,2,9,2,1,1,1,3
+4,7,2,11,2,1,1,1,3
+1,8,3,10,2,1,1,1,3
+2,8,2,11,2,1,1,1,3
+3,8,2,10,2,1,1,1,3
+4,8,2,10,2,1,1,1,3
+5,8,2,11,2,1,1,1,3
+1,9,20,50,4,1,1,1,3
+2,9,2,10,2,1,1,1,3
+3,9,3,11,2,1,1,1,3
+4,9,3,12,2,1,1,1,3
+5,9,3,11,2,1,1,1,3
+1,10,2,10,2,1,1,1,3
+2,10,2,10,2,1,1,1,3
+3,10,3,11,2,1,1,1,3
+4,10,3,12,2,1,1,1,3
+5,10,3,9,2,1,1,1,3
+2,11,2,11,2,1,1,1,3
+3,11,3,12,2,1,1,1,3
+4,11,2,10,2,1,1,1,3
+7,7,2,11,2,1,1,1,4
+7,7,15,55,5,1,3,3,4
+7,7,9,45,4,1,3,2,4
+11,7,3,11,2,2,1,1,4
+7,11,2,12,2,1,1,1,4
+11,11,2,10,2,2,1,1,4
+15,8,2,11,2,2,1,1,5
+14,9,2,9,2,1,1,1,5
+16,9,2,11,2,2,1,1,5
+15,10,2,10,2,2,1,1,5
+4,14,3,11,2,1,1,1,6
+10,13,9,67,4,2,2,1,7
+9,14,23,78,5,1,3,1,7
+8,15,19,69,5,1,3,2,7
+9,16,25,88,5,2,3,3,7
+8,17,18,65,4,1,3,3,7
+11,17,28,83,5,2,1,1,7
+16,14,3,9,2,2,1,1,8
+13,15,2,11,2,2,1,1,8
+15,15,15,51,4,2,2,2,8
+15,17,3,12,2,2,1,1,8

raster/r.in.pdal/testsuite/data/res_base_adj.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   float
+
+0 7.6666 9.25
+16.7142 9.5 0
+0.5 n 2

raster/r.in.pdal/testsuite/data/res_coeff_var_z.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   float
+
+0 30.0049 93.149
+114.7035 89.535 0
+33.3333 n 47.1404

raster/r.in.pdal/testsuite/data/res_ev1_z.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   float
+
+n 48.7647 38.3416
+15.1179 31.2019 0.6666
+2.2731 n 8.0

raster/r.in.pdal/testsuite/data/res_ev2_z.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   float
+
+n 2.3426 1.7489
+1.7000 3.3726 0.6666
+0.2923 n 5.6666

raster/r.in.pdal/testsuite/data/res_ev3_z.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   float
+
+n 1.0926 0.2023
+1.4963 3.0587 0.0
+0.0 n 0.0

raster/r.in.pdal/testsuite/data/res_filter_z_int_source.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   int
+
+1 n n
+1 1 n
+1 n 1

raster/r.in.pdal/testsuite/data/res_max_z.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   int
+
+3 28 15
+20 15 2
+2 n 5

raster/r.in.pdal/testsuite/data/res_mean_intensity.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   int
+
+11 75 20
+12 24 10
+10 n 33

raster/r.in.pdal/testsuite/data/res_mean_z.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   float
+
+3.0 20.3333 5.75
+3.2857 5.5 2.0
+1.5 n 3.0

raster/r.in.pdal/testsuite/data/res_median_z.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   float
+
+3.0 21.0 3.0
+2.0 2.5 2.0
+1.5 n 2.0

raster/r.in.pdal/testsuite/data/res_min_z.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   int
+
+3 9 2
+2 2 2
+1 n 2

raster/r.in.pdal/testsuite/data/res_mode_cellid.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   int
+
+6 7 8
+3 4 5
+0 n 2

raster/r.in.pdal/testsuite/data/res_mode_z.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   int
+
+3 9 3
+2 2 2
+1 n 2

raster/r.in.pdal/testsuite/data/res_n_all.ascii

@@ -0,0 +1,11 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+type:   int
+
+1 6 4
+21 6 4
+8 0 3

raster/r.in.pdal/testsuite/data/res_n_class_2.ascii

@@ -0,0 +1,11 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+type:   int
+
+1 0 3
+20 4 4
+6 0 2

raster/r.in.pdal/testsuite/data/res_range_z.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   int
+
+0 19 13
+18 13 0
+1 n 3

raster/r.in.pdal/testsuite/data/res_sidnmax_all.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   int
+
+1 3 4
+21 4 3
+8 0 2

raster/r.in.pdal/testsuite/data/res_sidnmin_all.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   int
+
+1 3 4
+21 2 1
+8 0 1

raster/r.in.pdal/testsuite/data/res_stddev_z.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   float
+
+0 6.101 5.356
+3.7688 4.9244 0
+0.5 n 1.4142

raster/r.in.pdal/testsuite/data/res_sum_z.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   int
+
+3 122 23
+69 33 8
+12 0 9

raster/r.in.pdal/testsuite/data/res_variance_z.ascii

@@ -0,0 +1,12 @@
+north: 18
+south:  0
+east:  18
+west:   0
+rows:   3
+cols:   3
+null:   n
+type:   float
+
+0 37.2222 28.6875
+14.2040 24.25 0
+0.25 n 2

raster/r.in.pdal/testsuite/test_r_in_pdal_binning.py

@@ -0,0 +1,363 @@
+"""
+Name:      r.in.pdal binning method test
+Purpose:   Validates output of various binning methods
+
+Author:    Maris Nartiss
+Copyright: (C) 2021 by Maris Nartiss and the GRASS Development Team
+Licence:   This program is free software under the GNU General Public
+           License (>=v2). Read the file COPYING that comes with GRASS
+           for details.
+"""
+
+import os
+import pathlib
+import unittest
+from tempfile import TemporaryDirectory
+
+from grass.script import core as grass
+from grass.script import shutil_which
+from grass.gunittest.case import TestCase
+from grass.gunittest.main import test
+
+
+class BinningTest(TestCase):
+    """Test binning methods
+
+    This test requires pdal CLI util to be available.
+    This tests expects r.in.ascii to work properly.
+    """
+
+    @classmethod
+    @unittest.skipIf(shutil_which("pdal") is None, "Cannot find pdal utility")
+    def setUpClass(cls):
+        """Ensures expected computational region and generated data"""
+        cls.use_temp_region()
+        cls.runModule("g.region", n=18, s=0, e=18, w=0, res=6)
+
+        cls.data_dir = os.path.join(pathlib.Path(__file__).parent.absolute(), "data")
+        cls.point_file = os.path.join(cls.data_dir, "points.csv")
+        cls.tmp_dir = TemporaryDirectory()
+        cls.las_file = os.path.join(cls.tmp_dir.name, "points.las")
+        grass.call(
+            [
+                "pdal",
+                "translate",
+                "-i",
+                cls.point_file,
+                "-o",
+                cls.las_file,
+                "-r",
+                "text",
+                "-w",
+                "las",
+                "--writers.las.format=0",
+                "--writers.las.extra_dims=all",
+                "--writers.las.minor_version=4",
+            ]
+        )
+
+    @classmethod
+    def tearDownClass(cls):
+        """Remove the temporary region and generated data"""
+        cls.tmp_dir.cleanup()
+        cls.del_temp_region()
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def tearDown(self):
+        """Remove the outputs created by the import
+
+        This is executed after each test run.
+        """
+        self.runModule("g.remove", flags="f", type="raster", name=self.bin_raster)
+        self.runModule("g.remove", flags="f", type="raster", name=self.ref_raster)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_method_n(self):
+        """Test binning with n method"""
+        self.bin_raster = "bin_n"
+        self.ref_raster = "ref_n"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.bin_raster,
+            flags="o",
+            quiet=True,
+            method="n",
+        )
+        self.assertRasterExists(self.bin_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_n_all.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.bin_raster, self.ref_raster, 0)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_method_min(self):
+        self.bin_raster = "bin_min"
+        self.ref_raster = "ref_min"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.bin_raster,
+            flags="o",
+            quiet=True,
+            method="min",
+        )
+        self.assertRasterExists(self.bin_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_min_z.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.bin_raster, self.ref_raster, 0)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_method_max(self):
+        self.bin_raster = "bin_max"
+        self.ref_raster = "ref_max"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.bin_raster,
+            flags="o",
+            quiet=True,
+            method="max",
+        )
+        self.assertRasterExists(self.bin_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_max_z.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.bin_raster, self.ref_raster, 0)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_method_range(self):
+        self.bin_raster = "bin_range"
+        self.ref_raster = "ref_range"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.bin_raster,
+            flags="o",
+            quiet=True,
+            method="range",
+        )
+        self.assertRasterExists(self.bin_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_range_z.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.bin_raster, self.ref_raster, 0)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_method_sum(self):
+        self.bin_raster = "bin_sum"
+        self.ref_raster = "ref_sum"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.bin_raster,
+            flags="o",
+            quiet=True,
+            method="sum",
+        )
+        self.assertRasterExists(self.bin_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_sum_z.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.bin_raster, self.ref_raster, 0)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_method_mean(self):
+        self.bin_raster = "bin_mean"
+        self.ref_raster = "ref_mean"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.bin_raster,
+            flags="o",
+            quiet=True,
+            method="mean",
+        )
+        self.assertRasterExists(self.bin_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_mean_z.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.bin_raster, self.ref_raster, 0.0001)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_method_stddev(self):
+        self.bin_raster = "bin_stddev"
+        self.ref_raster = "ref_stddev"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.bin_raster,
+            flags="o",
+            quiet=True,
+            method="stddev",
+        )
+        self.assertRasterExists(self.bin_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_stddev_z.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.bin_raster, self.ref_raster, 0.0001)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_method_variance(self):
+        self.bin_raster = "bin_variance"
+        self.ref_raster = "ref_variance"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.bin_raster,
+            flags="o",
+            quiet=True,
+            method="variance",
+        )
+        self.assertRasterExists(self.bin_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_variance_z.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.bin_raster, self.ref_raster, 0.0001)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_method_coeff_var(self):
+        self.bin_raster = "bin_coeff_var"
+        self.ref_raster = "ref_coeff_var"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.bin_raster,
+            flags="o",
+            quiet=True,
+            method="coeff_var",
+        )
+        self.assertRasterExists(self.bin_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_coeff_var_z.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.bin_raster, self.ref_raster, 0.0001)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_method_median(self):
+        self.bin_raster = "bin_median"
+        self.ref_raster = "ref_median"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.bin_raster,
+            flags="o",
+            quiet=True,
+            method="median",
+        )
+        self.assertRasterExists(self.bin_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_median_z.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.bin_raster, self.ref_raster, 0.0001)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_method_mode(self):
+        self.bin_raster = "bin_mode"
+        self.ref_raster = "ref_mode"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.bin_raster,
+            flags="o",
+            quiet=True,
+            method="mode",
+        )
+        self.assertRasterExists(self.bin_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_mode_z.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.bin_raster, self.ref_raster, 0)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_method_sidnmax(self):
+        self.bin_raster = "bin_sidnmax"
+        self.ref_raster = "ref_sidnmax"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.bin_raster,
+            flags="o",
+            quiet=True,
+            method="sidnmax",
+        )
+        self.assertRasterExists(self.bin_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_sidnmax_all.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.bin_raster, self.ref_raster, 0)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_method_sidnmin(self):
+        self.bin_raster = "bin_sidnmin"
+        self.ref_raster = "ref_sidnmin"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.bin_raster,
+            flags="o",
+            quiet=True,
+            method="sidnmin",
+        )
+        self.assertRasterExists(self.bin_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_sidnmin_all.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.bin_raster, self.ref_raster, 0)
+
+
+if __name__ == "__main__":
+    test()

raster/r.in.pdal/testsuite/test_r_in_pdal_selection.py

@@ -0,0 +1,187 @@
+"""
+Name:      r.in.pdal filter and selection test
+Purpose:   Validates output of input filtering and dimension selection
+
+Author:    Maris Nartiss
+Copyright: (C) 2020 by Maris Nartiss and the GRASS Development Team
+Licence:   This program is free software under the GNU General Public
+           License (>=v2). Read the file COPYING that comes with GRASS
+           for details.
+"""
+
+import os
+import pathlib
+import unittest
+from tempfile import TemporaryDirectory
+
+from grass.script import core as grass
+from grass.script import shutil_which
+from grass.gunittest.case import TestCase
+from grass.gunittest.main import test
+
+
+class SelectionTest(TestCase):
+    """Test input dimension selection and filtering
+
+    This test requires pdal CLI util to be available.
+    This tests expects r.in.ascii to work properly.
+    """
+
+    @classmethod
+    @unittest.skipIf(shutil_which("pdal") is None, "Cannot find pdal utility")
+    def setUpClass(cls):
+        """Ensures expected computational region and generated data"""
+        cls.use_temp_region()
+        cls.runModule("g.region", n=18, s=0, e=18, w=0, res=6)
+
+        cls.data_dir = os.path.join(pathlib.Path(__file__).parent.absolute(), "data")
+        cls.point_file = os.path.join(cls.data_dir, "points.csv")
+        cls.tmp_dir = TemporaryDirectory()
+        cls.las_file = os.path.join(cls.tmp_dir.name, "points.las")
+        grass.call(
+            [
+                "pdal",
+                "translate",
+                "-i",
+                cls.point_file,
+                "-o",
+                cls.las_file,
+                "-r",
+                "text",
+                "-w",
+                "las",
+                "--writers.las.format=0",
+                "--writers.las.extra_dims=all",
+                "--writers.las.minor_version=4",
+            ]
+        )
+
+    @classmethod
+    def tearDownClass(cls):
+        """Remove the temporary region and generated data"""
+        cls.tmp_dir.cleanup()
+        cls.del_temp_region()
+
+    def tearDown(self):
+        """Remove the outputs created by the import
+
+        This is executed after each test run.
+        """
+        self.runModule("g.remove", flags="f", type="raster", name=self.imp_raster)
+        self.runModule("g.remove", flags="f", type="raster", name=self.ref_raster)
+        try:
+            self.runModule("g.remove", flags="f", type="raster", name=self.base_raster)
+        except AttributeError:
+            pass
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_dimension(self):
+        """Test LAS dimension selection"""
+        self.imp_raster = "imp_intensity"
+        self.ref_raster = "ref_intensity"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.imp_raster,
+            flags="o",
+            quiet=True,
+            method="mean",
+            type="CELL",
+            dimension="intensity",
+        )
+        self.assertRasterExists(self.imp_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_mean_intensity.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.imp_raster, self.ref_raster, 0)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_user_dimension(self):
+        """Test PDAL user dimension selection"""
+        self.imp_raster = "imp_cellid"
+        self.ref_raster = "ref_cellid"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.imp_raster,
+            flags="o",
+            quiet=True,
+            method="mode",
+            type="CELL",
+            user_dimension="CellID",
+        )
+        self.assertRasterExists(self.imp_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_mode_cellid.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.imp_raster, self.ref_raster, 0)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_filter(self):
+        """Test input filtering"""
+        self.imp_raster = "imp_filtered"
+        self.ref_raster = "ref_filtered"
+
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.imp_raster,
+            flags="o",
+            quiet=True,
+            method="mode",
+            type="CELL",
+            dimension="source",
+            zrange=(2, 10),
+            irange=(10, 20),
+            drange=(1, 1),
+        )
+        self.assertRasterExists(self.imp_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_filter_z_int_source.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.imp_raster, self.ref_raster, 0)
+
+    @unittest.skipIf(shutil_which("r.in.pdal") is None, "Cannot find r.in.pdal")
+    def test_base_raster(self):
+        """Test Z adjustement by base raster"""
+        self.imp_raster = "imp_base_adj"
+        self.ref_raster = "ref_base_adj"
+        self.base_raster = "base_raster"
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_mean_z.ascii"),
+            output=self.base_raster,
+        )
+        self.assertModule(
+            "r.in.pdal",
+            input=self.las_file,
+            output=self.imp_raster,
+            flags="o",
+            quiet=True,
+            method="max",
+            base_raster=self.base_raster,
+        )
+        self.assertRasterExists(self.imp_raster)
+
+        self.runModule(
+            "r.in.ascii",
+            input=os.path.join(self.data_dir, "res_base_adj.ascii"),
+            output=self.ref_raster,
+        )
+        self.assertRastersEqual(self.imp_raster, self.ref_raster, 0.001)
+
+
+if __name__ == "__main__":
+    test()

raster/rasterintro.html

@@ -164,7 +164,7 @@ Furthermore, there are modules available for reinterpolation of "sparse"
 <li> Various vector modules for interpolation</li>
 </ul>
 
-For Lidar and similar data, <a href="r.in.lidar.html">r.in.lidar</a> and <a href="r.in.xyz.html">r.in.xyz</a>
+For Lidar and similar data, <a href="r.in.pdal.html">r.in.pdal</a> and <a href="r.in.xyz.html">r.in.xyz</a>
 support loading and binning of ungridded x,y,z ASCII data into a new raster map.
 The user may choose from a variety of statistical methods in creating the new raster map.
 <p>