LaDiCaoz and LiDARimager—MATLAB GUIs for LiDAR data handling and lateral displacement measurement

While many large and destructive earthquakes have occurred within the past decade alone, recurrence of those devastating events along the individual structures is fortunately rare, typically taking a hundred to several hundred years or more (e.g., Scholz, 2002). To gain a deeper understanding of the processes that control when and how earthquakes along individual structures recur requires statistical analysis of seismic records that span multiple earthquake cycles, constraining time and size of past large earthquake occurrences. Due to the prevailing absence of proper instrumental seismic data for such analysis, geologic and geomorphic evidence of seismic activity is commonly utilized to extend the seismic record into past centuries and millennia (e.g., Scholz, 2002; McCalpin, 2009), as sufficiently large earthquakes (>M6.0) disrupt and displace geologic and geomorphic units (e.g., Burbank and Anderson, 2001). These disruptions and displacements may be preserved in the geologic and geomorphic record, enabling size and age estimation for the causative earthquakes via geochronological methods and earthquake scaling relationships (e.g., Wells and Coppersmith, 1994; Sowers et al., 2000). This in turn enables rupture sequence reconstructions for surface rupturing events and the formulation of earthquake recurrence models that describe large earthquake recurrence along a given fault.

Surface-slip reconstructions for the great M7.8 Fort Tejon earthquake of 1857 along the San Andreas fault (SAF)—a structurally relatively simple and easily accessible major strike-slip fault that separates the North American and Pacific plates—and preceding seismic events along the 1857 surface rupture trace, have been of particular importance in the formulation of current earthquake recurrence models, shaping the understanding of seismic fault behavior along with earthquake forecasting and seismic hazard assessment (e.g., Sieh, 1978; Schwartz and Coppersmith, 1984; Sieh and Jahns, 1984; Sieh, 1996; WGCEP, 2008; Field et al., 2009; Zielke et al., 2010; Akciz et al., 2010; Zielke et al., 2012). Those surface-slip reconstructions along the 1857 rupture trace are largely based on offset measurements of ephemeral stream channels, laterally displaced as they cross the SAF trace (Fig. 1). The underlying assumption is that channel incision events in this region generally occur more frequently (decadal time scale) than large earthquakes recur (centennial time scale) (e.g., Wallace, 1968; Sieh, 1978; Sieh and Jahns, 1984; Grant and Sieh, 1993; Lienkaemper, 2001; Grant Ludwig et al., 2010; Zielke et al., 2010; Zielke et al., 2012). Consequently, the smallest observable offsets correspond to the most recent earthquake (here the 1857 Fort Tejon event), and successive larger offset groups record the cumulative slip of prior events (Fig. 1).

Investigations concerning fault zone characterization and coseismic surface-slip distribution commonly rely on aerial photography and field investigations as the main data sources (e.g., Wallace, 1968; Sieh, 1978; Sieh and Jahns, 1984; Lienkaemper, 2001; Rockwell et al., 2002; Haeussler et al., 2004). High-resolution digital elevation models (DEMs), generated from airborne laser swath mapping (ALSM) or terrestrial laser scanning (TLS), also known as light detection and ranging (LiDAR), have recently become an additional powerful data set, available to such studies (e.g., McCalpin, 2009). During data acquisition, the scanner (aboard an aircraft or mounted on a tripod) emits laser pulses that travel to the surface, where they are reflected. Part of the reflected radiation returns to the scanner, is detected, and stops the time counter that was started when the laser pulse was sent out. The intensity of the returning pulse is then recorded along with its time-of-flight.

The latter is converted to a distance between scanner (source) and surface (reflector), by utilizing the speed of light. For airborne systems, a combination of onboard global positioning system (GPS) and inertial navigation system (INS) allows further conversion of this distance to provide absolute geographic coordinates (latitude, longitude, and elevation) for the reflector (the surface). Note that each outgoing pulse may generate multiple returning pulses that differ in travel time and/or intensity, for example due to partial pulse reflection by vegetation cover. Classification of those returning pulses, for example by travel time and/or intensity, permits virtual deforestation of the scanned area, resulting in “bare earth” DEM (Fig. 2). Generally, shot densities (pulse returns per m²) and thus DEM resolution vary as a function of pulse frequency (200–400 kHz), scanner velocity (when airborne), and source-reflector distance. For example, average shot return densities of 3–4 m⁻² in the “B4” airborne LiDAR data set (Bevis et al., 2005; available at www.opentopography.org) that covers the southern San Andreas fault (SAF) allow generation of DEM with <0.5 m grid size, permitting identification of submeter-scale tectonogeomorphic features (Fig. 1). Thus, the possibility of virtual deforestation and the generally very high spatial resolution of LiDAR data makes them a valuable asset for geologic, geomorphic, and paleoseismic investigations, exemplified by studies along the SAF, the eastern California shear zone, and other fault zones around the world (e.g., Hudnut et al., 2002; Sherrod et al., 2004; Baum et al., 2005; Schulz, 2007; Oskin et al., 2007; Hilley et al., 2010; Arrowsmith and Zielke, 2009; Gold et al., 2009; Baran et al., 2010; Ganev, 2010, personal commun.; Zielke et al., 2010, 2012).

Along with this article, we provide two MATLAB-based graphical user interfaces (GUIs)—LiDARimager and LaDiCaoz—for LiDAR data processing and visualization. Both GUIs perform well on current standard desktop computers, and only a MATLAB license (no additional MATLAB toolboxes needed) is required. Naturally, the maximum DEM size that LiDARimager and LaDiCaoz can process depends on availability of memory. Digital elevation models with <10⁷ grid points work well in both GUIs for computers with ≥2Gb of memory. While most of the available visualization options of LiDARimager and LaDiCaoz may also be found in standard GIS packages (e.g., ESRI's ArcGIS), other GUI features present unique contributions: a special feature of LiDARimager allows importing DEM at lower-than-original resolution, enabling processing of large (≥10⁸ grid points) DEMs (we further discuss the motivation and possible application of this option in the next section). A key feature of LaDiCaoz, also not available in standard GIS packages, is DEM slicing and (laterally) back slipping for visual offset reconstruction assessment, improving measurement accuracy and precision. The presented GUIs therefore enable comprehensive LiDAR data processing, manipulation, and visualization as well as lateral offset measurement in a single software package. Recent studies have used those tools to reevaluate the surface-slip distribution along the 1857 Fort Tejon surface rupture trace (Zielke et al., 2010; Zielke et al., 2012) and to measure surface displacements along the Garlock fault (Ganev, 2010, personal commun.). While both GUIs were developed for LiDAR data processing, other gridded DEM data sets (e.g., U.S. Geological Survey [USGS] DEM, Space Shuttle Radar Topography Mission [SRTM] DEM data) may be processed as well when provided in the correct input format (see the tutorial [the AOM.doc file] in the Supplemental File¹), further adding to the functionality and range of potential applications of the provided tools. Following, we describe the functionality of LiDARimager and LaDiCaoz and present a workflow for channel offset measurements that utilizes these tools. We then discuss the repeatability and uniqueness of channel offset measurements, made with the LaDiCaoz and LiDARimager. The GUIs (LiDARimager.p and LaDiCaoz.p), the manuals (AOM.doc), and sample LiDAR data (Sample1.asc and Sample 2.asc) may be found in the Supplemental File (see ^{footnote 1}). Additionally, tutorial videos may be downloaded from http://stockdale.sese.asu.edu/Geosphere_AOM.

LiDARimager (Fig. 3) was created to facilitate the identification of fault-zone structure and offset geomorphic markers for further analysis in LaDiCaoz. For this purpose, a range of LiDAR-based DEM derivates may be generated, namely elevation, hillshade, surface-aspect, and surface-slope plots (Fig. 4). Those plots may be combined (e.g., draping a semitransparent slope plot over an opaque hillshade plot), overlain with elevation contours and UTM grid, and presented either in map view or freely adjustable oblique view (Fig. 5). Visualizations may be stored either as*.jpg or *.kmz files, the latter for import into GoogleEarth (Fig. 6). LiDARimager may further be used as a cropping tool: user-defined sections of DEM (imported at lower-than-original resolution) may be extracted at original resolution for further processing and visualization—a helpful option when dealing with large (≥10⁸ grid points) DEMs (also see Fig. S2 in the AOM.doc file in the Supplemental File [see ^{footnote 1}]). For example, you may import a large LiDAR scene at lower-than-original resolution and search for potentially offset features (lower resolution is chosen because of computer memory constraints; it is generally sufficient to identify sites that potentially exhibit an offset feature). For those potential sites, smaller sections of the LiDAR scene may then be extracted at original resolution for further inspection and potentially offset measurement of the feature. The quick and easy generation of DEM visualizations with LiDARimager has proven effective during identification of fault trace and potentially offset features (Fig. 7) in recent LiDAR-based surface-slip reconstructions along the south-central SAF (Zielke et al., 2010, 2012).

The following LiDARimager workflow reflects the adopted procedures for fault trace and offset identification. Note, however, that application of LiDARimager is neither limited to tectonogeomorphic studies nor southern California, but may contribute to a wide range of studies and study areas that utilize LiDAR data visualizations. Initially, we generated a number of hillshade plots of the fault zone (using different illumination angles) to identify primary and secondary fault traces (e.g., Figs. 7A and 7B). The distinction between primary and secondary traces was based on the number of subparallel fault strands and their relative geomorphic expression. Depending on the absolute level of surface expression, we assigned ratings of certain, uncertain, inferred, and queried to individual sections of the fault trace(s) (Fig. 7B). We then searched hillshade plots, contour plots, as well as slope-shade plots (slope draped over hillshade) along the fault traces for offset geomorphic markers (e.g., offsets and bends in sublinear features such as stream channels and alluvial fan edges as they cross the fault zone). For clarity of language, the remainder of this article generally refers to stream channels when addressing offset geomorphic markers (although sublinear features, other than stream channels, may be investigated in essentially the same manner). For each identified offset channel, we assigned a quality rating of either high, high-moderate, moderate, moderate-low, or low, depending on structural and geomorphic complexity of the surrounding fault zone and topography (e.g., Table 1; Sieh, 1978; Lienkaemper, 2001). This rating scheme was used to emphasize reliable measurements and vice versa. Unreliable features were not included in the record. To guide the identification of these unreliable sites, we followed the criteria suggested by Sieh (1978): “The principal categories for unreliable sites are (1) stream channels deflected around uphill-facing scarps; (2) irregular channels displaced across several-meter-wide fault zones; and (3) possibly offset features at localities where the fault trace cannot be precisely located.” Elevation, surface-aspect, surface-slope, and contour plots additionally facilitated the differentiation between reliable and unreliable offset markers. For sufficiently reliable offset markers, we cropped the DEM in LiDARimager (see the AOM.doc file in the Supplemental File [see ^{footnote 1}]) for further processing in LaDiCaoz.

LaDiCaoz—which stands for lateral displacement calculator—was developed to allow quick and easy-to-reproduce measurements of tectonically offset, sublinear geomorphic features (e.g., fluvial channels) along the 1857 Fort Tejon earthquake surface rupture trace (Fig. 8; Zielke et al., 2010; Zielke et al., 2012). It has been the main tool in these investigations, and the following LaDiCaoz workflow with corresponding figures reflects the adopted procedures. The presented workflow example is largely based on the northwestern channel at Bidart Fan, Carrizo Plain (Fig. 1A). We want to remind the reader that application of LaDiCaoz is not limited to this region or fault but may be used to measure practically any type of lateral displacement or deflection in a gridded data set, including offsets along other strike-slip fault systems. Offset measurements with LaDiCaoz consist of four steps that will be addressed in the following subsections, namely (1) the manual mapping of fault and offset channel, (2) an automated offset calculation, (3) back slipping to reconstruct pre-earthquake topography for visual offset measurement assessment, and (4) an automated production of output data.

For each identified offset feature, we created multiple base maps (contour plots and hillshade plots with varying illumination parameters and slope-shade plots) in LaDiCaoz to gain further understanding of the site's morphology, including the precise position of fault and channel trace (Figs. 9A–9C). This approach guided us in further sorting out unreliable sites (see respective criteria in previous section). If an offset feature was considered sufficiently reliable, we mapped the fault trace and determined the location (distance normal to fault trace) of upstream and downstream topographic profile (Fig. 9D). The profile positions were chosen to be close to the fault trace but outside of the geomorphically overprinted fault zone (e.g., Wallace, 1968; Sieh, 1978; Lienkaemper, 2001; Ouchi, 2004). Figure 9E shows the topographic profiles along the red and blue lines in Figure 9D. We typically crop one profile (the blue one) to approximately the cross-sectional extent of the channel that will be reconstructed. Cropping is justified because the goal is to determine the offset of the corresponding geomorphic feature (the stream channel cross section) and not the topography surrounding it. Cropping also gives meaning to the goodness of fit estimate (GoF—a cross-correlation metric that will be introduced in a following section) that defines how well two cross-sectional profiles match as a function of back slip.

The utilized offset calculation procedure requires that both profiles are projected onto the fault plane. Depending on channel segment obliquity (relative to the fault trace), channel segment thalweg gradient, and fault trace–profile distance, this projection results in an apparent horizontal and vertical displacement (Fig. 10). To account for apparent horizontal displacements, we traced the orientation of upstream and downstream channel segment at the respective profile locations (Figs. 9D and 10C). The apparent vertical displacements were discarded because LaDiCaoz is only measuring lateral displacements.

We want to emphasize that a careful identification of (1) fault position and orientation, as well as (2) profile position and respective channel segment orientation is crucial for meaningful offset calculations and the main source of measurement uncertainty. Similarly important is a sound understanding (i.e., assumption) of the pre-earthquake channel morphology that is to be reconstructed.

LaDiCaoz iteratively determines the optimal offset, by incrementally changing (1) vertical shift of the blue profile, (2) vertical stretch (z-factor) of the blue profile, and (3) horizontal position of the blue profile: (a) Naturally, the thalweg elevation at upstream and downstream profile location will be different (otherwise channel gradient would be zero, and no flow would occur; Fig. 10B). To account for that and thus allow for a better correlation of both cross-sectional profiles, LaDiCaoz iteratively changes the elevation of the blue profile by shifting it upwards or downwards (shift direction depends on whether the blue profile crosses the upstream or downstream section of the channel; Fig. 11A). (b) Upstream and downstream channel sections may evolve morphologically in different ways, for example due to tectonic activity along the fault that alters channel gradient and cross-sectional profile around the fault zone. The downstream segment may become abandoned and start to degrade diffusively, lowering the profile's cross-sectional relief. Alternatively, the downstream segment may be captured by another stream, causing morphological adjustment of its cross-sectional profile to the new hydrological conditions. To account for these scenarios, LaDiCaoz iteratively changes the z-factor to vertically stretch or unstretch the blue profile as a proxy for different morphologic evolution (Fig. 11B): First, the minimum elevation within the blue profile is subtracted from the blue profile. The resulting profile is then multiplied with the z-factor. Lastly, the minimum elevation of the original blue profile is added to the stretched profile. Note that it is possible to enter negative stretch values as well, enabling inversion of the (blue) profile. This approach might be taken when aiming to match a gully with a corresponding alluvial fan apex. (c) The third iterative dimension is the relative horizontal displacement—the value of interest. LaDiCaoz iteratively changes the horizontal position of the blue profile relative to the red profile (Fig. 11C). Note that the increment size for horizontal shift and the along-profile topographic resolution are set to be equal to simplify the offset calculation (see the AOM.doc file in the Supplemental File [see ^{footnote 1}]). To increase computational efficiency, it is possible to limit the range of these morphologic parameters over which LaDiCaoz is iterating.

We introduce a simple cross-correlation metric to determine the optimal offset estimate. LaDiCaoz iterates over the aforementioned channel morphology and position parameters (Fig. 12B, middle) and determines the summed absolute elevation difference Σ[Δ(elevation)] between both profiles—exemplified by gray area in Figures 11A–11C—for each parameter combination. We further define a goodness of fit parameter (GoF) that equals the inverse of the summed absolute elevation difference 1/Σ[Δ(elevation)] (Figs. 11 and 12B, bottom). The optimal horizontal offset is defined by the parameter combination that results in the least mismatch between both profiles, in other words the parameter combination for which the summed elevation difference between the profiles (area between profiles) has its minimum (Fig. 11), and GoF has its maximum. The blue profile (Fig. 12B, middle) is then horizontally back slipped by the optimal offset amount (also applying the corresponding vertical shift and stretch to the blue profile) and overlain with the red profile to visually assess the cross-sectional quality of the correlation (Figs. 11C and 12B, top).

We then back slipped contour plots and hillshade plots of the topography by the optimal horizontal slip estimate to visually assess the quality of the slip reconstruction—closely inspecting the reconstructed channel trace (its thalweg and both riser edges) for potential bends (Figs. 12C–12F).Their existence in the back-slipped topography indicates that an either too high or too low offset amount was used for reconstruction, while the real offset (back slip) amount should have removed those bends. We generally used contour plots for this offset quality assessment because the potentially existing bends in thalweg and riser edges are easier to identify here. Also very good are slope-shade plots with a contour-plot overlay. Figures 12D and 12F, respectively, present examples for marginally too low and too high back-slip values. In both examples only one of the riser edges is aligned across the fault zone, while the other one is not. A common cause for mismatches in the back-slipped topography—while the cross-sectional profiles show a very good match—is an incorrect representation of fault and/or channel segment trend and/or position and therefore an incorrect channel profile projection into the fault plane (Fig. 10). In such a scenario, we reexamined actual and mapped fault location and trend as well as actual and mapped channel trend and adjusted the respective values parameters accordingly.

It is therefore of crucial importance to have a good understanding of the fault trace location and the channel morphology prior to the earthquake(s) so that appropriate channels sections are matched. Topographic back slipping presents an important quality control to test whether the automatically determined offset amount is actually comparable to the real offset that has occurred along the fault. Back slipping is thus a key feature of LaDiCaoz.

Aside from defining the optimal offset estimate, we also determine minimum and maximum offset value, capable of reasonably well reconstructing the pre-earthquake topography. This is done by trial and error. We back slipped the topography by values bracketing the optimal slip amount and closely inspected how well these offset values reconstructed the initial topography (Figs. 12C–12F).

While measuring offsets of laterally displaced channels is methodologically simple, difficulties in estimating the exact fault zone position and orientation as well as estimating the pre-earthquake channel morphology may severely complicate this procedure. These difficulties are enhanced when the displaced channels exhibit high-frequency sinuosity, have only short sublinear channel segments, experience distinct differences in geomorphic evolution of upstream and downstream segment, cross the fault zone at a low angle, or have severe post-earthquake geomorphic overprinting.

Because difficulties in identifying pre-earthquake morphology and fault-zone structure find their expression in the offset estimate (via profile projection onto the fault trace), the question arises of how repeatable those surface-slip measurements are. How important is the human factor in such an approach to surface offset measurement? Different interpretations of pre-earthquake channel morphology and orientation will be reflected in different offset estimates. Depending on channel complexity and operator, the resulting offset estimates may have wide range and deviate distinctly from the actual tectonic displacement. We briefly addressed this problem in a pilot study, involving seven colleagues who were asked to measure the offsets of 11 selected channels (Fig. 13), using LaDiCaoz and the tutorial and videos provided with it. Note that LaDiCaoz and LiDARimager did not feature the draping option at the time of the pilot study (this option was included following a suggestion in the peer-review process). Everyone provided their best offset estimate as well as maximum and minimum offset values that were considered to reasonably well reconstruct the pre-rupture morphology. Also provided was the quality rating that was assigned as well as channel and fault trace orientation. Figure 14 and Table 2 present those corresponding results. When comparing the results for the 11 channels, we find good to very good correlation for some features (e.g., Sieh31, ZA6650a, and ZA8817a), but rather large discrepancies for others (e.g., Sieh34 and Sieh77b). Comparison of channel complexity and measured offset variability (Fig. 13 and Table 2) shows that the assigned quality rating (and respective variability) is reflected in standard deviation of offset estimates (channels with a distinctly good quality rating are highlighted in gray in Table 2). In other words, for channels with a lower quality rating, offset estimates varied notably more than for channels with higher quality rating, reflecting different assumptions of the pre-earthquake channel morphology. Additionally, offset variability also exhibits a dependency on total offset amount as well as feature scale—a small, discrete channel that is offset by a small amount will generally result in smaller variability (e.g., ZA6650a), compared to a larger, less discrete channel that was offset by a larger amount (e.g., ZA6917b). Then, interpretations of what the pre-rupture morphology may deviate significantly from person to person. Figure 15 provides map and oblique views of the back-slipped topography for channel “Sieh34” (the northwestern of the two displaced channels). Minimum and maximum reported offset values were 3.2 m and 9.1 m, respectively (Fig. 15A and 15C). Back slipping by 3.2 m aligns the southeastern of the two channels quite well, but leaves a distinct bend in the Sieh34 channel. This leaves us to suggest that the individual performing this measurement was accidentally measuring the wrong channel offset. Applying a back slip of 9.1m also leaves a distinct bend in channel Sieh34, at least under the assumption that the pre-earthquake flow direction near the fault trace was essentially perpendicular to the fault trace. Inspection of the oblique views of this back slip leaves us to suggest that the individual who made the 9.1 m measurement assumed a different pre-earthquake flow direction—distinctly oblique to the fault trace—and matched channel profiles relatively farther away from the fault trace. With this assumption of pre-earthquake morphology, channel alignment may appear reasonable. However, we suggest that the 5.8 m back slip better aligns both channel risers and channel thalweg across the fault zone (Fig. 15B), as we (along with most of the other colleagues who performed the back slipping) assume that flow direction near the fault trace is essentially perpendicular to it. The presented example, however, demonstrates the difficulties that may arise in those offset measurements. Nonetheless, for most cases back slipping clearly reveals which offset estimates are too high or too low and therefore increases both measurement accuracy and precision.

With this publication, we provided two computational tools (the MATLAB-based GUIs LiDARimager and LaDiCaoz) to process and visualize LiDAR-derived DEM data and to measure the lateral displacements of offset geomorphic markers. These tools were initially developed for identification and measurement of tectonically displaced geomorphic features along the rupture trace of the 1857 M7.8 Fort Tejon earthquake (southern California). They are not, however, limited to this study region and corresponding scientific questions, but may find application in a range of studies that utilize LiDAR data visualizations. They may further be used to allow repeatable and self-consistent measurements of laterally displaced geomorphic markers.

While the presented methodology utilizes a semiautomatic offset calculation algorithm, we want to stress the importance of the operator in the measurement process. Carefully placing fault and projection lines as well as carefully inspecting the back-slipped topography are crucial aspects that rely on the tectonogeomorphic experiences of the operator. Back slipping the topography by the computed amount helps to ensure that the computed offset values are as close as possible to the actual offset amount. The latter option enables operators with only a basic understanding of tectonic geomorphology to make meaningful measurements.

We would like to thank the participants of the small sensitivity study on offset repeatability for their efforts. The “B4 Project” (www.earthsciences.osu.edu/b4) collected LiDAR point-cloud data along the southern San Andreas and San Jacinto faults. Data acquisition and processing were performed by the National Center for Airborne Laser Mapping (NCALM). The project was led by Ohio State University and USGS with funding from the Division of Earth Sciences Geophysics program at the National Science Foundation (NSF). Optech International contributed the ALTM3100 laser scanner. University Navstar Consortium (UNAVCO) and Southern California Integrated GPS Network assisted in GPS ground control. This work was supported by the NSF Tectonics program and the Southern California Earthquake Center (SCEC). SCEC is funded by NSF Cooperative Agreement EAR-0529922 and U.S. Geological Survey Cooperative Agreement 07HQAG0008. We would also like to thank the editors and reviewers for their helpful contributions. We will be thankful for further suggestions on potential improvements for tools and documentation. The source code of both GUIs may be provided upon request.