How can I compute a raster showing the number of return points per square meter using SAGA GIS or other open source software?
My goal is to obtain a mask which would show me the places (cells) where there were no hits and would help me in interpreting LiDAR-derived data.
I develop a free and open-source GIS called Whitebox Geospatial Analysis Tools (can be downloaded here) that has extensive analysis functionality for processing LiDAR data. Whitebox contains a tool specifically for calculating the point-density of LiDAR LAS files called Point Density LiDAR.
The tool is highly specific to LiDAR, taking one or more LAS files as input. It can also calculate the density of first-return points only, last-return points only or all-returns. Additionally, if your LiDAR data contains classification information, you can optionally exclude points based on their classification type (e.g. buildings or vegetation points).
There is a second tool, called Hex-binning, which will also perform a density-type analysis on LAS files (Point-type ShapeFiles as well), but will output a vector (ShapeFile) hexagonal grid. This can be a very effective way of visualizing high-density point coverage (see blog on hex-binning here).
You might like to use Fusion. It is a free software for point cloud processing and visualization .
1- Look into the manual for a command program called Catalog. It returns descriptive statistics from the point cloud.
What you want is the Catalog's switch density:area,min,max. The manual description says:
Creates an image for all data files that shows the return density for the area represented by each pixel. area is the pixel area, min is the minimum acceptable point density per unit area, and max is the upper limit for the acceptable density range. Cells with point densities falling within the min-max range are colored green, cells with point densities below the minimum are colored red, and cells with densities above the maximum are colored blue.
You need to set min and max accordingly so it will highlight the information you want. For example: if you want to see pixels where there were less than 1 return, set min equal to 1. Pixels with less than one return per square meter will show up red.
The command line would be similar to the following, considering pixels of 5 meters (25 m²).
Catalog /density:25,1,5 pointcloud.las
If you need more help to run this command line, here I provided a detailed example about it.
This is the catalog /density output, taken from one sample data provided by Fusion:
The RGB picture from the scene above can be seen here.
2- There is another option, more flexible than Catalog.
Use the command GridMetrics to extract the number of returns per pixel. Set the pixel size the way is more suitable for analysis. The output will be a comma separated value (csv) file where each line will correspond to a pixel.
Example (with pixel size equal to 5):
gridmetrics /noground /nointensity 2 5 4800K_gridmetrics lda_4800K_data.las
/nogroundis a switch used together with normalized datasets, /nointensity is a switch for not generating LiDAR intensity statistics (will speed up processing), and 2 is the height break for calculating cover statistics.
Then, use the command CSV2Grid to transform the .csv file into an ASCII raster format file (.asc), and you will be able to open it in different GIS software.
Example (where the variable of interest 'total number of returns' was in column 5 from previous gridmetrics output):
csv2grid 4800K_gridmetrics_all_returns_elevation_stats.csv 5 4800K_return_per_pixel.asc
To check the number of points per grid cell from a LIDAR point cloud in GRASS you can use the r.in.xyz module. This module creates a grid from the point cloud using a "method" parameter by which you choose how to aggregate the points when importing. If you choose method=n then the resulting raster will contain a count of points for each cell. By setting method=mean you get the average value of points in each grid cell.
More details on the GRASS Wiki LIDAR processing page.
