v.in.pdal: LiDAR import tool based on PDAL [news]

Contains most of the v.in.lidar filters/selection tools.
Supports two advaced PDAL filters.

Tests for basic functonality.

Doesn't build topology. Stores
return, class and color info in predefined
layers as categories.


Tests and other .c files are similar
to those in v.in.lidar.

See also https://trac.osgeo.org/grass/ticket/2732 and r67293.


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

vector/Makefile

@@ -28,6 +28,7 @@ SUBDIRS = \
 	v.in.ascii \
 	v.in.db \
 	v.in.dxf \
+	v.in.pdal \
 	v.in.region \
 	v.kcv \
 	v.kernel \

vector/v.in.pdal/Makefile

@@ -0,0 +1,19 @@
+MODULE_TOPDIR = ../..
+
+PGM=v.in.pdal
+
+LIBES = $(GPROJLIB) $(VECTORLIB) $(DBMILIB) $(GISLIB) $(MATHLIB) $(PDALLIBS)
+DEPENDENCIES = $(GPROJDEP) $(VECTORDEP) $(DBMIDEP) $(GISDEP)
+
+EXTRA_INC = $(VECT_INC) $(PROJINC) $(PDALINC)
+EXTRA_CFLAGS = $(VECT_CFLAGS) $(PDALCPPFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+LINK = $(CXX)
+
+ifneq ($(strip $(CXX)),)
+ifneq ($(USE_PDAL),)
+default: cmd
+endif
+endif

vector/v.in.pdal/filters.c

@@ -0,0 +1,156 @@
+/*
+ * v.in.lidar filtering functions
+ *
+ * Copyright 2011-2015 by Markus Metz, and The GRASS Development Team
+ * Authors:
+ *  Markus Metz (r.in.lidar)
+ *  Vaclav Petras (move code to a separate files)
+ *
+ * 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 zrange_filter_from_option(struct Option *option, double *zmin,
+                              double *zmax)
+{
+    if (option->answer != NULL) {
+        if (option->answers[0] == NULL || option->answers[1] == NULL)
+            G_fatal_error(_("Invalid zrange <%s>"), option->answer);
+        sscanf(option->answers[0], "%lf", zmin);
+        sscanf(option->answers[1], "%lf", zmax);
+        /* for convenience, switch order to make valid input */
+        if (*zmin > *zmax) {
+            double tmp = *zmax;
+
+            *zmax = *zmin;
+            *zmin = 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;
+}

vector/v.in.pdal/filters.h

@@ -0,0 +1,47 @@
+/*
+ * v.in.lidar filtering functions
+ *
+ * Copyright 2011-2015 by Markus Metz, and The GRASS Development Team
+ * Authors:
+ *  Markus Metz (r.in.lidar)
+ *  Vaclav Petras (move code to a separate files)
+ *
+ * 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 zrange_filter_from_option(struct Option *option,
+                              double *zmin, double *zmax);
+
+int return_filter_create_from_string(struct ReturnFilter *return_filter,
+                                     const char *name);
+int return_filter_is_out(struct ReturnFilter *return_filter, int return_n,
+                         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__ */

vector/v.in.pdal/lidar.c

@@ -0,0 +1,77 @@
+
+/****************************************************************************
+ *
+ * MODULE:       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_default_layers(struct GLidarLayers *layers)
+{
+    layers->id_layer = 0;
+    layers->return_layer = G_RETURN_LAYER;
+    layers->n_returns_layer = G_NUM_RETURNS_LAYER;
+    layers->class_layer = G_CLASS_LAYER;
+    layers->rgb_layer = G_RGB_LAYER;
+    layers->red_layer = 0;
+    layers->green_layer = 0;
+    layers->blue_layer = 0;
+}
+
+void GLidarLayers_set_all_layers(struct GLidarLayers *layers)
+{
+    layers->id_layer = G_ID_LAYER;
+    layers->return_layer = G_RETURN_LAYER;
+    layers->n_returns_layer = G_NUM_RETURNS_LAYER;
+    layers->class_layer = G_CLASS_LAYER;
+    layers->rgb_layer = G_RGB_LAYER;
+    layers->red_layer = G_RED_LAYER;
+    layers->green_layer = G_GREEN_LAYER;
+    layers->blue_layer = G_BLUE_LAYER;
+}
+
+void GLidarLayers_set_no_layers(struct GLidarLayers *layers)
+{
+    layers->id_layer = 0;
+    layers->return_layer = 0;
+    layers->n_returns_layer = 0;
+    layers->class_layer = 0;
+    layers->rgb_layer = 0;
+    layers->red_layer = 0;
+    layers->green_layer = 0;
+    layers->blue_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}
+};

vector/v.in.pdal/lidar.h

@@ -0,0 +1,94 @@
+
+/****************************************************************************
+ *
+ * MODULE:       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_LAST 2
+#define LAS_MID 3
+
+#define G_ID_LAYER 1
+#define G_RETURN_LAYER 2
+#define G_NUM_RETURNS_LAYER 3
+#define G_CLASS_LAYER 4
+#define G_RGB_LAYER 5
+#define G_RED_LAYER 6
+#define G_GREEN_LAYER 7
+#define G_BLUE_LAYER 8
+
+struct GLidarLayers
+{
+    int id_layer;
+    int return_layer;
+    int n_returns_layer;
+    int class_layer;
+    int rgb_layer;
+    int red_layer;
+    int green_layer;
+    int blue_layer;
+};
+
+void GLidarLayers_set_default_layers(struct GLidarLayers *layers);
+void GLidarLayers_set_all_layers(struct GLidarLayers *layers);
+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[];
+
+#endif /* GRASS_LIDAR_H */

vector/v.in.pdal/main.cpp

@@ -0,0 +1,547 @@
+#include <pdal/PointTable.hpp>
+#include <pdal/PointView.hpp>
+#include <pdal/StageFactory.hpp>
+#include <pdal/LasReader.hpp>
+#include <pdal/LasHeader.hpp>
+#include <pdal/Options.hpp>
+#include <pdal/ReprojectionFilter.hpp>
+
+extern "C"
+{
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/gprojects.h>
+#include <grass/glocale.h>
+}
+
+extern "C"
+{
+#include "lidar.h"
+#include "projection.h"
+#include "filters.h"
+}
+
+#ifdef HAVE_LONG_LONG_INT
+typedef unsigned long long gpoint_count;
+#else
+typedef unsigned long gpoint_count;
+#endif
+
+
+void pdal_point_to_grass(struct Map_info *output_vector,
+                         struct line_pnts *points, struct line_cats *cats,
+                         pdal::PointViewPtr point_view, pdal::PointId idx,
+                         struct GLidarLayers *layers, int cat,
+                         pdal::Dimension::Id::Enum dim_to_use_as_z)
+{
+    Vect_reset_line(points);
+    Vect_reset_cats(cats);
+
+    using namespace pdal::Dimension::Id;
+    double x = point_view->getFieldAs<double>(X, idx);
+    double y = point_view->getFieldAs<double>(Y, idx);
+    double z = point_view->getFieldAs<double>(dim_to_use_as_z, idx);
+
+    if (layers->id_layer) {
+        Vect_cat_set(cats, layers->id_layer, cat);
+    }
+    if (layers->return_layer) {
+        Vect_cat_set(cats, layers->return_layer,
+                     point_view->getFieldAs<int>(ReturnNumber, idx));
+    }
+    if (layers->n_returns_layer) {
+        Vect_cat_set(cats, layers->n_returns_layer,
+                     point_view->getFieldAs<int>(NumberOfReturns, idx));
+    }
+    if (layers->class_layer) {
+        Vect_cat_set(cats, layers->class_layer,
+                     point_view->getFieldAs<int>(Classification, idx));
+    }
+    if (layers->rgb_layer) {
+        int red = point_view->getFieldAs<int>(Red, idx);
+        int green = point_view->getFieldAs<int>(Green, idx);
+        int blue = point_view->getFieldAs<int>(Blue, idx);
+        int rgb = red;
+        rgb = (rgb << 8) + green;
+        rgb = (rgb << 8) + blue;
+        Vect_cat_set(cats, layers->rgb_layer, rgb);
+    }
+    if (layers->red_layer) {
+        Vect_cat_set(cats, layers->red_layer,
+                     point_view->getFieldAs<int>(Red, idx));
+    }
+    if (layers->green_layer) {
+        Vect_cat_set(cats, layers->green_layer,
+                     point_view->getFieldAs<int>(Green, idx));
+    }
+    if (layers->blue_layer) {
+        Vect_cat_set(cats, layers->blue_layer,
+                     point_view->getFieldAs<int>(Blue, idx));
+    }
+
+    Vect_append_point(points, x, y, z);
+    Vect_write_line(output_vector, GV_POINT, points, cats);
+}
+
+
+int main(int argc, char *argv[])
+{
+    G_gisinit(argv[0]);
+
+    GModule *module = G_define_module();
+    G_add_keyword(_("vector"));
+    G_add_keyword(_("import"));
+    G_add_keyword(_("LIDAR"));
+    module->description =
+        _("Converts LAS LiDAR point clouds to a GRASS vector map with PDAL.");
+
+    Option *in_opt = G_define_standard_option(G_OPT_F_INPUT);
+    in_opt->label = _("LAS input file");
+    in_opt->description =
+        _("LiDAR input files in LAS format (*.las or *.laz)");
+
+    Option *out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+
+    Option *spatial_opt = G_define_option();
+    spatial_opt->key = "spatial";
+    spatial_opt->type = TYPE_DOUBLE;
+    spatial_opt->multiple = NO;
+    spatial_opt->required = NO;
+    // TODO: does this require multiple or not?
+    spatial_opt->key_desc = "xmin,ymin,xmax,ymax";
+    spatial_opt->label = _("Import subregion only");
+    spatial_opt->description =
+        _("Format: xmin,ymin,xmax,ymax - usually W,S,E,N");
+    spatial_opt->guisection = _("Selection");
+
+    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->description = _("Filter range for z data (min,max)");
+    zrange_opt->guisection = _("Selection");
+
+    Option *filter_opt = G_define_option();
+    filter_opt->key = "return_filter";
+    filter_opt->type = TYPE_STRING;
+    filter_opt->required = NO;
+    filter_opt->label = _("Only import points of selected return type");
+    filter_opt->description = _("If not specified, all points are imported");
+    filter_opt->options = "first,last,mid";
+    filter_opt->guisection = _("Selection");
+
+    Option *class_opt = G_define_option();
+    class_opt->key = "class_filter";
+    class_opt->type = TYPE_INTEGER;
+    class_opt->multiple = YES;
+    class_opt->required = NO;
+    class_opt->label = _("Only import points of selected class(es)");
+    class_opt->description = _("Input is comma separated integers. "
+                               "If not specified, all points are imported.");
+    class_opt->guisection = _("Selection");
+
+    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");
+
+    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");
+
+    // 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");
+
+    Option *max_ground_window_opt = G_define_option();
+    max_ground_window_opt->key = "max_ground_window_size";
+    max_ground_window_opt->type = TYPE_DOUBLE;
+    max_ground_window_opt->required = NO;
+    max_ground_window_opt->answer = "33";
+    max_ground_window_opt->description =
+        _("Maximum window size for ground filter");
+    max_ground_window_opt->guisection = _("Ground filter");
+
+    Option *ground_slope_opt = G_define_option();
+    ground_slope_opt->key = "ground_slope";
+    ground_slope_opt->type = TYPE_DOUBLE;
+    ground_slope_opt->required = NO;
+    ground_slope_opt->answer = "1.0";
+    ground_slope_opt->description = _("Slope for ground filter");
+    ground_slope_opt->guisection = _("Ground filter");
+
+    Option *max_ground_distance_opt = G_define_option();
+    max_ground_distance_opt->key = "max_ground_distance";
+    max_ground_distance_opt->type = TYPE_DOUBLE;
+    max_ground_distance_opt->required = NO;
+    max_ground_distance_opt->answer = "2.5";
+    max_ground_distance_opt->description =
+        _("Maximum distance for ground filter");
+    max_ground_distance_opt->guisection = _("Ground filter");
+
+    Option *init_ground_distance_opt = G_define_option();
+    init_ground_distance_opt->key = "initial_ground_distance";
+    init_ground_distance_opt->type = TYPE_DOUBLE;
+    init_ground_distance_opt->required = NO;
+    init_ground_distance_opt->answer = "0.15";
+    init_ground_distance_opt->description =
+        _("Initial distance for ground filter");
+    init_ground_distance_opt->guisection = _("Ground filter");
+
+    Option *ground_cell_size_opt = G_define_option();
+    ground_cell_size_opt->key = "ground_cell_size";
+    ground_cell_size_opt->type = TYPE_DOUBLE;
+    ground_cell_size_opt->required = NO;
+    ground_cell_size_opt->answer = "1";
+    ground_cell_size_opt->description =
+        _("Initial distance for ground filter");
+    ground_cell_size_opt->guisection = _("Ground filter");
+
+    Flag *region_flag = G_define_flag();
+    region_flag->key = 'r';
+    region_flag->guisection = _("Selection");
+    region_flag->description = _("Limit import to the current region");
+
+    Flag *nocats_flag = G_define_flag();
+    nocats_flag->key = 'c';
+    nocats_flag->label =
+        _("Store only the coordinates");
+    nocats_flag->description =
+        _("Do not add any categories to points");
+    nocats_flag->guisection = _("Categories");
+
+    Flag *idcat_flag = G_define_flag();
+    idcat_flag->key = 'd';
+    idcat_flag->label =
+        _("Store generated unique ID for each point");
+    idcat_flag->guisection = _("Categories");
+
+    Flag *extract_ground_flag = G_define_flag();
+    extract_ground_flag->key = 'j';
+    extract_ground_flag->label =
+        _("Classify and extract ground points");
+    extract_ground_flag->description =
+        _("This assignes class 2 to the groud points");
+    extract_ground_flag->guisection = _("Ground filter");
+
+    // TODO: by inverting class filter and choosing 2 we can select non-groud points
+    // this can be done as a separate flag (generally useful?)
+    // or this flag can be changed (only ground is classified anyway)
+    // and it would Classify ground and extract non-ground
+    // probably better if only one step is required to get ground and non-ground
+    Flag *classify_ground_flag = G_define_flag();
+    classify_ground_flag->key = 'k';
+    classify_ground_flag->description = _("Classify ground points");
+    classify_ground_flag->guisection = _("Ground filter");
+
+    Flag *height_filter_flag = G_define_flag();
+    height_filter_flag->key = 'h';
+    height_filter_flag->label =
+        _("Compute height for points as a difference from ground");
+    height_filter_flag->description =
+        _("This requires points to have class 2");
+    height_filter_flag->guisection = _("Transform");
+
+    Flag *approx_ground_flag = G_define_flag();
+    approx_ground_flag->key = 'm';
+    approx_ground_flag->description =
+        _("Use approximate algorithm in ground filter");
+    approx_ground_flag->guisection = _("Ground filter");
+
+    G_option_exclusive(spatial_opt, region_flag, NULL);
+    G_option_exclusive(reproject_flag, over_flag, NULL);
+    G_option_exclusive(nocats_flag, idcat_flag, NULL);
+    G_option_exclusive(extract_ground_flag, classify_ground_flag, NULL);
+
+    if (G_parser(argc, argv))
+        return EXIT_FAILURE;
+
+    if (access(in_opt->answer, F_OK) != 0) {
+        G_fatal_error(_("Input file <%s> does not exist"), in_opt->answer);
+    }
+
+    // we use full qualification because the dim ns conatins too general names
+    pdal::Dimension::Id::Enum dim_to_use_as_z = pdal::Dimension::Id::Z;
+
+    double xmin = 0;
+    double ymin = 0;
+    double xmax = 0;
+    double ymax = 0;
+    bool use_spatial_filter = false;
+    if (spatial_opt->answer)
+        use_spatial_filter = spatial_filter_from_option(spatial_opt,
+                                                        &xmin, &ymin,
+                                                        &xmax, &ymax);
+    else if (region_flag->answer)
+        use_spatial_filter = spatial_filter_from_current_region(&xmin,
+                                                                &ymin,
+                                                                &xmax,
+                                                                &ymax);
+
+    double zrange_min, zrange_max;
+    bool use_zrange = zrange_filter_from_option(zrange_opt, &zrange_min,
+                                                &zrange_max);
+    struct ReturnFilter return_filter_struct;
+    bool use_return_filter =
+        return_filter_create_from_string(&return_filter_struct,
+                                         filter_opt->answer);
+    struct ClassFilter class_filter;
+    bool use_class_filter =
+        class_filter_create_from_strings(&class_filter, class_opt->answers);
+
+    pdal::StageFactory factory;
+    std::string pdal_read_driver = factory.inferReaderDriver(in_opt->answer);
+    if (pdal_read_driver.empty())
+        G_fatal_error("Cannot determine input file type of <%s>",
+                      in_opt->answer);
+
+    pdal::Option las_opt("filename", in_opt->answer);
+    pdal::Options las_opts;
+    las_opts.add(las_opt);
+    // TODO: free reader
+    // using plain pointer because we need to keep the last stage pointer
+    pdal::Stage * reader = factory.createStage(pdal_read_driver);
+    if (!reader)
+        G_fatal_error("PDAL reader creation failed, a wrong format of <%s>",
+                      in_opt->answer);
+    reader->setOptions(las_opts);
+
+    pdal::Stage * last_stage = reader;
+    pdal::ReprojectionFilter reprojection_filter;
+
+    // we reproject when requested regardless 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(*reader);
+        last_stage = &reprojection_filter;
+    }
+
+    if (extract_ground_flag->answer || classify_ground_flag->answer) {
+        if (extract_ground_flag->answer)
+            G_message(_("Extracting ground points"));
+        if (classify_ground_flag->answer)
+            G_message(_("Classifying ground points"));
+        pdal::Options groundOptions;
+        groundOptions.add<double>("max_window_size",
+                                  atof(max_ground_window_opt->answer));
+        groundOptions.add<double>("slope",
+                                  atof(ground_slope_opt->answer));
+        groundOptions.add<double>("max_distance",
+                                  atof(max_ground_distance_opt->answer));
+        groundOptions.add<double>("initial_distance",
+                                  atof(init_ground_distance_opt->answer));
+        groundOptions.add<double>("cell_size",
+                                  atof(ground_cell_size_opt->answer));
+        groundOptions.add<bool>("classify",
+                                classify_ground_flag->answer);
+        groundOptions.add<bool>("extract",
+                                extract_ground_flag->answer);
+        groundOptions.add<bool>("approximate",
+                                approx_ground_flag->answer);
+        groundOptions.add<bool>("debug", false);
+        groundOptions.add<uint32_t>("verbose", 0);
+
+        // TODO: free this or change pointer type to shared
+        pdal::Stage * ground_stage(factory.createStage("filters.ground"));
+        if (!ground_stage)
+            G_fatal_error(_("Ground filter is not available"
+                            " (PDAL probably compiled without PCL)"));
+        ground_stage->setOptions(groundOptions);
+        ground_stage->setInput(*last_stage);
+        last_stage = ground_stage;
+    }
+
+    if (height_filter_flag->answer) {
+        // TODO: we should test with if (point_view->hasDim(Id::Classification))
+        // but we don't have the info yet
+        // TODO: free this or change pointer type to shared
+        pdal::Stage * height_stage(factory.createStage("filters.height"));
+        if (!height_stage)
+            G_fatal_error(_("Height above ground filter is not available"
+                            " (PDAL probably compiled without PCL)"));
+        height_stage->setInput(*last_stage);
+        last_stage = height_stage;
+    }
+
+    pdal::PointTable point_table;
+    last_stage->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 = reader->getSpatialReference();
+        if (spatial_reference.empty())
+            G_fatal_error(_("The input dataset has undefined projection"));
+        std::string dataset_wkt =
+            spatial_reference.
+            getWKT(pdal::SpatialReference::eHorizontalOnly);
+        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..."));
+    pdal::PointViewSet point_view_set = last_stage->execute(point_table);
+    pdal::PointViewPtr point_view = *point_view_set.begin();
+
+    // 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
+
+    struct GLidarLayers layers;
+    GLidarLayers_set_default_layers(&layers);
+
+    // update layers we are writting based on what is in the data
+    // update usage of our filters as well
+    if (point_view->hasDim(pdal::Dimension::Id::ReturnNumber) &&
+        point_view->hasDim(pdal::Dimension::Id::NumberOfReturns)) {
+        use_return_filter = true;
+    }
+    else {
+        layers.return_layer = 0;
+        layers.n_returns_layer = 0;
+        use_return_filter = false;
+    }
+
+    if (point_view->hasDim(pdal::Dimension::Id::Classification)) {
+        use_class_filter = true;
+    }
+    else {
+        layers.class_layer = 0;
+        use_class_filter = false;
+    }
+
+    if (!(point_view->hasDim(pdal::Dimension::Id::Red) &&
+            point_view->hasDim(pdal::Dimension::Id::Green) &&
+            point_view->hasDim(pdal::Dimension::Id::Blue) &&
+            (layers.rgb_layer || layers.red_layer
+             || layers.green_layer || layers.blue_layer))) {
+        layers.rgb_layer = 0;
+        layers.red_layer = 0;
+        layers.green_layer = 0;
+        layers.blue_layer = 0;
+    }
+    if (idcat_flag->answer)
+        layers.id_layer = G_ID_LAYER;
+
+    // we just force it without any checking (we rely on the options)
+    if (nocats_flag->answer) {
+        GLidarLayers_set_no_layers(&layers);
+    }
+
+    G_important_message(_("Scanning points..."));
+    struct Map_info output_vector;
+
+    // the overwrite warning comes quite late in the execution
+    // but that's good enough
+    if (Vect_open_new(&output_vector, out_opt->answer, 1) < 0)
+        G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
+    Vect_hist_command(&output_vector);
+
+    // height is stored as a new attribute
+    if (height_filter_flag->answer) {
+        dim_to_use_as_z = point_view->layout()->findDim("Height");
+        if (dim_to_use_as_z == pdal::Dimension::Id::Unknown)
+            G_fatal_error(_("Cannot indentify the height dimension"
+                            " (probably something changed in PDAL)"));
+    }
+
+    // this is just for sure, we test the individual dimensions before
+    // TODO: should we test Z explicitly as well?
+    if (!point_view->hasDim(dim_to_use_as_z))
+        G_fatal_error(_("Dataset doesn't have requested dimension '%s'"
+                        " with ID %d (possibly a programming error)"),
+                      pdal::Dimension::name(dim_to_use_as_z).c_str(),
+                      dim_to_use_as_z);
+
+    struct line_pnts *points = Vect_new_line_struct();
+    struct line_cats *cats = Vect_new_cats_struct();
+
+    gpoint_count n_outside = 0;
+    gpoint_count zrange_filtered = 0;
+    gpoint_count n_filtered = 0;
+    gpoint_count n_class_filtered = 0;
+
+    int cat = 1;
+    bool cat_max_reached = false;
+
+    for (pdal::PointId idx = 0; idx < point_view->size(); ++idx) {
+        // TODO: avoid duplication of reading the attributes here and when writing if needed
+        double x = point_view->getFieldAs<double>(pdal::Dimension::Id::X, idx);
+        double y = point_view->getFieldAs<double>(pdal::Dimension::Id::Y, idx);
+        double z = point_view->getFieldAs<double>(dim_to_use_as_z, idx);
+
+        if (use_spatial_filter) {
+            if (x < xmin || x > xmax || y < ymin || y > ymax) {
+                n_outside++;
+                continue;
+            }
+        }
+        if (use_zrange) {
+            if (z < zrange_min || z > zrange_max) {
+                zrange_filtered++;
+                continue;
+            }
+        }
+        if (use_return_filter) {
+            int return_n =
+                point_view->getFieldAs<int>(pdal::Dimension::Id::ReturnNumber, idx);
+            int n_returns =
+                point_view->getFieldAs<int>(pdal::Dimension::Id::NumberOfReturns, idx);
+            if (return_filter_is_out
+                (&return_filter_struct, return_n, n_returns)) {
+                n_filtered++;
+                continue;
+            }
+        }
+        if (use_class_filter) {
+            int point_class =
+                point_view->getFieldAs<int>(pdal::Dimension::Id::Classification, idx);
+            if (class_filter_is_out(&class_filter, point_class)) {
+                n_class_filtered++;
+                continue;
+            }
+        }
+        pdal_point_to_grass(&output_vector, points, cats, point_view,
+                            idx, &layers, cat, dim_to_use_as_z);
+        if (layers.id_layer) {
+            // TODO: perhaps it would be better to use the max cat afterwords
+            if (cat == GV_CAT_MAX) {
+                cat_max_reached = true;
+                break;
+            }
+            cat++;
+        }
+    }
+    // not building topology by default
+    Vect_close(&output_vector);
+}

vector/v.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;
+}

vector/v.in.pdal/projection.h

@@ -0,0 +1,35 @@
+
+/****************************************************************************
+ *
+ * MODULE:       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 */

vector/v.in.pdal/testsuite/basic_test.py

@@ -0,0 +1,71 @@
+"""
+Name:      basic_test
+Purpose:   v.in.pdal basic functionality test
+
+Author:    Vaclav Petras
+Copyright: (C) 2015 by Vaclav Petras 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
+from grass.gunittest.case import TestCase
+from grass.gunittest.main import test
+
+
+class BasicTest(TestCase):
+    """Test case for watershed module
+
+    This tests expects v.random and v.out.lidar to work properly.
+    """
+
+    # Setup variables to be used for outputs
+    vector_generated = 'vinlidar_basic_generated'
+    vector_points = 'vinlidar_basic_original'
+    imported_points = 'vinlidar_basic_imported'
+    las_file = 'vinlidar_basic_points.las'
+    npoints = 300
+
+    @classmethod
+    def setUpClass(cls):
+        """Ensures expected computational region and generated data"""
+        cls.use_temp_region()
+        cls.runModule('g.region', n=20, s=10, e=25, w=15, res=1)
+        cls.runModule('v.random', flags='zb', output=cls.vector_generated,
+            npoints=cls.npoints, zmin=200, zmax=500, seed=100)
+        cls.runModule('v.category', input=cls.vector_generated,
+            output=cls.vector_points, option='del', cat=-1)
+        cls.runModule('v.out.lidar', input=cls.vector_points,
+            output=cls.las_file)
+
+    @classmethod
+    def tearDownClass(cls):
+        """Remove the temporary region and generated data"""
+        cls.runModule('g.remove', flags='f', type='vector',
+            name=(cls.vector_points, cls.vector_generated))
+        if os.path.isfile(cls.las_file):
+            os.remove(cls.las_file)
+        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='vector',
+            name=self.imported_points)
+
+    def test_same_data(self):
+        """Test to see if the standard outputs are created"""
+        self.assertModule('v.in.pdal', input=self.las_file, flags='c',
+            output=self.imported_points)
+        self.assertVectorExists(self.imported_points)
+        self.assertVectorEqualsVector(
+            actual=self.imported_points,
+            reference=self.vector_points,
+            digits=2, precision=.01)
+
+
+if __name__ == '__main__':
+    test()

vector/v.in.pdal/testsuite/filter_test.py

@@ -0,0 +1,208 @@
+"""
+Name:      filter_test
+Purpose:   v.in.pdal test if various filters and selections
+
+Author:    Vaclav Petras
+Copyright: (C) 2015 by Vaclav Petras 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
+from grass.gunittest.case import TestCase
+from grass.gunittest.main import test
+
+
+POINTS = """\
+17.46938776,18.67346939,143,1,1,2
+20.93877551,17.44897959,125,1,1,2
+18.89795918,14.18367347,130,1,1,3
+15.91836735,10.67346939,126,1,1,3
+21.26530612,11.04081633,128,1,2,3
+22.24489796,13.89795918,123,2,2,3
+23.79591837,17.12244898,151,1,2,3
+17.2244898,16.34693878,124,2,2,4
+17.14285714,14.10204082,134,1,3,4
+19.87755102,11.81632653,146,2,3,4
+18.48979592,11.48979592,140.6,2,3,4
+21.26530612,15.73469388,147,3,3,5
+21.18367347,19.32653061,138,1,3,5
+23.91836735,18.83673469,144,2,3,5
+23.51020408,13.65306122,143,3,3,5
+23.55102041,11.32653061,123,1,4,5
+18.41009273,14.51618034,140.4,2,4,5
+22.13996161,17.2278263,147,3,4,5
+21.41013052,11.05432488,132,4,4,5
+"""
+
+
+class FilterTest(TestCase):
+    """Test case for filter and selection options
+
+    This tests expects v.random and v.out.lidar to work properly.
+    """
+
+    # Setup variables to be used for outputs
+    vector_points = 'vinlidar_filters_original'
+    imported_points = 'vinlidar_filters_imported'
+    las_file = 'vinlidar_filters_points.las'
+    npoints = 300
+
+    @classmethod
+    def setUpClass(cls):
+        """Ensures expected computational region and generated data"""
+        cls.use_temp_region()
+        cls.runModule('g.region', n=20, s=10, e=25, w=15, res=1)
+        cls.runModule('v.in.ascii', input='-', stdin_=POINTS,
+                      flags='z', z=3, cat=0, separator='comma',
+                      output=cls.vector_points,
+                      columns="x double precision, y double precision,"
+                              " z double precision, return_n integer,"
+                              " n_returns integer, class_n integer")
+        cls.runModule('v.out.lidar',
+                      input=cls.vector_points, layer=1,
+                      output=cls.las_file,
+                      return_column='return_n',
+                      n_returns_column='n_returns',
+                      class_column='class_n')
+
+    @classmethod
+    def tearDownClass(cls):
+        """Remove the temporary region and generated data"""
+        cls.runModule('g.remove', flags='f', type='vector',
+            name=cls.vector_points)
+        if os.path.isfile(cls.las_file):
+            os.remove(cls.las_file)
+        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='vector',
+            name=self.imported_points)
+
+    def test_no_filter(self):
+        """Test to see if the standard outputs are created
+
+        This shows if the inpute data are as expected.
+        """
+        self.assertModule('v.in.pdal', input=self.las_file,
+            output=self.imported_points)
+        self.assertVectorExists(self.imported_points)
+        self.assertVectorFitsTopoInfo(
+            vector=self.imported_points,
+            reference=dict(points=19))
+
+    def return_filter(self, name, npoints):
+        """Mid return filter test"""
+        self.assertModule('v.in.pdal', input=self.las_file,
+            output=self.imported_points,
+            return_filter=name)
+        self.assertVectorExists(self.imported_points)
+        self.assertVectorFitsTopoInfo(
+            vector=self.imported_points,
+            reference=dict(points=npoints))
+
+    def test_first_return_filter(self):
+        """First return filter test"""
+        self.return_filter('first', 9)
+
+    def test_mid_return_filter(self):
+        """Mid return filter test"""
+        self.return_filter('mid', 5)
+
+    def test_last_return_filter(self):
+        """Last return filter test"""
+        self.return_filter('last', 5)
+
+    def class_filter(self, class_n, npoints):
+        """Actual code for testing class filter"""
+        self.assertModule('v.in.pdal', input=self.las_file,
+            output=self.imported_points,
+            class_filter=class_n)
+        self.assertVectorExists(self.imported_points)
+        self.assertVectorFitsTopoInfo(
+            vector=self.imported_points,
+            reference=dict(points=npoints))
+
+    def test_class_2_filter(self):
+        """Test to filter classes"""
+        self.class_filter(2, 2)
+
+    def test_class_3_filter(self):
+        """Test to filter classes"""
+        self.class_filter(3, 5)
+
+    def test_class_4_filter(self):
+        """Test to filter classes"""
+        self.class_filter(4, 4)
+
+    def test_class_5_filter(self):
+        """Test to filter classes"""
+        self.class_filter(5, 8)
+
+    def return_and_class_filter(self, return_name, class_n, npoints):
+        """Return and class filter combined test code"""
+        self.assertModule('v.in.pdal', input=self.las_file,
+            output=self.imported_points,
+            return_filter=return_name, class_filter=class_n)
+        self.assertVectorExists(self.imported_points)
+        self.assertVectorFitsTopoInfo(
+            vector=self.imported_points,
+            reference=dict(points=npoints))
+
+    def test_first_return_and_class_filter(self):
+        """Combined test for return and class"""
+        self.return_and_class_filter('first', 2, 2)
+
+    def test_last_return_and_class_filter(self):
+        """Combined test for return and class"""
+        self.return_and_class_filter('last', 5, 3)
+
+    def zrange_filter(self, zrange, npoints):
+        """Actual code for zrange option test"""
+        self.assertModule('v.in.pdal', input=self.las_file,
+            output=self.imported_points,
+            zrange=zrange)
+        self.assertVectorExists(self.imported_points)
+        self.assertVectorFitsTopoInfo(
+            vector=self.imported_points,
+            reference=dict(points=npoints))
+
+    def test_zrange_filter(self):
+        """Test zrange option"""
+        self.zrange_filter((130.1, 139.9), 3)
+
+    def test_non_int_zrange_filter(self):
+        """Test zrange option with float number
+
+        One test point has z right under and one other right above the min.
+        """
+        self.zrange_filter((140.5, 900), 8)
+
+    def test_zrange_and_class_filter(self):
+        """zrange and class_filter option combined test"""
+        self.assertModule('v.in.pdal', input=self.las_file,
+            output=self.imported_points,
+            zrange=(141, 900), class_filter=5)
+        self.assertVectorExists(self.imported_points)
+        self.assertVectorFitsTopoInfo(
+            vector=self.imported_points,
+            reference=dict(points=4))
+
+    def test_zrange_and_return_filter(self):
+        """zrange and class_filter option combined test"""
+        self.assertModule('v.in.pdal', input=self.las_file,
+            output=self.imported_points,
+            zrange=(141, 900), return_filter='last')
+        self.assertVectorExists(self.imported_points)
+        self.assertVectorFitsTopoInfo(
+            vector=self.imported_points,
+            reference=dict(points=2))
+
+
+if __name__ == '__main__':
+    test()

vector/v.in.pdal/v.in.pdal.html

@@ -0,0 +1,39 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.in.pdal</em> supports the following PDAL filters:
+
+<ul>
+<li>groud detection with Point Cloud Library (PCL) Progressive Morphological Filter (Zhang 2003)
+    (<a href="http://www.pdal.io/stages/filters.ground.html">filters.ground</a>)
+<li><a href="http://www.pdal.io/stages/filters.height.html">filters.height</a>
+    (requires PCL)
+</ul>
+
+<em>v.in.pdal</em> supports the following filters natively:
+
+<ul>
+<li>2D region filter
+<li>Z coordinates filter
+<li>return filter
+<li>class filter
+</ul>
+
+<h2>EXAMPLES</h2>
+
+Import only XYZ coordinates of points, limit the import to the current
+computational region and reproject to the Location projection during
+the import:
+<div class="code"><pre>
+v.in.pdal input=points.las output=points -c -r -w
+</pre></div>
+
+
+<h2>REFERENCES</h2>
+
+<ul>
+<li>Zhang, Keqi, et al.
+    <em>A progressive morphological filter for removing nonground
+    measurements from airborne LIDAR data.</em>
+    Geoscience and Remote Sensing,
+    IEEE Transactions on 41.4 (2003): 872-882.
+</ul>