Analysis of terrestrial laser scanning data for Landers 1992 earthquake fault scarp

by Ramón Arrowsmith, School of Earth and Space Exploration, Arizona State University

Updated March 15, 2010.

We (UT Dallas team and me) scanned the fault scarp that formed in the 1992 Landers earthquake in the Mojave Desert California in May 2008. This is a site that Dallas Rhodes and I have studied since immediately after the earthquake (see Arrowsmith and Rhodes, 1994 and Arrowsmith and Rhodes, 200?). Here is a link to some of other pictures from the site.

Here is an AGU 2009 presentation which includes analysis of these data: AGU09Arrowsmith_TLS.ppt

Data overview

We used the UTD Riegl terrestrial laser scanner. The data were processed by the UTD team. They gave me a 9.54 million point 3 column file of the laser returns that looks like this:

E,N,Z
540969.011,3822032.969,890.643
540969.011,3822032.97,890.643
540969.092,3822032.91,890.633
540970.323,3822032.053,890.932
540970.405,3822031.995,890.942
540970.487,3822031.938,890.952
540970.569,3822031.88,890.962

Here are the original data to download. These data were bad. When the UTD guys transformed them, there were a bunch of truncation and other errors that were introduced. Here is the old version of these pages: old page.

Here are the updated data to download in compressed ascii format: EF83_2.zip.

Here are the updated data to download in compressed LAS format: Landers_TLS_data_UTM_LAS.zip.

Figure 1 shows an overview of the site and figures 2-4 zoom to the channels the changes in which we have monitored erosional change since 1992. Note that in the areas of high shot density in the areas of the channels, the 10 cm DEMs are very nice and show the channel landforms well. However, in the areas of lower shot density, the resulting DEM exhibits corduroy and contour parallel grooving which are probably artifacts of the IDW local binning algorithm, rather than the scan stitching.

Figure 1. View of a hillshade made from the optimal DEM (0.1 resolution and 0.5 m search radius; see below). Color overlay shows illustrates shots per square meter. Note the channels that we emphasized have more than 35,000 shots per square meter (4 orders of magnitude more than a typical airborne laser scan!).

Figure 2. NW channel zoom. Look here for a balloon photo of the site from 1998 of the site.

Figure 3. Center channel zoom. Look here for a balloon photo of the site from 1998 of the site.

Figure 4. SE channel zoom.

Analysis of shot density and optimal DEM local binning search radius

The problem we faced was to analyze the data with the standard digital elevation model processing tools which we use for analysis of airborne laser scanning data and also for the geomorphic investigations typical for DEMs.

For those analyses, I have found that the inverse distance weighted local binning approach (see IDW local binning algorithm) works fairly well. In this simple DEM processing algorithm, the DEM grid node is estimated by computing the inverse distance weighted (1/r²) mean of the points within a specified search radius. If there are no points, a null value is assigned. The tradeoff comes then from having a large enough search radius to find at least one point for the grid node, and not using a point so far away from the grid node as to provide an unacceptable representation of the elevation. For a given DEM grid spacing, one can compute DEMs with progressively increasing search radii, and compute the number of nodes versus search radius. This will plot as a negatively sloping curve which becomes asymptotic to some finite number of nulls. An inflection point can be qualitatively determined, and that provides the optimal minimum search radius with the maximum number of acceptable nulls (Figure 5). For this dataset and a 0.1 m DEM, the optimal search radius is about 0.5 m.

Figure 5. # nulls versus search radius in meters. Optimal search radius is about 0.5 m for this 0.1 m DEM analysis; although it is not as sensitive as it is for some airborne LiDAR cases.

Figures 6-8 show the effect of varying search radius on the resulting hillshade of a 0.1 m DEM. With a search radius of 0.1 m (Figure 2), much of the scene is null. The lower portion (to the west) was illuminated more easily from the station locations on the west side of the study area. The upper parts of the watershed were not scanned densely.

Figure 6. Hillshade map of 0.1 m DEM topography with 0.1 m search radius.

Figure 7. Hillshade map of 0.1 m DEM topography with 0.5 m search radius. This is optimal.

Figure 8. Hillshade map of 0.1 m DEM topography with 2 m search radius.

Figure 9. Hillshade map of 1 m DEM topography with 1 m search radius. Looks ok if you are far enough away.

Filling in the nulls

The figures above with 0.1 m DEM topography with 0.5 m search radius (Figures 1-4 and 7) look pretty good and are optimal as roughly defined in the last section. Nevertheless, the gaps (nulls) are annoying. What can we do? There are two ways to handle this. First is to use a new null filling module of OpenTopography Points2Grid. The second is to try to progressively build the dem out with coarser resolution data filling the gaps. I will present both approaches below.

Null filling with OpenTopography Points2Grid

The basic idea of null filling in this implementation uses a focal IDW mean - similar to ESRI's FOCALMEAN ()http://webhelp.esri.com/arcgisdesktop/9.3/index.cfm?TopicName=focalmean con(isnull([DEM]),focalmean([DEM], rectangle, 3, 3, DATA),[DEM]) but in this case we are using an inverse distance weight if the window is greater than 3x3. Fill window size options are 3, 5 and 7. I wrote more about this in this document: Exploring the new Open Topography Points 2 Grid with null filling.

I ran it on the Landers point cloud and the results are shown below.

Figure 10. Overview of hillshades of the 0.1 m DEM with 0.5 m search radius and progressive increases in the null filling by OpenTopography new Points2Grid.

Figure 11. Zoom to area with major nulls on hillshades of the 0.1 m DEM with 0.5 m search radius and progressive increases in the null filling by OpenTopography new Points2Grid.

Figure 12. Null reduction as a function of the fill window size. In this case, as we go from no filling to a 7x7 window size for filling, we reduce the number of nulls by about 170,000.

Null filling with by grid combination

The 2 m search radius (Figure 8) covers a lot of ground, but probably is not representing the topography acceptably, so now, we can take the optimal 0.1 m DEM (Figures 1 and 7) and combine it with a broader DEM (2 m with a 2 m search radius; Figure 8) to make a blended DEM?

Simple overlay

Figure 13. Hillshade map of 2 m DEM topography with 2 m search radius to be used as a background filler for the finer 0.1 optimal DEM.

Figure 14. Hillshade map of 2 m DEM topography with 2 m search radius underneath the finer 0.1 optimal DEM.

Figure 15. Hillshade map of 2 m DEM topography with 2 m search radius underneath the finer 0.1 optimal DEM zoomed to show combination.

Map algebra to compute a composite DEM

Given all of the individual analyses from above, is it possible to make a composite DEM, taking advantage of the increasing shot density towards the core of the survey?
I did this Spatial Analyst in ArcMap, with the options set as cell size (same as finest resolution DEM; 0.1 m in this case). My basic logic was to build a series of grids of the elevations (at the same resolution), starting with the interior resolved portion with the 0.5 m search radius, then take the 2 m DEM computed with the 2 m search radius and use its values where the 0.1 m DEM was null. Then add in the 5 m DEM computed with the 5 m radius beyond where the 2 m DEM was null. In Spatial Analyst--Raster Calculator these are the commands:
twombandb=IsNull([idwpt5pt1]) & [idw2_2] > 0
fivembandb=IsNull([idw2_2]) & [idw5_5] > 0
Zero out the base dem where null:
pt1pt5z = con(IsNull([idwpt5pt1]),0,[idwpt5pt1])

add together newpt1dem=[pt1pt5z] + [twombandb] * [idw2_2] + [fivembandb] * [idw5_5]

[idwpt5pt1] is the 0.1 m DEM with 0.5 m search radius.
[idw2_2] is the 2 m DEM with the 2 m search radius.
[idw5_5] is the 5 m DEM with the 5 m search radius.

Figure 12. Sequence of grids used to produced the combined final DEM. Upper plots show processed elements which were used to control the combination.

Other processing

From within ArcMap, I could barely get the points loaded (9.54 million point 3 column file). I could not compute any surfaces or a TIN.