r.in.lidar: use base raster with its resolution

Computational region is used by default for input and output, when flag -d is enabled,
input follows base raster resolution (and also extent, this is suboptimal).

Adding also Bash test with small artificial data for different resolutions
in combination with base raster.

This commit follows https://trac.osgeo.org/grass/changeset/66094 and r66151.


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

raster/r.in.lidar/main.c

@@ -152,6 +152,7 @@ int main(int argc, char *argv[])
     struct Option *method_opt, *base_raster_opt, *zrange_opt, *zscale_opt;
     struct Option *trim_opt, *pth_opt, *res_opt;
     struct Flag *print_flag, *scan_flag, *shell_style, *over_flag, *extents_flag, *intens_flag;
+    struct Flag *base_rast_res_flag;
 
     /* LAS */
     LASReaderH LAS_reader;
@@ -310,6 +311,11 @@ int main(int argc, char *argv[])
     intens_flag->description =
         _("Import intensity values rather than z values");
 
+    base_rast_res_flag = G_define_flag();
+    base_rast_res_flag->key = 'd';
+    base_rast_res_flag->description =
+        _("Use base raster actual resolution instead of computational region");
+
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
@@ -660,15 +666,25 @@ int main(int argc, char *argv[])
     /* using row-based chunks (used for output) when input and output
      * region matches and using segment library when they don't */
     int use_segment = 0;
-    if (base_raster_opt->answer && res_opt->answer)
+    int use_base_raster_res = 0;
+    if (base_rast_res_flag->answer)
+        use_base_raster_res = 1;
+    if (base_raster_opt->answer && (res_opt->answer || use_base_raster_res))
         use_segment = 1;
     if (base_raster_opt->answer) {
+        if (use_base_raster_res) {
+            /* TODO: how to get cellhd already stored in the open map? */
+            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 */
+        }
         /* TODO: do we need to test existence first? mapset? */
         base_raster = Rast_open_old(base_raster_opt->answer, "");
         base_raster_data_type = Rast_get_map_type(base_raster);
     }
     if (base_raster_opt->answer && use_segment) {
-        Rast_get_input_window(&input_region);  /* we have split window */
+        if (!use_base_raster_res)
+            Rast_get_input_window(&input_region);  /* we have split window */
         /* TODO: these numbers does not fit with what we promise about percentage */
         int segment_rows = 64;
         /* writing goes row by row, so we use long segments as well */
@@ -915,6 +931,7 @@ int main(int argc, char *argv[])
             }
             else if (use_segment) {
                 double base_z;
+                /* 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

raster/r.in.lidar/testsuite/data/points.txt

@@ -0,0 +1,16 @@
+21,19,30,1,2
+21,16,30,1,2
+24,18,35,2,5
+26,19,35,2,5
+23,16,35,2,5
+27,17,35,2,5
+29,19,40,1,5
+28,14,40,1,5
+21,14,32,2,5
+24,13,37,1,5
+27,14,37,1,5
+29,16,40,2,5
+21,11,32,2,5
+24,11,32,3,3
+26,11,32,3,3
+29,11,40,2,5

raster/r.in.lidar/testsuite/test_base_resolution.sh

@@ -0,0 +1,36 @@
+#!/usr/bin/env bash
+
+set -e
+set -x
+
+basename="test_rinlidar_"
+
+g.region n=20 s=10 e=30 w=20 res=2.5
+r.in.xyz data/points.txt output=${basename}base x=1 y=2 z=3 separator=comma
+
+echo "With base raster resolution matching current region"
+
+r.in.lidar input=data/points.las output=${basename}with_region \
+    base_raster=${basename}base method=min -o
+echo "Almost all in the following r.univar output should be zero"
+r.univar ${basename}with_region
+echo "Automatic test if there are only allowed values..."
+r.univar ${basename}with_region -g \
+    | grep -ve "=0$" | grep -ve "=-nan$" | grep -e "=[^2-9][^12345789]$"
+
+echo "With base raster resolution different from current region"
+
+g.region res=5
+
+r.in.lidar input=data/points.las output=${basename}with_base \
+    base_raster=${basename}base method=min -o -d
+echo "Almost all in the following r.univar output should be zero"
+r.univar ${basename}with_base
+echo "Automatic test if there are only allowed values..."
+r.univar ${basename}with_base -g \
+    | grep -ve "=0$" | grep -ve "=-nan$" | grep -e "=[^12356789]$"
+
+echo "Test successful"
+echo "When running manually maps can be now removed with:"
+echo "  g.remove type=rast pattern='test_rinlidar_*' -f"
+echo "However, the region was changed to whatever the test needed."