This blog is written by Clayton Crawford, a Product Engineer in the Software Products Group’s 3D Team in Redlands.
Lidar Solutions in ArcGIS_part2: Creating raster DEMs and DSMs from large lidar point collections
Raster, or gridded, elevation models are one of the most common GIS data types. They can be used in many ways for analysis and are easily shared. Lidar provides us with the opportunity to make high quality elevation models of two distinct flavors: first return and ground. A first return surface includes tree canopy and buildings and is often referred to as a Digital Surface Model (DSM). The ground, or bare earth, contains only the topography and is frequently called a Digital Elevation Model (DEM).
ArcGIS provides tools to take large lidar point collections, and optionally other surface related information like photogrammetric breaklines, and use them to produce high quality raster surfaces.
Coming up with a plan
Before continuing some basic factors need to be evaluated:
Consideration of these factors will help determine whether you’ll be producing one raster or a collection. Part of this process requires you figure out how many rows and columns you’re willing to have in one raster. This depends on what you intend to do with the raster in terms of analysis, display, and potential sharing or distribution of the data. The desire to work off one dataset for analysis can be in conflict with practical constraints associated with dataset size. Another thing to consider is the amount of lidar you have. Trying to process 10 billion lidar points as one dataset, while possible, is likely to prove unwieldy. It’s pretty much a given you’ll be making multiple rasters from this amount of lidar, so consider splitting up the lidar processing as well. Not only does this keep individual datasets at reasonable sizes, it keeps process duration on those datasets shorter. The longer a process takes to execute, the more likely something’s going to go wrong (e.g., power outage). We all want to reduce risk, right?
If you’ve determined you need to split up your data, the next question is how. Will it be based on a regular grid system, political boundaries, or division based on an anticipated application? Since lidar collections tend to have multiple uses, splitting them up based on a regular grid system or political divisions like county boundaries makes the most sense. An engineer can mosaic the different pieces he or she needs for an individual project. If the intended use is weighted heavily for one type of application, such as hydrology, then use divisions logical for the application. For example, in the case of hydrology, watershed boundaries are a good candidate.
ArcGIS supports many raster formats, so you have a choice of what format to write to. This decision is best based on the intended use of the product. If it’s to be shared with the general public you might think about distributing in either TIFF or JPEG format. For analysis on the ArcGIS platform, the ESRI GRID format is best. This is the native format for many functions so to improve I/O efficiency, and therefore processing time, using GRID is the way to go.
The first step in getting from lidar points to raster is loading the points into a GDB. If you haven’t already done so, review the first part of this series. Use the link at the beginning of this blog to get to it.
Using the Point to Raster tool
If your only source of data is the lidar you can use the Point to Raster tool to produce a raster elevation models. While this does not produce the highest quality result possible lidar tends to be so dense that for many applications the accuracy is good enough and the convenience and speed of this tool make it worthwhile.
The Point to Raster tool produces gridded elevation models from lidar point sets.
While Point to Raster offers the easiest and fastest way to produce a raster from lidar, it has a significant drawback. You can end up with many NoData cells since it only returns values for cells that have one or more points in them. The problem is exacerbated when only using ground points to make a DEM because gaps in the data occur where there’s vegetation and buildings obscuring the ground. To reduce the salt & pepper effect of NoData vs. data cells you can increase the output cellsize relative to the average point spacing. You can also reduce the number of NoData cells after execution of Point to Raster by using this expression in Spatial Analyst calculator on the output (in this example the output from Point to Raster is called ‘point2ras):
You can run the expression multiple times to fill in larger NoData areas, but I wouldn’t recommend running it more than a couple times. It’s better to just accept larger void areas as a consequence of using this approach.
Using the Terrain Dataset
If you have photogrammetric breaklines to go along with your lidar, or need higher quality results than can be produced with the Point to Raster tool, use the terrain dataset. For an overview of the terrain dataset and related help topics check out the online help here.
On the left is a surface made without breaklines along the river banks. The image on the right has breakline enforcement. Breaklines are important for maintaining the definition of water related features in the elevation model.
Breaklines are used to capture linear discontinuities in the surface. The most common types are edge of pavement, lake shorelines, single line drains for small rivers and double line drains for large rivers. Sometimes breaklines are also collected to help define and sculpt the surface without necessarily representing discontinuities. Examples of these include contour-like form lines and the crests of rounded ridges.
Breaklines, while frequently used in bare earth models, tend to be detrimental when used with first return surfaces because they can be in conflict with the above ground points. For example, breaklines capturing road edge of pavement can be coincident in XY but different in Z to points in tree canopy overhanging the road. Because of this, consider excluding breaklines from your first return surface or at least those where you know there’s potential conflict.
The most efficient means of organizing breaklines for use in a terrain dataset – see table below- is to separate them into different feature classes based on Surface Feature Type (SFType). SFType controls how the features are enforced in the model and how the natural neighbor interpolator, used during rasterization, interprets the surface as it crosses over these features. A distinct break in slope will occur across ‘hard’ features but not across ‘soft’ features.
Common input measurement types and recommend feature class storage type and SFType settings for terrain dataset definition.
It’s best for the sake of terrain performance to place all hardlines together in one feature class. It’s understood this might not be possible, for example, if you have the need to keep road and water features separate. It’s OK, just keep in mind the fewer feature classes used to define a terrain the better.
The ‘Replace’ SFType deserves special mention. This type is used to force everything inside a polygon to be set flat at a constant height. It’s used mostly for lakes when there’s inadvertently other data inside them, such as lidar points, whose heights are not exactly the same as the shoreline and therefore prevent the water bodies from being flat. Use of the Replace SFType does incur more processing cost than normal hard or softlines so it’s best to avoid. Ideally there ought not be lidar samples in your water bodies (consider adding this as a stipulation in the contract with your data provider), but if you do, you can either use the Replace SFType to handle them or get rid of the offending points before building your terrain using the Erase geoprocessing tool.
If you’ll be producing both bare earth and first return surfaces via terrain datasets, load the lidar points into two different multipoint feature classes, a feature class for the ground points and a feature class for the above ground points. Your bare earth terrain is defined with a reference to just the ground points. Your first return terrain references the same ground point feature class as the bare earth terrain and has the additional reference to the above ground points. Yes, this means two different terrains can reference the same feature class.
Starting with ArcGIS 9.3, terrain datasets can be pyramided using one of two point thinning filters: z-tolerance and windowsize. For DEM production you can use either. If you intend to rasterize from the full resolution point set, then use the windowsize filter for terrain construction because it’s significantly faster. If you’re willing to use thinned data for analysis, which is reasonable if the lidar is oversampled for your needs, use the z-tolerance filter. While more time consuming, it’s most appropriate because it provides an estimate of vertical accuracy of the thinned representation. For DSM production use the windowsize filter with the MAX option.
Use the Terrain To Raster tool to produce your rasterized elevation model. This provides options for interpolation, output cellsize, and which pyramid level to use from the terrain.
The Terrain to Raster tool produces gridded elevation models from terrain datasets.
For interpolation, the natural neighbors options is your best bet. While not as fast as linear interpolation, it generally produces better results both in terms of aesthetics and accuracy. Set the output cellsize relative to the lidar point sample density. You won’t gain anything by using a cellsize that’s substantially smaller than the average point spacing. Also, make sure to set the analysis extent, as set through the environment, for the extraction of subsets where appropriate. The use of a snap raster can also be of use for the sake of alignment of raster outputs.
Using either Point To Raster or Terrain To Raster geoprocessing tools you can process hundreds of millions, even billions, of lidar points into hi-resolution gridded DEMs and DSMs. These can then be used with the large collection of raster tools available in ArcGIS for analysis. They’re also great for making maps (see graphic below) and, due to their simple data structure, easy to share.
Color hillshade for a DSM of Jasper County, South Carolina. Made from a terrain dataset built on top of 1.7 billion lidar points for the county.
That’s it for the second part of this series on Lidar Solutions in ArcGIS. Check back for the next part: Data Area Delineation from Lidar Points.