1

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
  • As an alternative approach you can create a raster from non ground LAS points. Then you can use raster calculator to create cross sections. I hope your LIDAR data was taken when the leaves on trees. – Amadeus Feb 15 '18 at 21:36
  • Linking related question: gis.stackexchange.com/questions/183924/… – Andre Silva Jan 17 at 13:18
3

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
  • Hey Jae! Too funny! Looking at this, I see a lot of things I am not familiar with, but DensityMetrics looks to be exactly what I need. I am trying to follow the manual but am stumped at how to make a 'DTM'. I have only ever dealt with GRID format. – D_C Feb 15 '18 at 4:20
  • 1
    I’ve edited my original answer with detail of each step. – Jae Feb 15 '18 at 12:47
  • do you know why the lowest cell size and height slice I can go using densitymetrics is 1 and 1? I am trying to get 0.5 and 0.5 but keep getting memory errors. I have tried busting up the LAS file with polyclipdata by 4, 9, and 16 squares but still no go. – D_C Feb 27 '18 at 5:26
1
  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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.