When you have a big pile of LiDAR data, sometimes things are a bit harder to get to the result you need. You can convert the LiDAR data to for instance DEM/DTM and DSM high resolution rasters (see LAS Dataset To Raster—Help | ArcGIS for Desktop  ), but any product you derive from these rasters will yield another large raster. If you are not directly interested in having a physical result you may want to consider using raster chain functions (see What are the functions used by a raster or mosaic dataset?—Help | ArcGIS for Desktop  there are lot to choose from).

Raster chain functions can be very useful for deriving products from multispectral imagery, but they are also very useful for all kinds of raster analysis. Take the case mentioned in a previous post Creating a 3D tree with Arcpy for 3D Analysis  where I created 3D geometries to be able to determine the distance between the tree crown and the nearest overhead power line. What if we would like to apply this in a raster environment using DEM/DTM and DSM information and a raster representation of the “surface” composed by the lowest overhead power lines?

The first problem may be Translating the lowest overhead power lines into a continuous surface . 

So, let’s assume that we have three raster datasets for the same area.

• Digital Elevation Model / Digital Terrain Model
• Digital Surface Model
• Cable Surface

We can create composite bands raster containing the 3 rasters (per tile), where band 1 is the DEM, band 2 is the DSM and band 3 is the Cable surface.

Below a snippet of Python code that I used to generate the composite bands:

 def main():
    import arcpy
    import os
    arcpy.env.overwriteOutput = True

    ruta_tif = r'D:\Xander\_LiDAR\ElSalto'
    lst_prefix = ['DTM_AMALFI_', 'DSM_AMALFI_']
    ext = '.tif'
    ws = r'D:\Xander\_LiDAR\ElSalto\tmp.gdb'

    for b in range(2, 18):
        dtm = os.path.join(ruta_tif, '{0}{1}{2}'.format(lst_prefix[0], b, ext))
        dsm = os.path.join(ruta_tif, '{0}{1}{2}'.format(lst_prefix[1], b, ext))

        in_ras = '{0};{1}'.format(dtm, dsm)
        out_ras_name = 'comp_bands_v%02d' % (b,)
        out_ras = os.path.join(ws, out_ras_name)

        # composite bands [2-17]
        arcpy.CompositeBands_management(in_rasters=in_ras, out_raster=out_ras)

if __name__ == '__main__':
    main()

The next step is to create a Mosaic dataset with three bands and 32_BIT_FLOAT as Pixel Type (if this is the pixel type of the input models). To create and update Mosaic Datasets you will need a Standard license level or higher:

@

Once you have the empty Mosaic dataset created, you can load the composite bands into the Mosaic:

… and build the overviews:

Now we have a fast performing mosaic, but when we add it to the TOC, ArcMap will create a color composite of the 3 bands and that doesn’t make much sense. This is where the Raster Chain Functions come in.

Suppose that we want to use a visualization of the Digital Terrain Model using a shaded relief as default display, we can configure this in the Mosaic.

In the Catalog window, double click on the Mosaic and activate the “Processing Templates” tab. In this case I have already configured a number of processing templates and the default processing template is set to a Shaded Relief of the elevation:

If I add the mosaic dataset to the TOC, I will display very fast and show this:

If I look at the properties of the Mosaic it reveals the size of the data contained:

This is still a small example, but you could have TB of data underneath and the display would still be fast. How does that work? The Mosaic dataset has all kinds of intelligence to yield fast rendering. When you apply a processing template with 1 or more Raster Chain Functions these functions are not applied to the entire dataset, but only to those pixels displayed. You can calculate all kinds or products, without creating large datasets for each products, but you still have the possibility to work with them as physical rasters. You can publish them as service and apply the raster chain functions to the service, to provide the users with different kinds of visualizations in the browser. The mosaic dataset and even the service when added to the TOC will be recognized by ArcMap as valid raster layer that can be used in the Raster Calculator and other tools.

So, how can we configure the processing templates to extract the information we need? Let’s look at a very simple example. Suppose, we would like to be able to extract only the Digital Terrain Model and work with that. This can be done using the “Extract Band Function”:

And configure that it should extract the first band which in our case corresponds to the Digital Terrain Model:

… and in the General tab configure the output pixel type (in our case it is 32 bits float).

You can do all kinds of stuff with these functions like visualizing the flow of water over the surface with arrows, using vector field renderer:

Which will look like this:

Let’s get back the topic of this post. How can we use this to validate the design of an overhead power line? For that I would need to know the distance between the cable and the Digital Surface Model. You extract band 3 (the Cable surface) and band 2 (the DSM), execute a Minus (Spatial Analyst required) and apply the Statistics and Histogram Function to enhance the output result:

This will look like this:

If we compare this result with the Canopy Height Model (CHM) or normalized DSM (nDSM), you will see that not all trees are near to the overhead power lines:

So, using a multiband Mosaic datasets with the DTM, DSM and cable surface enables us to publish this analytical functionality to an image service and share this in a viewer (like the Web AppBuilder) with users that don’t have access to Desktop software.

Please note that you will need Spatial Analyst activated in Desktop to load a processing template that uses functions that require Spatial Analyst.

Spatial Analyst Lidar and other point clouds Imagery and Remote Sensing Elevation Data python snippets