Illuminating geology in areas of limited exposure using texture shading of lidar digital terrain models Open Access

The mapped distribution of rock units and orientation data is essential for establishing the geologic structure of any area. Without it, we cannot define a fold’s geometry, the slip on a fault, or the map thickness of stratigraphic units. Many geologists focus their studies in arid regions because good exposures facilitate the collection of field data. Furthermore, satellite images in such areas can commonly be used to establish continuity from one outcrop to the next. Of course, interesting geology is not limited to arid regions. In areas of soil cover, vegetation, and even anthropogenic surface disruption, what modern technologies can we use to extract the maximum amount of structure data? Judicious use of remote-sensing data not only facilitates the collection of large amounts of data but can also guide field geologists to areas and structures where their limited time on the ground will be best spent.

Hillshade images have become a standard method for displaying digital terrain model (DTM) data, yet the interpretation of simple hillshading is plagued by the bias produced when a single artificial illumination direction is used. In this paper, we demonstrate how to use DTMs collected with lidar and post-processed with a novel, directionally independent filtering method known as texture shading, not yet widely used in geology, to image subtle surface variations related to compositional stratification and metamorphic foliations beneath heavily vegetated regions. Texture shading has found use in cartography (e.g., Patterson, no date) and is one of several advanced DTM visualization schemes developed to overcome the limitations of hillshading (see list in Kennelly et al., 2021). To date, these techniques have found geological application primarily in geomorphology and activetectonics studies. Here, we emphasize their application to bedrock geologic mapping.

Several case studies are presented for areas in the northeastern United States (Fig. 1). The first is located at the front of the Valley and Ridge fold- and-thrust belt in north-central Pennsylvania (Gwinn, 1964), immediately south of the city of Williamsport. This area is well mapped geologically, has complete lidar coverage with 1 m resolution, and presents a range of rock types. The second case study focuses on a region of metamorphic rocks in the Berkshire massif of western Massachusetts. The third is from the eastern margin of the Green Mountain massif in southern Vermont. Texture shading can also provide additional insight into surficial, depositional features as illustrated in the fourth case study focused on glacial features near Ithaca, New York. Each of these examples illustrates the value of texture shading for analysis of structures or surficial features in vegetated regions.

Topographic data collected with airborne lidar sensors have revolutionized many parts of geology, especially those concerned with the evolution of the Earth’s surface such as geomorphology and active tectonics (e.g., DeLong et al., 2010; Meigs, 2013). To date, applications of lidar to bedrock geologic mapping, the focus of this paper, have lagged behind. Light pulses from an airborne laser are reflected by near-surface features, such as vegetation, and the last returning pulse is taken to represent the bare earth. The ability to extract a bare-earth model, or DTM, is what makes lidar so valuable to the geologist and what distinguishes lidar from other high-resolution remotesensing methods such as structure from motion. Lidar data have pixel sizes reflective of the scale of the layering and subtle surface features that are of interest to geologists.

The most common method of presentation of DTMs is hillshading. This method, however, has several limitations. The terrain features perceived by the viewer depend on the azimuth and vertical angle of the artificial illumination, even though it is the illumination and artificial shadows that give the image their three-dimensional appearance. Slope and relief dominate a hillshade image at the expense of more subtle features that may be the geologist’s primary interest. In addition, features oriented at a low angle to the artificial illumination azimuth are difficult to discern.

To overcome these limitations, a variety of imaging modifications and post-processing enhancements have been applied. Multiple sources of illumination can be used (e.g., Kennelly and Stewart, 2014), and hillshading can be combined with planimetric or profile curvature to further enhance subtle details (Kennelly, 2008). Additional image-processing techniques include slope-shade images and red relief image maps, which combine topographic slope (in red) with an openness metric where the sign indicates convexity (negative) or concavity (positive) (Chiba et al., 2008; Zielke et al., 2015).

In this paper, we use a cartographic post-processing technique, described in more detail below, known as “texture shading.” Rather than process lidar point-cloud data ourselves, we use the 1-m-resolution lidar DTMs that are available from the U.S. Geological Survey (USGS) 3D Elevation Program (3DEP) and the state government GIS data websites. The USGS data are of level 2 quality, which corresponds to ≥2 pulses/m2 and pulse spacing of ≥0.71 pulse/m. The root-mean-square difference in the vertical direction is ≤0.06 m.

Lidar DTMs are large—the one used for the north-central Pennsylvania study is 2.35 gigabytes in size covering an area of 470 km2—and the images of the topography produced from those data are likewise large, close to 500 million pixels. Thus, we convert these large images into MBTiles files (Fischer et al., 2018), which tile the maps so that the viewing program never needs to load the entire, full-resolution image into memory at any one time. In this study, we tile the base-map image with zoom levels between 12 (~76 m/pixel) and 17 (~1 m/pixel) using the program MapTiler (https://www.maptiler.com); the resulting MBTiles SQLite database contains almost 12,000 individual, 256 × 256 pixel tiles.

As described in more detail below, we use the program GMDE (Allmendinger, 2020) to visualize and analyze these data. GMDE can read both raster base-map images (e.g., of geologic maps of an area) and MBTiles files and can extract elevations at any point from a GridFloat- or ESRI BIL-formatted DTM (Allmendinger, 2020). GMDE does not hold the entire DTM in memory; instead, it reads the necessary elevation from the DTM file on disk using a direct binary read, a process that is essentially as fast and much less RAM intensive than placing the entire DTM in memory at once.

Leland Brown (2010, 2014, 2019) introduced a novel approach to imaging a DTM known as texture shading. Unlike simple hillshade images, texture shading has no directional significance because it does not rely on artificial illumination. The technique uses a linear, fractional Laplacian operator, which is a more complicated relative of the standard Laplacian operator used in edge-detection algorithms. In the spatial domain, the fractional Laplacian operator is (Brown, 2014):

where the power d ranges from 0 to 2 and k_d is a constant. In the above equation, x and y represent the horizontal position and u and v the components of a displacement vector from (x,y). The power law in the denominator produces a type of scale invariance that Brown (2014, p. 5) described as: “smaller and smaller details should have lower and lower contrast, where the ratio of contrast reduction should have a fixed relationship to the ratio of feature sizes. For example, if two similar features with size ratio 2:1 have a contrast ratio of 1.5:1, the same contrast ratio should hold regardless of whether the features are two large mountains or two small bumps on a mountain.” This scale invariance enhances the continuity of the bands of light and dark color in that a wide band, related to a thicker, more highly erodible rock unit, will persist as a wide band across the image as long as the underlying rock layer responsible for the banding maintains its thickness and lithology. As one zooms in, narrower bands, corresponding to thinner layers, become more apparent, just as in a fractal image. This means that a single texture-shade image can be used for any desired scale of magnification. The value of d weights the importance of nearby and distant features; larger values of d place greater emphasis on details. To overcome the difficulties posed by applying this continuous function to discrete data (i.e., a DTM), in Brown’s (2014) algorithm, a Fourier transform is applied to render the data continuous in the frequency domain, the operator applied, and the results transformed back into the spatial domain.

The practical effect of this type of processing is that convex-up edges or curves in the DTM are light in color whereas concave-up surfaces are dark in color. Thus, local ridges are light and valleys are dark (Fig. 2). Whereas hill-shading emphasizes slope and relief, texture shading enhances changes in slope, which are commonly more pertinent to the underlying structure and composition.

In this study, we use texture-shade images produced with the commercial program Natural Scene Designer (https://naturalgfx.com); the inventor of texture shading also provides Python code for free on his web page, https://app.box.com/v/textureshading/. The images presented here are rendered with a d = 1.3–1.4, values that were empirically derived as best suited to capturing bedding in the lidar DTMs. Two detailed examples from our study area in north-central Pennsylvania compare photographic, hillshade, and texture-shade images for the same areas. Figure 3 shows an area of gently dipping bedding. Some of the layering is evident in the hillshade image, but in most other areas it is washed out because of the illumination direction; using multiple artificial light sources might provide additional detail. The texture shade provides an exquisite image of changes in slope that correspond to changes in the erodibility of the underlying rocks. An area of more complicated structure is depicted in Figure 4. Again, some structure is apparent in the hillshade image, but much more complicated fold and fault patterns are revealed by the texture-shade image.

With an image of a planar surface crossing topography, the orientation of that feature can be determined from three or more non-collinear points on the surface. We use GMDE to calculate orientations (Allmendinger and Judge, 2013; Allmendinger, 2020) from texture-shade images. Because the lidar DTM has a grid spacing of 1 m, the three (or more) points need not be widely separated, although the points should have separation considerably greater than the width of the light or dark band being used. If the elevation differences between points are small, then the spacing should not be large. The larger the horizontal spacing, the more nearly collinear the points are likely to be, especially for beds of significant dip.

Though GMDE can find the best fit to a plane using more than three points, an over-constrained problem that allows one to estimate uncertainties, our workflow takes advantage of another capability of GMDE, that is, the ability to calculate the intersection of a planar feature of known orientation with the Earth’s surface, which we call contact projection (Allmendinger, 2020). First, the user selects three points along a light or dark band, consistently locating the points to sample the same part of the band. GMDE uses these three points to calculate an orientation. The user can then project the surface trace of the plane based on that orientation. If the calculated trace matches the observed trace, the orientation is accepted (Fig. 5). If, however, the calculated trace does not match the observed trace, the user can adjust the position of any of the three points with a click in the map window and GMDE will recalculate the surface trace automatically. This interactive approach will rapidly converge on an acceptable solution, if there is one, or will make it clear that the observed trace does not correspond to a planar feature (Fig. 5).

Three-point calculations are complementary to field measurement of orientations in several ways. Most obviously, even with the resolution of lidar data, the area of a contact surface sampled in a three-point calculation is considerably larger than that captured by a strike and dip of a bedding plane measured in the field, so the measurements are not identical. Three-point calculations can determine the dip of a low-angle plane (dip <5°) much more accurately than a field measurement, although small errors in either method for sub-horizontal strata result in large strike errors. Conversely, small errors in locations of the three points in steeply dipping strata are likely to result in larger dip errors.

The Pennsylvania Valley and Ridge province and the adjoining Appalachian Plateau expose one of the world’s classic thin-skinned fold-and-thrust belts (Gwinn, 1964; Faill, 1998; Pohn, 2000; Sak et al., 2012; Mount, 2014). Although most major thrust faults in this part of the Appalachian orogen are blind, a distinct morphological and structural boundary, the Allegheny structural front, separates very low-amplitude folds that overlie a décollement in Silurian Salina Group evaporites beneath the plateau from the high-amplitude folds overlying the deeper décollement of the Valley and Ridge.

The frontal structure of the Valley and Ridge is the Nittany anticline. The northeastern termination of the Nittany anticline south of Williamsport, Penn-sylvania, makes an excellent case study (Fig. 6). The region is well mapped geologically, with copious field measurements of strikes and dips (Faill et al., 1977; Faill, 1979; Inners, 1993). It presents a range of orientations and different lithologic sequences, while exhibiting vegetative and agricultural coverage typical of the north-central and northeastern Appalachian Mountains.

The major east-plunging anticlines are beautifully outlined by the massive, resistant quartzite of the Silurian Tuscarora Formation (Fig. 7) (Faill et al., 1977; Faill, 1979). The cores of the folds expose rocks as old as Middle Ordovician and as young as the Upper Devonian. Overall, this part of the Paleozoic sequence is dominated by clastic strata with a few notable carbonate units including the Tonoloway, Keyser, and Onondaga Formations. The clastic rocks vary in grain size from shales and siltstones to coarse-grained sandstone.

The Nittany anticline has a steep northern limb with beds locally overturned and a gentle southern limb with dips rarely exceeding 30°S. The White Deer syncline immediately to the south (Fig. 6) is a broad open structure that in its extreme southern limb attains dips as high as 60°–70°. Much of the syncline is obscured by the broad flood plain of the Susquehanna River, though outcrops exist in the Muncy Hills east of the river.

In a region of ~450 km2, we calculated the bedding orientation using three-point calculations and surface contact projection at numerous localities on the texture-shade base map (Fig. 6). In all, 235 orientations were determined by three-point calculation over a time span of two days. In the same area, 364 orientations measured by field geologists, probably measured over a time span of many months, were reported on published maps (Faill et al., 1977; Faill, 1979; Inners, 1993).

Clearly, not all stratigraphic units are equally well represented in the orientation data collected from the texture-shade image (Fig. 8). There are two reasons for this: (1) the spatial distribution of some units coincides with the flood plain of the Susquehanna River and are largely covered, or, more interesting for our study; and (2) the bedding characteristics of some formations simply are not well imaged by texture shading. For example, massive shales and siltstones of the Harrell and Mahantango Formations (units Dh, Dmt, and Dmh; Fig. 6) yielded no texture-shade orientations at all despite the latter constituting 9.5% of the total field geologist–measured orientations in the same area. In contrast, formations characterized by interbedding of different lithologies such as the Trimmers Rock, Willis Creek, Rose Hill, and Bald Eagle Formations constitute the majority (55%) of the total orientations measured on the texture-shaded lidar DTM.

It is a bit surprising that 30% of the texture-shade orientations were measured in the Tuscarora and Juniata Formations, both dominated by thick-bedded quartzite highly resistant to weathering. Apparently, there are just enough thin interbeds of shale and siltstone to provide the necessary changes in slope to be captured by texture shading. It should also be noted that many of the measurements on the texture-shade image in these units were judged qualitatively to be of fair or poor quality, compared to a unit like the Trimmers Rock Formation where a significant majority of the measurements were judged good or better.

A comparison of orientations determined with three-point calculations on the texture-shade image with geologist-measured orientations on the outcrop shows that the former captures very well the structure of the anticline (Fig. 9). In fact, the cylindrical best fit, a type of principal-components analysis, shows that orientations measured on the texture-shade images show less dispersion and a tighter great circle than the outcrop measurements.

Geologic mapping in metamorphic rocks is challenging because the stratigraphy is commonly not well characterized, and poor exposure inhibits the ability to establish continuity between sparse, and commonly small, outcrops in regions where orientation and lithology can vary over short distances. A lack of fossils, unrecognized faulting and folding, and poor preservation of primary sedimentary structures, including bedding, are responsible for the many stratigraphic-structural controversies for which the New England Appalachians are justly famous. The orientation of bedding is not always discernible in metamorphosed sedimentary rocks, and in some cases bedding measurements must be regarded as interpretations rather than observations. In cases where bedding can be reliably identified, orientation measurements commonly vary widely due to small-scale folding and thus may not constrain map-scale lithologic boundaries. Glacial deposits commonly obscure bedrock exposure in northern North America, adding to the difficulties in geologic mapping.

The Berkshire massif is a north-south–oriented structure in western Massachusetts in a region affected by both the Ordovician Taconic and Devonian Acadian orogenies (Zen et al., 1983). The western frontal thrust of the Berk-shire massif transported Mesoproterozoic basement rocks and unconformably overlying Ediacaran to Cambrian clastic rocks over Cambrian to Ordovician carbonate rocks in Massachusetts (Zen et al., 1983; Fig. 10). Exposures of deformed and metamorphosed gneisses of the Mesoproterozoic basement rocks and the clastic cover are locally excellent but widely scattered in the hanging wall of the Day Mountain thrust (Fig. 10). Dolomitic and calcitic marble of the Cambrian to Ordovician Stockbridge Formation is poorly exposed in the footwall of the Day Mountain thrust (Ratcliffe, 1984).

Ratcliffe (1984) and Karabinos and Pierce (2017) interpreted the stratigraphy and structure in the hanging wall of the Day Mountain thrust sheet differently. Karabinos and Pierce (2017) mapped the Dalton Formation with fewer subunits than Ratcliffe (1984) did. Also, Karabinos and Pierce (2017) found no evidence to support the existence of isoclinal recumbent folds in the Dalton Formation that rests unconformably above Mesoproterozoic basement gneisses. Both Ratcliffe (1984) and Karabinos and Pierce (2017) observed that bedding in the Dalton Formation dips shallowly toward the northeast and southwest, and that gneissic foliation in the basement rocks commonly dips steeply and strikes approximately east-west (Fig. 11). We used a lidarderived texture-shade image and GMDE to calculate the strike and dip of bedding and gneissic foliation in this ~25 km2 area (Figs. 10 and 11) using three-point calculations. The results of orientation measurements made using GMDE and the texture-shade image agree well with field measurements reported by Ratcliffe (1984) and Karabinos and Pierce (2017), although the number of measurements is smaller than from the field-based studies (Fig. 11). To address the approximate linearity of three-point combinations, it was helpful to plot topographic contours on top of the texture-shade image to avoid selecting three points at similar elevations or along uniformly sloping surfaces.

The texture-shade image is helpful in showing areas with excellent bed-rock exposure compared to regions dominated by glacial cover. The contrast between the steeply dipping foliation in basement rocks and the shallowly dipping foliation in the clastic cover rocks is clear in the texture-shade image (Fig. 10B) and helps locate the trace of the unconformity between basement gneisses and the Dalton Formation. Further, the gentle anticlinal folding of cover rocks in the thrust sheet is visible; beds on the western limb dip to the west, and beds on the eastern limb dip to the east. Field work remains necessary, however, to ascertain what kind of planar feature is being measured with the texture-shade image and GMDE.

Jamaica, Vermont, is located on the eastern flank of the Green Mountain massif and west of the Chester dome (Karabinos, 1984; Ratcliffe, 1997). Deformed and metamorphosed Mesoproterozoic basement rocks are unconformably overlain by Ediacaran Tyson Formation (quartzite, marble, quartz-rich schist), the Hoosac Formation (pelitic schist, rare marble, amphibolite), and the Pinney Hollow Formation (pelitic schist, amphibolite). Acadian deformation and garnet-grade metamorphism affected this part of Vermont and obscured the effects of Taconic deformation (Doll et al., 1961; Karabinos, 1984; Ratcliffe et al., 2011). Karabinos (1984) mapped thrust faults and the Jamaica anticline in this ~40 km2 area (Fig. 12). It is rare to find well-preserved bedding except in the conglomerate, quartzite, and marble lithologies of the Tyson Formation (Fig. 12).

Figure 12A shows a texture-shade image for this area and the geologic interpretation of Karabinos (1984). The field work for this area was completed in 1978–1980, and the only available base maps were two 15-minute quadrangle maps published in the 1950s. As a result, the contacts shown in Figure 12A had to be transferred from an old base map onto a modern lidar-derived image manually. Every effort was made to match topographic features between the original base map and the texture-shade image.

Some contacts shown in Figure 12A correspond well with features in the texture-shade image. In particular, the basement-cored, overturned Jamaica anticline stands out in the texture-shade image because of the contrast in erosion between the resistant basement gneisses and the recessive marble in the thin overlying Tyson Formation quartzite-marble unit. The contact between mafic rocks in the Hoosac Formation (unit CZhtm, the Turkey Mountain Member) and the albite schist of the Hoosac Formation (unit CZhas) is also readily apparent. Portions of the Jamaica thrust (JT in Fig. 12A) correlate well with features in the texture-shade image, but the Cobb Brook thrust (CBT in Fig. 12A), as mapped, does not follow notable features in the texture-shade image. The Cobb Brook thrust separates two kinds of schist belonging to the Hoosac Formation, and the contrast between them may not be great enough to create distinctive topographic breaks. Two areas shown by white rectangles in Figure 12A are enlarged in Figures 12C and 12D. Each area contains a strike-and-dip value from a three-point calculation using GMDE, the texture-shade image, and the lidar DTM. The blue lines are contacts projected for 500 m from the strike-and-dip symbol on the surface of the DTM. The correlation between the projected contacts and the features in the texture-shade image is good for ~200 m from the point where the strike-and-dip symbol is placed but decreases at greater distances, suggesting that the contacts are not planar beyond 200 m.

The applications of texture shading to regions of metamorphic and igneous rocks tend to be more subtle than in areas of deformed sedimentary layers. There is, however, one class of structures that are exceptionally well imaged in metamorphic and igneous basement: distribution of brittle fractures. Such fractures would also be imaged on traditional hillshade images if the user positioned the artificial illumination direction at a high angle to the fracture strike. Although these features are seldom mapped, they can be of considerable importance given their role in channeling groundwater, the potential for induced seismicity in areas of fluid injection (e.g., Keranen et al., 2013), and, locally, their roles in slope instability.

Figure 13A depicts the fractures in metamorphic basement at the top of Mount Greylock in the Berkshire Hills of northwestern Massachusetts. A set of north-northeast–striking fractures with subsidiary, more northerly striking fractures are clearly defined by the narrow, linear, concave-upward depressions that appear dark on the image. South of the area shown, some of these structures appear to be moderately east-dipping reverse faults, one of the few places where one can distinguish between faults and joints. The geologic map of the same area (Ratcliffe et al., 1993; Fig. 13B) shows none of these fractures because their definition was not the intent of the mapping. Although the mappers locally measured the orientations of minor faults and brittle fractures, it is impossible to assess their potential role in groundwater or induced seismicity without also knowing something of the fractures’ lateral extent. The texture-shade image provides this information more completely than field work alone could.

Texture shading also enhances the imaging of surficial features including stream networks, river terraces, paleoshorelines, and glacial features. Figure 14 shows glacial features northeast of Ithaca, New York. The majority of the features in this image were initially recognized on lidar hillshade images, but texture shading enhances the images and improves the continuity, especially for the small recessional moraines, which have a maximum relief of 5 m above the surrounding terrain, and the “mega–glacial lineaments” (MGLs) with local reliefs of 2–10 m. In hillshading, one must use a low artificial illumination angle to image such subtle features, but in texture shading all that matters is that these linear features of molded glacial till have relatively sharp changes in slope with a convexupward geometry and are ubiquitous and systematically varying in their orientation throughout the image.

The well-resolved MGLs show how the ice flow was diverted by the Portage escarpment, a set of hills 150–250 m higher than the terrain to the north. The high terrain is underlain by Upper Devonian Sonyea Group siltstone and sandstone whereas the low terrain to the north is underlain by Upper Devonian Genesee Group shales. Although it can be difficult in this area to distinguish surficial and bedrock features on the textureshade image, Figure 14B shows the exquisite interplay between the two. There, a small recessional moraine drapes a bedrockcored bluff of nearly horizontal Genesee Group clearly visible despite the local housing development in the area. However, elsewhere in this area, especially along the walls of the north-south glaciated valleys, differentiating surficial and bedrock features can be quite difficult.

The light and dark banding in the texture-shade images appears to mimic stratification and compositional layering. As described in the Methods section, convex-up surfaces appear as light gray, and concave-up surfaces are dark (Fig. 2). Thus, texture shading is capable of imaging subtle changes in slope, referred to as “profile curvature” in many GIS programs, at the scale of resolution of the lidar DTM. Because of the scale invariance of texture shading, finer-scale changes in slope correspond to thinner, more erodible rock units at higher zoom levels. These slope variations could be due to a variety of natural and anthropogenic processes, but a common one would be the variation in resistance to weathering and erosion of the underlying layers. Even where mantled by thin to moderate soil development that completely covers the bedrock, the convex- and concave-upward changes in slope would persist. At the scale of a few meters typical of lidar DTMs, it is likely that texture shading is, most commonly, not imaging slope variations related to individual layers but packages of layers. As the thickness of the unconsolidated cover increases, we suspect that only proportionally thicker bedrock packages will appear as light and dark bands on a texture-shade image; this point, however, remains to be demonstrated with careful additional study.

Because texture shading is equally good at highlighting geomorphic features, such as river terrace edges and glacial moraines (case study near Ithaca, New York), as bedrock features, it can be difficult to distinguish variations produced by surface deposits from subtler variations produced by bedrock features. Many geomorphic surfaces are relatively horizontal and planar and thus appear similar to nearly horizontal bedding.

In metamorphic terrains, texture-shade images are valuable in assessing the continuity of geologic units over short distances, for calculating the orientation of the contacts between adjacent units, and for distinguishing dip slopes from escarpments. The unconformity present on Day Mountain in the Berkshire massif is apparent from the contrast between shallowly dipping bedding in the cover rocks and the steeply dipping gneissic foliation in the basement gneisses. Texture-shade images in conjunction with GMDE can help determine the orientation of foliation, but it is commonly difficult to obtain many such measurements. Because the topography appears to be highly correlated with geologic units, sets of three non-collinear points to use in the orientation calculation can be elusive. In addition, bedding is poorly preserved in most metamorphic rocks, so field work is essential to discover what kind of foliation is being imaged with lidar. Texture-shade images would be valuable in planning field work and in testing interpretations during mapping campaigns. They are also useful in assessing the validity of published maps. In these areas, the greatest practical utility of texture-shade images may be the mapping of brittle fracture sets, which are commonly not depicted on published maps because they cause little or no offset of map units and/or because the mappers’ primary focus is on the bedrock geology. To capture fracture sets of varying orientation using hillshading, one would require images with multiple illumination angles.

A field geologist measures the orientation of bedding based on a sample of the bedding surface as small as a few square centimeters, when the compass is placed directly on the outcrop, to perhaps a few square meters where sighting along a bedding plane is employed. Thus, even when the compass is in the hands of the most expert field geologist, stochastic variability of bedding surfaces is likely to introduce considerable scatter in orientation measurements. In contrast to compass measurements, orientations determined from three-point calculations on lidar data typically span tens to hundreds of meters, smoothing out more detailed variability in bedding.

Like an earlier, smaller study in southeastern Idaho (Allmendinger, 2020), the current study of the Nittany anticline bears this out. Comparing the π-diagram of 242 poles to lidar-based bedding orientations to 282 field geologist–measured bedding poles, the latter clearly has more scatter than the former (Fig. 9). A cylindrical best fit to each data set confirms the visual conclusion. Though the fold axis given by each data set is statistically indistinguishable from the other, the normalized magnitude of the smallest eigenvector for the lidar orientations is three times smaller than that for the field geologist measurements.

Of course, one person’s noise is another person’s signal: there are undoubtedly situations where meter-scale variability of bedding surfaces is the primary target of the research. For regional mapping at the scale of 1:24,000 or smaller, however, outcrop-scale bedding variability is noise. Such variability seldom helps to explain the distribution of rock units or aids in drawing cross sections. This does not mean, however, that outcrop measurements with compass are inferior or lack value. There are a few places in the current study where the field geologist identified small-scale structural or sedimentary structures not resolvable with the lidar data that may help in understanding the folding and faulting in the region.

Lidar data, enhanced by texture shading, brings a new dimension to imaging geology in areas of poor exposure and complements existing advanced techniques of image processing. Just as contact boundaries are defined by changes in lithologic properties, changes in slope reflect changes in erosional properties that reflect the underlying changes in lithology, even where the rocks do not crop out at all. The image is not perfect: units with uniform erodibility do not reveal well the underlying structure and orientation.

We imagine that producing and studying a texture-shade image will become a common first step before going to the field. The texture-shade image will provide observations and orientations to be field tested and will help the field geologist to identify and focus on critical areas requiring more attention. In the field, the geologist with a texture-shade image as a base map on their mobile device can use the image to assess the continuity of key rock units between sparse and isolated outcrops. The continuity of units in the texture-shade images limits the degrees of interpretational freedom otherwise available to the mapper in complex and poorly exposed regions. Texture-shade images, especially when draped over hillshaded DTMs, make for striking presentation graphics where, otherwise, a satellite image would just show dense tree cover with perhaps a few ledges poking out (e.g., Figs. 7, 10).

Texture shading is in some ways akin to an X-ray in that it reveals the underlying structure of a region based on subtle changes in how the skin (vegetation and soil) drapes the underlying structure. More resistant rock units produce convex-upward surfaces of the overlying regolith and soil whereas less-resistant layers produce concave-up surfaces. These subtle changes in surface slope exist even where bedrock outcrop is completely lacking. Texture-shade images robustly reveal both surficial and bedrock features because they are unaffected by artificial illumination direction bias and because of the inherent scale invariance, which allows resolution of small-scale variations corresponding to the layering of bedrock or subtle geomorphic features. We have concentrated on lidar DTMs here, but texture shading has proven use-ful even with lower-resolution DEMs such as the worldwide Advanced Land Observing Satellite (ALOS) 30 m DEM.

The technique described here will enable field geologists to extract much more bedrock detail in regions of sparse outcrop. In areas of deformed sed-imentary strata (e.g., Nittany anticline case study), substantial amounts of orientation data can be obtained in a short period of time. In metamorphosed terrains, the rock units may be less planar, but the continuity afforded by texture shading can prevent “flights of fantasy,” which sparse outcrops in otherwise highly vegetated terrains might otherwise permit. In all areas, but especially in metamorphic and igneous terrains, texture shading provides exquisite imaging of bedrock fracture sets, which can be critically important for any project involving subsurface fluid flow as well as natural and induced seismicity.

We thank Leland Brown, the originator of the texture-shading method, and Brett Casebolt for comments on an early version of this manuscript and for clarifications about the algorithm. Chelsea Scott provided a particularly helpful review of the current manuscript, and Ramon Arrowsmith pointed us toward key earlier references. Allmendinger is grateful to Teresa Jordan and Dan Karig for comments on lidar texture shading in the Ithaca area. We thank Patrick Kennelly, Ramon Arrowsmith, Mark A. Helper, and Associate Editor Barbara Tewksbury for their expert comments on earlier versions of the manuscript. Orientations and contacts mapped in this study will be made publicly available through the StraboSpot database (https://www.strabospot.org/). The authors declare no financial interest in any of the software mentioned.