Processing Lidar with standard ArcPy

Photo by Christopher Burns on Unsplash

I’m involved in a project that identifies the location of fire hydrants with Light Detection and Ranging (lidar) data. ArcGIS Pro has a special toolbar ribbon dedicated to lidar visualization and analysis. I have learned some new techniques that are necessary to incorporate those tools in a standalone python script.

Store output in memory

I’ve been using arcpy’s memory-based workspace since the Python 2.7 days when it was called in_memory/. It stores the output of python methods in the computer’s RAM instead of on disk. I love that it keeps my geodatabases and file folders free of “intermediate data” and usually speeds up processing.

[embed]Write geoprocessing output to memory-ArcGIS Pro | Documentation ArcGIS provides a system memory-based workspace where output feature classes, tables, and raster datasets can be…pro.arcgis.com

Previously, I would delete the data from memory at the end of my script. But now, I have learned that I should instead delete it in the finally clause of a try-catch block that I place around all of the code. That way, the memory will be cleared at the end, even if the script crashes. I just need to add a little extra code to make sure that if the script crashes early in the execution, it doesn’t try to delete data it hasn’t created yet (resulting in another crash!). It works like this:

Starting_Point_Layer = None
try:
  Starting_Point_Layer = arcpy.management.MakeFeatureLayer
    (in_features="path_on_disk", out_layer="memory/Starting_Point_Layer")
catch: 
  //error processing
finally:
  if Starting_Point_Layer is not None:
    arcpy.Delete_management(in_data=Starting_Point_Layer)

I have also learned there are some methods that will not accept a memory/ path as an output and/or an input. These include:

• arcpy.management.Merge()
• arcpy.gapro.DissolveBoundaries()
• arcpy.analysis.Erase()

and a few others that I will talk about later.

Export Mosaic Dataset Geometry

Our lidar data is sorted into tiles and stored in a “Mosaic Dataset”.

[embed]Mosaic datasets-ArcGIS Pro | Documentation Mosaic datasets are used to manage, display, serve, and share raster data.pro.arcgis.com

When I add the mosaic dataset to a map, it looks like a grid. Each tile has a border called a “Footprint.” We wanted to run the algorithm that we had developed in one tile at a time. The first step in our workflow is to select all of the fire hydrant points that fall within a footprint.

I knew I needed to convert the footprints into a feature layer before I could select the hydrants that fall within them by location. But to my surprise, the standard arcpy.management.MakeFeatureLayer() command that I’d used many times before did not work. These footprints look like they are normal polygons, but they are not!

I was finally successful when I used arcpy.management.ExportMosaicDatasetGeometry() with an input parameter of geometry_type="FOOTPRINT". Just like MakeFeatureLayer(), it has a where_clause parameter, so I could export each footprint individually by name.

Mosaic_Footprint = arcpy.management.ExportMosaicDatasetGeometry
  (in_mosaic_dataset=DEMPath, out_feature_class="memory\Mosaic_Footprint",
  where_clause="Name = '" + DEMName + "'", geometry_type="FOOTPRINT")

Now the footprint will act like an actual GIS feature that I can do typical GIS things with.

arcpy.management.SelectLayerByLocation
  (in_layer=Hydrant_Layer, overlap_type="WITHIN_CLEMENTINI",
  select_features=Mosaic_Footprint, search_distance="",
  selection_type="NEW_SELECTION", invert_spatial_relationship="NOT_INVERT")

Count the number of features

Our algorithm does repeated selecting and clipping in order to narrow down our results. When running geoprocessing tools or using the python window in ArcGIS Pro, it’s easy to see what’s happening at each step along the way, because selected features are highlighted and clipped-out shapes disappear before my eyes. In a python IDE, all I have is text on the screen.

If I’m doing something wrong, it doesn’t always throw an exception with an error message. Sometimes it’ll just output a lot more or less features than I’m expecting. (Usually zero features or all the features!). So I can catch this early, I’ve resorted to using the GetCount() method on my feature layer after almost every command.

Number_Of_Hydrants = int(arcpy.management.GetCount(Hydrant_Layer).getOutput(0))

if Number_Of_Hydrants == 0:
  print("0 hydrants, exiting script.")
  sys.exit()
else:
  print(str(Number_Of_Hydrants) + " hydrants fall within the footprint.")

When I print messages to the console, I can follow along. It’s fun to watch long-running scripts with a lot of steps go!

Extract a clipped and filtered LAS file

Each tile in the mosaic dataset references a .las file. LAS is the standardized binary file format for storing 3D point cloud data. Our LAS data is classified by landcover type. ArcGIS Pro has tools that let you symbolize and filter by class.

We assumed that hydrants would have been classified as “1 — Unassigned” or “2 — Ground”, so we restricted our analysis to only those classes. Filtered lidar can be saved as a LAS Dataset (.lasd file)

[embed]What is a LAS dataset? A LAS dataset stores reference to one or more LAS files on disk, as well as to feature classes containing surface…desktop.arcgis.com

Next, we wanted to compare the elevation of each point with the elevation of the points around it, to find places where a single point stuck up 2–3 feet above the rest. I thought I’d buffer each point, then loop through all the points in the buffer.

It turns out, just like a mosaic dataset, LAS files aren’t made up of standard GIS point features. If I want ArcPy tools to be able to see their coordinates and attributes, I have to convert them to something else. This can be done with the LASToMultipoint() command, which will give the points a “SHAPE” attribute, followed by the AddZInformation() command, which will give them a “Z” elevation attribute.

Unfortunately, there’s a catch. A LAS Dataset cannot be used as the input to the LASToMultipoint() command. It only accepts .las files. After a lot of trial and error, I settled on this workflow as a solution:

1. Create a LAS Dataset from the original unfiltered .las file. arcpy.management.CreateLasDataset(input="original.las", out_las_dataset="unfiltered.lasd")
2. Create a LAS Dataset Layer in order to apply the class filter. LAS_Dataset_Layer = arcpy.management.MakeLasDatasetLayer(in_las_dataset="unfiltered.lasd", out_layer="memory/LAS_Dataset_Layer", class_code="1;2")
3. Use the Extract LAS tool to save the filtered layer back out to a .las file. I clipped the lidar by a polygon boundary with this tool as well. arcpy.ddd.ExtractLas(in_las_dataset=LAS_Dataset_Layer, target_folder=LASPath, extent="MINOF", boundary="path_to_polygon_clipper", process_entire_files="PROCESS_EXTENT", **name_suffix="_Extract"**, remove_vlr="MAINTAIN_VLR", rearrange_points="REARRANGE_POINTS", compute_stats="COMPUTE_STATS", **out_las_dataset=""**, compression="SAME_AS_INPUT")

The key here is to provide a name_suffix value and don’t provide an out_las_datasetvalue. This will force a .las file to be created instead of a .lasd file, and control what the file is named. Now I could run the LAStoMultipoint()tool:

arcpy.ddd.LASToMultipoint(input="original_Extract.las", 
out_feature_class="multi_point", average_point_spacing=0.0001, 
class_code=[1, 2])

At this point you might notice that I inputted aclass_code filter as a parameter directly into the LAStoMultipoint() tool. Couldn’t I have just done that with my original.las file and avoided all those other steps? The answer is yes, if I didn’t want to clip the lidar with a polygon boundary. Clipping was very important for us to narrow the results down and increase the processing speed.

Photo by Anna Savina on Unsplash

Later, I used the same las file in the LocateOutlier() tool. That tool provides another way to find lidar locations that are different from the others around them, and it made our algorithm more efficient. The LocateOutlier() tool does not accept a class_code filter as a parameter. So, I recommend steps 1–3 above as the only way I know to create a multipoint feature class and locate outliers on clipped and filtered las.

Create a short loop

The final version of my script loops through all of the outliers. There are thousands of them in a typical grid. While I was testing things, I only wanted to loop through the first few of them, so I used a range() function. But my script would keep crashing when I got to the end of the loop.

I fixed this by catching the StopIteration exception.

try:
  for iteration in range(0, 10):
  //do stuff with the outlier point
except StopIteration:
  pass

//do other stuff

Now I can continue doing things after 10 loops have finished.

Join to a Near Table

The final step of our algorithm uses a near table to determine which of the outliers are closest to water mains. I run the GenerateNearTable() method to get a NEAR_RANK for each outlier. Then I join the output back to the outliers layer.

Near_Table = arcpy.analysis.GenerateNearTable(
  in_features=Starting_Point_Layer, near_features=Outlier_Layer,
  out_table="path_to_Near_Table", search_radius="75 Feet",
  location="NO_LOCATION", angle="NO_ANGLE", closest="ALL",
  closest_count=3, method="PLANAR")

GenerateNearTable() is another one of those methods that will not allow you to store its output in memory. I have to write it to a file on disk. When I was testing things, I would often be writing the same file multiple times with small changes.

This image is from the ESRI Generate Near Table tool reference document

I have discovered that ArcPy behaves weirdly when overwriting an existing near table. Instead of removing all the old records in the table and replacing them with new ones, it creates new IN_FID and NEAR_FID fields followed by an underscore and a number (_1,_2, etc.). Then it populates those new fields with the new attributes.

Since my next step was to join based on the NEAR_FID field, the new records would not participate in the join!

arcpy.management.JoinField(in_data=Outlier_Layer, in_field="OID", 
  join_table="path_to_Near_Table", join_field="NEAR_FID", 
  fields="IN_FID;NEAR_DIST;NEAR_RANK")

This function produces an empty dataset when the join_field has spontaneously changed its name to NEAR_FID_1!

This was incredibly confusing to troubleshoot. But once I finally discovered the problem, the solution was easy. I just made sure to delete the near table from disk in the finally block at the same time I delete all the other intermediate data from memory.

This script has helped us improve the accuracy of the hydrant locations in our GIS database. The project was my first introduction to lidar, and I’ve enjoyed watching presentations about other cool things people are doing with it. I think it’s worthwhile to know how to use it. I hope my tips encourage you to try it too! You can get some data to experiment with here:

USGS LidarExplorer