Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP r.in.pdal - PDAL-based version of r.in.lidar #61

Closed
wants to merge 7 commits into from

Conversation

wenzeslaus
Copy link
Member

@wenzeslaus wenzeslaus commented Jul 24, 2019

r.in.pdal has most of the features from r.in.lidar and some from v.in.pdal. It uses PDAL streaming API as v.in.pdal and it uses segment lib to access the output raster, so it does not have the custom row cache and repeated looping over input point cloud (which was concern for combination of large output raster and large point cloud).

This is work in progress, but I'm not sure how far should this PR go and what to leave for the future. Features are there. CRS handling error, handling, and metadata printing is not finished. Tests are missing. Testing in general is welcome and needed. I was able to test only the basic functionality so far.

TODO:

  • a89c8e5 needs to be applied here.

@ninsbl
Copy link
Member

ninsbl commented Jan 14, 2020

liblas is deprecated on Ubuntu: https://trac.osgeo.org/grass/ticket/4020
So this enhancement can close an important gap on that platform......

@landam
Copy link
Member

landam commented May 10, 2020

@wenzeslaus What is missing to merge this PR to master? It would be nice to have r.in.pdal available for testing...

@landam
Copy link
Member

landam commented May 10, 2020

Related issue: What about moving liblas-based modules to addons? Than r|v.in.pdal could be renamed to r|v.in.lidar in master (addons to r|v.in.lidar2, r|v.in.lidar.old or so). What do you think?

@landam landam modified the milestones: 7.10, 8.0.0/7.10 May 10, 2020
@cmbarton
Copy link
Contributor

Just a note that v.in.pdal does not show up in the menu. Both r.in.pdal and v.in.pdal should be in the import menus at least.

@wenzeslaus
Copy link
Member Author

wenzeslaus commented May 14, 2020

The branch is now rebased and the code is available for testing:

git checkout -b wenzeslaus-pdal-binning master
git pull https://github.com/wenzeslaus/grass.git pdal-binning

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.
@neteler
Copy link
Member

neteler commented Oct 4, 2020

@wenzeslaus thanks for having written this new module.

I am currently testing it with the "open.NRW LiDAR" data (all LiDAR data of the federal state Nordrhein-Westfalen/DE are freely available, 2.2 TB of LAZ files).

While compiling, I got these warnings (r.in.pdal test results to be reported later):

c++    -I/home/mneteler/software/grass_master/dist.x86_64-pc-linux-gnu/include -I/home/mneteler/software/grass_master/dist.x86_64-pc-linux-gnu/include   -I/usr/include/gdal -I/usr/include -O2 -g -pipe -Wall -Werror=format-security -Wp,-D_FORTIFY_SOURCE=2 -Wp,-D_GLIBCXX_ASSERTIONS -fexceptions -fstack-protector-strong -grecord-gcc-switches -specs=/usr/lib/rpm/redhat/redhat-hardened-cc1 -specs=/usr/lib/rpm/redhat/redhat-annobin-cc1 -m64 -mtune=generic -fasynchronous-unwind-tables -fstack-clash-protection -fcf-protection -std=c++11 -DPACKAGE=\""grassmods"\"  -I/usr/include/pgsql  -I/usr/include -I/usr/include/gdal -I/usr/include/libxml2 -I/usr/include -I/home/mneteler/software/grass_master/dist.x86_64-pc-linux-gnu/include -I/home/mneteler/software/grass_master/dist.x86_64-pc-linux-gnu/include -DRELDIR=\"raster/r.in.pdal\" -o OBJ.x86_64-pc-linux-gnu/grasslidarfilter.o -c grasslidarfilter.cpp
In file included from grasslidarfilter.cpp:19:
grasslidarfilter.h: In constructor ‘GrassLidarFilter::GrassLidarFilter()’:
grasslidarfilter.h:145:18: warning: ‘GrassLidarFilter::n_class_filtered_’ will be initialized after [-Wreorder]
  145 |     gpoint_count n_class_filtered_;
      |                  ^~~~~~~~~~~~~~~~~
grasslidarfilter.h:128:10: warning:   ‘bool GrassLidarFilter::use_spatial_filter_’ [-Wreorder]
  128 |     bool use_spatial_filter_;
      |          ^~~~~~~~~~~~~~~~~~~
grasslidarfilter.h:40:5: warning:   when initialized here [-Wreorder]
   40 |     GrassLidarFilter() :
      |     ^~~~~~~~~~~~~~~~
grasslidarfilter.h:152:10: warning: ‘GrassLidarFilter::use_return_filter_’ will be initialized after [-Wreorder]
  152 |     bool use_return_filter_;
      |          ^~~~~~~~~~~~~~~~~~
grasslidarfilter.h:150:10: warning:   ‘bool GrassLidarFilter::use_class_filter_’ [-Wreorder]
  150 |     bool use_class_filter_;
      |          ^~~~~~~~~~~~~~~~~
grasslidarfilter.h:40:5: warning:   when initialized here [-Wreorder]
   40 |     GrassLidarFilter() :
      |     ^~~~~~~~~~~~~~~~
grasslidarfilter.h:155:14: warning: ‘GrassLidarFilter::base_segment_’ will be initialized after [-Wreorder]
  155 |     SEGMENT *base_segment_;
      |              ^~~~~~~~~~~~~
grasslidarfilter.h:147:12: warning:   ‘double GrassLidarFilter::zscale_’ [-Wreorder]
  147 |     double zscale_;
      |            ^~~~~~~
grasslidarfilter.h:40:5: warning:   when initialized here [-Wreorder]
   40 |     GrassLidarFilter() :
      |     ^~~~~~~~~~~~~~~~
grasslidarfilter.h:148:12: warning: ‘GrassLidarFilter::iscale_’ will be initialized after [-Wreorder]
  148 |     double iscale_;
      |            ^~~~~~~
grasslidarfilter.h:125:25: warning:   ‘pdal::Dimension::Id GrassLidarFilter::dim_to_use_as_z_’ [-Wreorder]
  125 |     pdal::Dimension::Id dim_to_use_as_z_;
      |                         ^~~~~~~~~~~~~~~~
grasslidarfilter.h:40:5: warning:   when initialized here [-Wreorder]
   40 |     GrassLidarFilter() :
      |     ^~~~~~~~~~~~~~~~
c++    -I/home/mneteler/software/grass_master/dist.x86_64-pc-linux-gnu/include -I/home/mneteler/software/grass_master/dist.x86_64-pc-linux-gnu/include   -I/usr/include/gdal -I/usr/include -O2 -g -pipe -Wall -Werror=format-security -Wp,-D_FORTIFY_SOURCE=2 -Wp,-D_GLIBCXX_ASSERTIONS -fexceptions -fstack-protector-strong -grecord-gcc-switches -specs=/usr/lib/rpm/redhat/redhat-hardened-cc1 -specs=/usr/lib/rpm/redhat/redhat-annobin-cc1 -m64 -mtune=generic -fasynchronous-unwind-tables -fstack-clash-protection -fcf-protection -std=c++11 -DPACKAGE=\""grassmods"\"  -I/usr/include/pgsql  -I/usr/include -I/usr/include/gdal -I/usr/include/libxml2 -I/usr/include -I/home/mneteler/software/grass_master/dist.x86_64-pc-linux-gnu/include -I/home/mneteler/software/grass_master/dist.x86_64-pc-linux-gnu/include -DRELDIR=\"raster/r.in.pdal\" -o OBJ.x86_64-pc-linux-gnu/main.o -c main.cpp
In file included from main.cpp:38:
grasslidarfilter.h: In constructor ‘GrassLidarFilter::GrassLidarFilter()’:
grasslidarfilter.h:145:18: warning: ‘GrassLidarFilter::n_class_filtered_’ will be initialized after [-Wreorder]
  145 |     gpoint_count n_class_filtered_;
      |                  ^~~~~~~~~~~~~~~~~
grasslidarfilter.h:128:10: warning:   ‘bool GrassLidarFilter::use_spatial_filter_’ [-Wreorder]
  128 |     bool use_spatial_filter_;
      |          ^~~~~~~~~~~~~~~~~~~
grasslidarfilter.h:40:5: warning:   when initialized here [-Wreorder]
   40 |     GrassLidarFilter() :
      |     ^~~~~~~~~~~~~~~~
grasslidarfilter.h:152:10: warning: ‘GrassLidarFilter::use_return_filter_’ will be initialized after [-Wreorder]
  152 |     bool use_return_filter_;
      |          ^~~~~~~~~~~~~~~~~~
grasslidarfilter.h:150:10: warning:   ‘bool GrassLidarFilter::use_class_filter_’ [-Wreorder]
  150 |     bool use_class_filter_;
      |          ^~~~~~~~~~~~~~~~~
grasslidarfilter.h:40:5: warning:   when initialized here [-Wreorder]
   40 |     GrassLidarFilter() :
      |     ^~~~~~~~~~~~~~~~
grasslidarfilter.h:155:14: warning: ‘GrassLidarFilter::base_segment_’ will be initialized after [-Wreorder]
  155 |     SEGMENT *base_segment_;
      |              ^~~~~~~~~~~~~
grasslidarfilter.h:147:12: warning:   ‘double GrassLidarFilter::zscale_’ [-Wreorder]
  147 |     double zscale_;
      |            ^~~~~~~
grasslidarfilter.h:40:5: warning:   when initialized here [-Wreorder]
   40 |     GrassLidarFilter() :
      |     ^~~~~~~~~~~~~~~~
grasslidarfilter.h:148:12: warning: ‘GrassLidarFilter::iscale_’ will be initialized after [-Wreorder]
  148 |     double iscale_;
      |            ^~~~~~~
grasslidarfilter.h:125:25: warning:   ‘pdal::Dimension::Id GrassLidarFilter::dim_to_use_as_z_’ [-Wreorder]
  125 |     pdal::Dimension::Id dim_to_use_as_z_;
      |                         ^~~~~~~~~~~~~~~~
grasslidarfilter.h:40:5: warning:   when initialized here [-Wreorder]
   40 |     GrassLidarFilter() :
      |     ^~~~~~~~~~~~~~~~
main.cpp: In function ‘int main(int, char**)’:
main.cpp:565:28: warning: comparison of integer expressions of different signedness: ‘unsigned int’ and ‘int’ [-Wsign-compare]
  565 |     for (unsigned i = 0; i < infiles.num_items; i++) {
      |                          ~~^~~~~~~~~~~~~~~~~~~
main.cpp:56:17: warning: unused variable ‘base_raster’ [-Wunused-variable]
   56 |     int out_fd, base_raster;
      |                 ^~~~~~~~~~~
main.cpp:58:9: warning: unused variable ‘percent’ [-Wunused-variable]
   58 |     int percent;
      |         ^~~~~~~
main.cpp:70:15: warning: unused variable ‘last_rows’ [-Wunused-variable]
   70 |     int rows, last_rows, cols;  /* scan box size */
      |               ^~~~~~~~~
main.cpp:79:22: warning: unused variable ‘cellhd’ [-Wunused-variable]
   79 |     struct Cell_head cellhd, loc_wind;
      |                      ^~~~~~
main.cpp:430:12: warning: variable ‘zscale’ set but not used [-Wunused-but-set-variable]
  430 |     double zscale = 1.0;
      |            ^~~~~~
main.cpp:431:12: warning: variable ‘iscale’ set but not used [-Wunused-but-set-variable]
  431 |     double iscale = 1.0;
      |            ^~~~~~
main.cpp:770:21: warning: ‘infile’ may be used uninitialized in this function [-Wmaybe-uninitialized]
  770 |     Rast_set_history(&history, HIST_DATSRC_1, infile);
      |     ~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

@neteler
Copy link
Member

neteler commented Oct 4, 2020

I checked the parameters and flags, comparing them to the (Python) addon r.in.pdal which we use a lot:

The Python version of r.in.pdal additionally offers the following extra functionality which would be nice to see also in this new C implementation (which offers a lot! of new functionality):

Potentially wanted flags:

-s Scan data file for extent then exit
-g In scan mode, print using shell script style

... useful for getting an idea about the input LAS/LAZ file without import. Esp. nice to be able to set the computational region properly. However, perhaps it is solved differently in this new C implementation, I didn't check that yet.

Desired parameters:

  • raster_reference: Raster map to be used as pixel geometry reference
  • raster_file: External raster map to be used as pixel geometry reference

... useful for taking a GRASS raster maps (e.g. orthophoto) or an external file (e.g. DEM GeoTIFF) as the pixel geometry reference for the computational region.

  • footprint: Footprint of the data as vector map

... useful for generating a vector polygon (at least in my area the LAS/LAZ files are not delivered in square tiles but of rather arbitrary shape.

@wenzeslaus
Copy link
Member Author

Potentially wanted flags:

-s Scan data file for extent then exit
-g In scan mode, print using shell script style

These are one of the reason this is not finished; implementing these nicely requires some additional thought. But, in other words, I plan them. After GSoC, this PR coming up in my todo list, so there is hope.

One issue there is metadata if present versus going through all the points. Going through all the point is how it is in r.in.lidar/r.in.xyz, any bad experience with that?

Desired parameters:

  • raster_reference: Raster map to be used as pixel geometry reference
  • raster_file: External raster map to be used as pixel geometry reference

... useful for taking a GRASS raster maps (e.g. orthophoto) or an external file (e.g. DEM GeoTIFF) as the pixel geometry reference for the computational region.

I'm sure you have heard about computational region, right? :-) If you want that for more than one module, the question is why not for all?

  • footprint: Footprint of the data as vector map

... useful for generating a vector polygon (at least in my area the LAS/LAZ files are not delivered in square tiles but of rather arbitrary shape.

You mean as an output, right? v.in.lidar has that as an input (a mask/filter). Footprint as a vector map would be a vector polygon, no? How would a vector polygon match the raster cells? This sounds like a convex or concave hull. Is that what you mean? I don't know, but PDAL lib may require another run through the cloud to get that. GRASS has these algs, but only in modules, not lib, I think.

@neteler
Copy link
Member

neteler commented Oct 5, 2020

One issue there is metadata if present versus going through all the points. Going through all the point is how it is in r.in.lidar/r.in.xyz, any bad experience with that?

Not a bad experience but certainly not fast. However, in r.in.lidar/r.in.xyz the PDAL API isn't really used and hence it might be all faster in this "pure" PDAL implementation.

Desired parameters:

  • raster_reference: Raster map to be used as pixel geometry reference
  • raster_file: External raster map to be used as pixel geometry reference

... useful for taking a GRASS raster maps (e.g. orthophoto) or an external file (e.g. DEM GeoTIFF) as the pixel geometry reference for the computational region.

I'm sure you have heard about computational region, right? :-)

I think so, see my line above :-)

If you want that for more than one module, the question is why not for all?

Very easy to answer: not for all but only here as we do raster binning here which often depends on further (external) data. For example, we often use orthophotos along with LAS files. The pixel geometry (esp. position) needs to be aligned with the orthophotos in case you want to perform a classification on top of it, using both orthophoto and LiDAR data. However, this may happen in different parts of the data storage infrastructure and having the possibility to "grab" the pixel positioning from an external file (orthophoto, existing DEM) comes very handy.

  • footprint: Footprint of the data as vector map

... useful for generating a vector polygon (at least in my area the LAS/LAZ files are not delivered in square tiles but of rather arbitrary shape.

You mean as an output, right?

Yes. The footprint as a vector map would be a vector polygon. Example:

image

How would a vector polygon match the raster cells? This sounds like a convex or concave hull. Is that what you mean? I don't know, but PDAL lib may require another run through the cloud to get that. GRASS has these algs, but only in modules, not lib, I think.

This is the function, based on pdal info --boundary, then conversion of its output to a GRASS GIS vector map:

https://github.com/OSGeo/grass-addons/blob/5eee668aac4a604380dcd778a0e79416384c6253/grass7/raster/r.in.pdal/r.in.pdal.py#L168

@neteler
Copy link
Member

neteler commented Oct 8, 2020

Initial functionality test:

I was able to successfully generate a DSM from a LAZ point cloud:

# area: part of Bonn, Germany (UN Campus), EPSG:25832 (close to FOSS4G 2016 location)

# download LAZ file
wget https://www.opengeodata.nrw.de/produkte/geobasis/hm/3dm_l_las/3dm_l_las/3dm_32_367_5620_1_nw.laz
# 
# import with filters
# - first return
# - classes
#   - 2: ground
#   - 9: water
#   - 17: bridge deck
#   - 20: last return not soil
#
r.in.pdal input=3dm_32_367_5620_1_nw.laz output=3dm_32_367_5620_1_nw_FIRST_RETURN \
     return_filter=first class_filter=2,9,17,20  -o

The result looks cool!

image

One issue (which I don't know how to address): water surface (river Rhine here) in open.NRW is set to the class 41 which, however, comes with ReturnNumber=0 (synthetic). At time "first, mid, last" returns are supported but not yet "synthetic".

@marisn
Copy link
Contributor

marisn commented Oct 26, 2020

Please take a look at my request to remove ground filtering options from *.in.pdal at #1049 Reason – PDAL has already changed. This functionality could cause us to constantly change module parameters to catch-up with latest developments in PDAL. Support of pipelines could be a way out, but do we need to turn import module into a full filtering module?

@marisn
Copy link
Contributor

marisn commented Nov 6, 2020

Related issue: What about moving liblas-based modules to addons? Than r|v.in.pdal could be renamed to r|v.in.lidar in master (addons to r|v.in.lidar2, r|v.in.lidar.old or so). What do you think?

PDAL can work with other formats too and thus I would vote to keep it as r.in.pdal similarly to r.in.gdal. Removing existing liblas based modules for G8 seems like a good idea. I would even not move them to addons as anyone using liblas is on an old system (with old GRASS) anyway.

@marisn
Copy link
Contributor

marisn commented Nov 6, 2020

One issue (which I don't know how to address): water surface (river Rhine here) in open.NRW is set to the class 41 which, however, comes with ReturnNumber=0 (synthetic). At time "first, mid, last" returns are supported but not yet "synthetic".

Markus, please clarify your expectations. A separate filter to include/exclude points with Withheld/Keypoint/Synthetic flags set? A single point can have more than one flag set – how conflicts should be handled?

PS. For testing you can use https://github.com/marisn/grass/tree/wenzeslaus-pdal-binning as it contains some extra functionality.

@neteler
Copy link
Member

neteler commented Nov 17, 2020

One issue (which I don't know how to address): water surface (river Rhine here) in open.NRW is set to the class 41 which, however, comes with ReturnNumber=0 (synthetic). At time "first, mid, last" returns are supported but not yet "synthetic".

Markus, please clarify your expectations. A separate filter to include/exclude points with Withheld/Keypoint/Synthetic flags set? A single point can have more than one flag set – how conflicts should be handled?

No strong opinion here, perhaps an include_classes parameter and/or exclude_classes parameter?
For sure I am more in favour of seeing this merged rather than overcomplicating it.

PS. For testing you can use https://github.com/marisn/grass/tree/wenzeslaus-pdal-binning as it contains some extra functionality.

Cool, I see new "binning" methods in https://github.com/marisn/grass/tree/wenzeslaus-pdal-binning/raster/r.in.pdal
This is on top of this PR?

@marisn
Copy link
Contributor

marisn commented Nov 18, 2020

No strong opinion here, perhaps an include_classes parameter and/or exclude_classes parameter?
For sure I am more in favour of seeing this merged rather than overcomplicating it.

We have to thing more about how it should be presented to the user as class filter is already implemented and this is some kind of flag based filter.

Cool, I see new "binning" methods in https://github.com/marisn/grass/tree/wenzeslaus-pdal-binning/raster/r.in.pdal
This is on top of this PR?

Unfortunately I am just working on things I miss from the module and not features missing for final inclusion. Still I start to worry that soon my version will have diverged (from this PR) too much for an easy merge and will have to replace version of this PR.

@marisn marisn moved this from To do to In progress in OSGeo Virtual Community Sprint 2020 Nov 19, 2020
@marisn
Copy link
Contributor

marisn commented Dec 14, 2020

This is the function, based on pdal info --boundary, then conversion of its output to a GRASS GIS vector map:

https://github.com/OSGeo/grass-addons/blob/5eee668aac4a604380dcd778a0e79416384c6253/grass7/raster/r.in.pdal/r.in.pdal.py#L168

And thus it should go into an add-on. Markus, please open a feature request for such add-on. I could see this as an option for v.in.pdal, but certainly not in r.in.pdal.

@neteler
Copy link
Member

neteler commented Dec 15, 2020

This is the function, based on pdal info --boundary, then conversion of its output to a GRASS GIS vector map:
https://github.com/OSGeo/grass-addons/blob/5eee668aac4a604380dcd778a0e79416384c6253/grass7/raster/r.in.pdal/r.in.pdal.py#L168

And thus it should go into an add-on. Markus, please open a feature request for such add-on. I could see this as an option for v.in.pdal, but certainly not in r.in.pdal.

Here the use case: imagine a bunch of LAS/LAZ files (like many thousand) which partially overlap and are not rectangular as it happens with those provided for our federal state here (and no footprint map provided). With the possibility to write out the footprints of these files, you generate a tile index vector map while importing.

@marisn
Copy link
Contributor

marisn commented Dec 16, 2020

Here the use case: imagine a bunch of LAS/LAZ files (like many thousand) which partially overlap and are not rectangular as it happens with those provided for our federal state here (and no footprint map provided). With the possibility to write out the footprints of these files, you generate a tile index vector map while importing.

Markus, there are a lot of useful features that could be implemented in r.in.pdal but the main question remains – should they. r.in.pdal is already quite feature heavy.
As of footprints – there are two options:

  • use output of pdal boundary. Still if user filters out points during import, imported data could occupy smaller area.
  • create convex hull based on matching points. Someone has to come up with an algorithm of growing convex hull on-fly.

And still I think it would be best to have v.pdal.tileset module to implement such functionality to not make r.in.pdal even larger monster. KISS.

@lucadelu
Copy link
Contributor

lucadelu commented Jul 8, 2021

Could we close this since it is obsolete?

@wenzeslaus
Copy link
Member Author

Replaced by #1200.

@wenzeslaus wenzeslaus closed this Jul 14, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
No open projects
Development

Successfully merging this pull request may close these issues.

None yet

7 participants