Getting statistics in cross sections at varying heights in LAS files?

Question

I am looking to cross section my forest at varying heights. I have high return LAS data (8 points/m2 and average point spacing of 0.2m) and I am looking to do more than the simple 1, 2, 3 approach of vegetation classification within the LAS file.

I would like to have the flexibility to create multiple cross sections at specific height increments. Looking at, e.g., LASTools and ArcGIS tools, I can't do what I would like. So, I am looking at scripting something myself but I am running into personal knowledge roadblocks.

Firstly, when I create a LAS dataset in Arc each point that is created seems to represent 3500 LiDAR returns; when I run multi-point to single point this does not change.

Secondly, when I try and tile my problem using a fishnet approach the selected points for the current processing tile are not bounded by the tile and produce patterns of selected points well outside the tile (a multi-point property I'm sure).

How do I get at the data to perform my cross sections? Is there a way to play with LAS data outside of a LAS dataset?

UPDATE:

I have successfully sliced up my LAS data using some painful arcpy methods (see code). I am in the process now of running Jae's answer and the Fusion toolbox looks like the way to go in the future.

for file in os.listdir(in_las_folder):
    if file.endswith(".las"):
        print '...Working on LAS file {} of {}'.format(current, count)
        process = 2
        las = os.path.join(in_las_folder, file)
        las_name = os.path.splitext(file)[0]
        arcpy.CreateLasDataset_management(las, out_location_las + '\\ld_{}.lasd'.format(las_name), "NO_RECURSION", "", spatial_reference)
        lasd = out_location_las + '\\ld_{}.lasd'.format(las_name)
        arcpy.ChangeLasClassCodes_3d(lasd,"3 0;4 0;5 0")

        print '.....Creating Cross Sections'
        while process <= 38:
            print '........Working on heights between {}m and {}m'.format(process - 2, process)
            raster = out_location_height + '\\h{}_{}'.format(las_name,  process)
            if overwrite_height_model == 'NO':
                if arcpy.Exists(raster):
                    process += 2
                    continue
            print '.........Classify By Height'
            if process == 2:
                value = 20
                arcpy.ClassifyLasByHeight_3d(lasd, "GROUND", "19 1;20 2;21 4;22 6;23 8;24 10;25 12;26 14;27 16;28 18;29 20;30 22;31 24", "NONE", "false", "DEFAULT", "false")
            elif process == 26:
                arcpy.ChangeLasClassCodes_3d(lasd,"19 0;20 0;21 0;22 0;23 0;24 0;25 0;26 0;27 0;28 0;29 0;30 0;31 0")
                arcpy.ClassifyLasByHeight_3d(lasd, "GROUND", "19 24;20 26;21 28;22 30;23 32;24 34;25 36;26 50", "NONE", "false", "DEFAULT", "false")
                value = 20
            arcpy.MakeLasDatasetLayer_management(lasd, 'las_layer', [value])
            print '.......Creating raster'
            arcpy.LasDatasetToRaster_conversion('las_layer', raster, "ELEVATION", "BINNING MAXIMUM NONE", "FLOAT", "CELLSIZE", "0.5", "1")

            value += 1
            process += 2
        current += 1

        sys.exit()
print 'Creating full rasters'
process = 2
current = 1
while process <= 36:
    print '...Working on {} of {}'.format(current, count)
    for ras in arcpy.ListRasters():
        if int(ras.split('_')[1]) == process:
            height.append(ras)
    arcpy.MosaicToNewRaster_management(height, out_location_height, 'h{}'.format(process), spatial_reference, '32_BIT_FLOAT', '1', '1')
    for f in height:
        arcpy.DeleteFeature_management(f)
    process += 2
    current += 1

Answer 1

If you’re looking for point density at specific height increments, you could look at the FUSION command line tool DensityMetrics

Assuming your data are in units of metres:

Step 1: Create a ground file from class 2 (ground) and class 8 (model key) LAS points, at 2 metre spatial resolution.

gridsurfacecreate /class:2,8 /minimum YourDrive:\Path\groundfile.dtm 2 M M 0 0 0 0 YourDrive:\Path\*.las

Step 2: Create point density slices at 1 metre height intervals using a 20 metre cell size.

densitymetrics YourDrive:\Path\groundfile.dtm 20 1 YourDrive:\Path\densities.csv YourDrive:\Path\*.las

Step 3: (optional) Rasterize output from the CSV file. Here we will rasterize the 5th column representing densities from the 0-1 metre slice.

csv2grid YourDrive:\Path\densities_all_returns.csv 5 return_density_0to1m.asc

Answer 2

1. Profiles: You can work with cross-sections/profiles for LAS data in other software, e.g. CloudCompare has a section tool and it's open source. Some other commercial software, like Global Mapper and Quick Terrain Modeler have some superb profile tools, if you can obtain a licence.

2. Play with LAS data: ArcGIS doesn't permit field export when creating a shapefile from a LAS file. You can extract the field info from LAS by exporting them in shapefile format (e.g. Global Mapper exports all available fields) or another format (convert to XYZ, CSV or ASCII using CloudCompare).

3. Your ~3500 returns per point conversion is obviously not correct. I guess you are using LAS to Multipoint Tool. Try to enter the correct values for Average Point Spacing and Coordinate System, these could be a source of error in your conversion.