## Abstract

GMDE is a program, available on the desktop for Apple Macintosh, Microsoft Windows, and Linux platforms and in a mobile version for iOS, that enables geologists to extract quantitative structural information from geologic maps and satellite images. The program facilitates the digitizing of strikes and dips or calculating them from three-point problems, calculation of stratigraphic map thickness, determination of piercing points on faults, and the construction of down-plunge projections and vertical cross sections with projected apparent dips, contacts, and cylindrical folds. The program also enables the automatic plotting of planar contacts across topography based on orientation calculated from three clicked points, which can be carried out in the field for immediate hypothesis testing. Error propagation is built into many of the calculations. Maps and satellite images require no projection or datum information, just four points with known latitude and longitude information. Alternatively, the user can enter the map or image in MBTiles format. Users can easily extract *X*-*Y*-*Z* data for any clicked or calculated point or polygon, enabling them to make their own calculations.

## INTRODUCTION

Geologic maps are seen in some quarters as arcane graphical devices that are rapidly losing their significance in a world awash with digital data, easily accessible online with the click of a mouse. This impression is inadvertently reinforced in the typical undergraduate structural geology laboratory where students are taught a variety of seemingly unconnected graphical constructions. We now have the tools to do twenty-first century map interpretation. All of the calculations typically taught in structural geology and performed by professional geologists can be carried out precisely and accurately with simple vector calculations and matrix transformations from linear algebra. In the past, however, it was so tedious to extract *X*-*Y*-*Z* data from paper copies of geologic maps that few people took advantage of the quantitative information contained in these maps. Those data can be extracted from full-fledged GIS programs, but the learning curve of those programs is typically steep and the cost high. While many of those systems can be scripted to carry out specialized structural geology calculations, the functions required are seldom built in.

Three trends now make it feasible to extract information from a geologic map quickly and easily: (1) the widespread availability of free raster images of geologic maps from national and state geological surveys such as the U.S. Geological Survey (USGS) and the Geological Survey of Canada (GSC); (2) the ability to access the elevation of most points on Earth with a simple query over the internet or from free digital elevation models (DEMs); and (3) easy access to high-resolution satellite imagery from Google and other providers.

This paper describes a computer program called GMDE (for Geologic Map Data Extractor; Fig. 1)^{1} that is available free of charge for Apple Macintosh, Microsoft Windows, and Linux platforms. A companion program, GMDE Mobile for Apple iPads (but not iPhones), is freely available for use in the field. GMDE is written in Xojo, a modern compiled, object-oriented language with a Basic-like syntax that enables rapid development and multi-platform deployment. The menu and window–driven program is highly interactive, making it easy for those with little prior experience to make accurate measurements on geologic maps and satellite images, digitize contacts and strikes and dips already on the map, prepare topographic profiles, construct down-plunge projections, determine orientations from three-point problems, calculate piercing points on a fault, project planar contacts across topography, and more. If the user specifies uncertainties on horizontal and vertical position, the program propagates those errors through to the calculated values. In addition to allowing visualization on satellite imagery within the program, the program can also export KML files of the features digitized and plot three-dimensional orientation symbols in Google Earth (Blenkinsop, 2012).

## PREPARING THE MAP AND TOPOGRAPHY

### Base Maps

GMDE can use base maps in two general formats. The desktop version can read single raster images at full scale and then display them at different zoom percentages. It can also read maps in MBTiles format (Fischer et al., 2018). The mobile version can only ready MBTiles format because of the small amounts of random-access memory (RAM) in most mobile devices.

#### Raster Images

GMDE can read geological maps and satellite images in several raster image formats—JPG, TIFF, PNG—as well as Acrobat PDF format. Many published maps contain incomplete datum and projection information. Thus, the program does not expect there to be any geocoding information attached to existing files nor does it use georeferencing if it does exist. Instead, immediately after the program reads in the file, the user clicks on four points of known longitude and latitude to provide the necessary georeferencing. Georeferencing only needs to be done once because the information is automatically saved as a text file in the same directory as the map or image and is read automatically the next time the user opens the map. Maps and images are not required to be oriented with north at the top; the transformation that GMDE performs is also valid for rotated maps.

*x*and

*y*are the column and row number of each pixel in the image, and

*x*′ and

*y*′ are the Universal Transverse Mercator (UTM) Eastings and Northings, respectively, that most closely match the image; the latter are calculated from the longitude and latitude coordinates of the four clicked points. The transformation (

*A*,

*B*;

*D*,

*E*) and translations (

*C*,

*F*) are not known

*a priori*. To determine them, the program does a least-squares inversion by rearranging the above equation as follows and solving for the column vector of unknowns on the right-hand side:

*A*–

*F*, but GMDE uses four points, which allows us to estimate uncertainties. Those uncertainties are expressed as a χ

^{2}value; a moderately good value would be close to the number of degrees of freedom (Press et al., 1986), which in the case of UTM positions at four known points would be two. The user can assess the goodness of fit of a UTM projection to their image by examining the χ

^{2}value (Fig. 2). Of course, the starting image may be in any projection, which is commonly not specified fully, if at all, but if the area covered by the image is small enough, the error of assuming a UTM projection is small, and structural calculations can be made accurately.

Most elevations, either retrieved from an internet elevation server or in downloaded DEMs (next section), are referenced to the World Geodetic System 1984 (WGS84) geodetic datum (although some older DEMs in the United States are referenced to the North American Datum of 1927 [NAD27]). To use elevations in calculations, the difference between the datum used by an elevation server and that of the image can be critical. If the datum of the image is known, then the user can select it from a drop-down menu. GMDE then calculates the difference using the multiple linear regression equations (National Geospatial-Intelligence Agency, 2014). This is a critical step because the difference can be significant. For example, in the part of southeastern Idaho (USA) used in some examples below, the difference between the North American Datum of 1983 (NAD83) and NAD27 is >60 m (Fig. 2).

Unfortunately, the datum of a map or satellite image is commonly not specified at all or not in sufficient detail to allow precise conversions between it and WGS84. To circumvent this problem, GMDE offers a simple procedure to determine a constant offset between the image and WGS84. The user finds a point or feature—e.g., a road intersection, stream junction, topographic summit, or building—that is easily and precisely identifiable on both the image and in Google Earth (which uses WGS84). The latitude and longitude of the feature in Google Earth is recorded. Then, in GMDE the user clicks on the same feature, and in the ensuing dialog box, enters the coordinates recorded from Google Earth. GMDE then calculates the linear offsets in *X* (east) and *Y* (north). After setting the coordinate system offset, the program reports all longitude and latitude values in WGS84 coordinates, meaning that when the user points to a corner point on the map, GMDE reports the position with a WGS84 datum, which would probably be slightly different from the value entered for that point.

All of these approximations are acceptable as long as the image area covered is not too large. How big is “too large”? There is no definitive answer to that question. USGS 1:24,000-scale quadrangle maps covering ∼1500 km^{2} and GSC 1:50,000-scale quadrangle maps covering ∼5500 km^{2} work very well with χ^{2} < 3 even though neither are typically found in a UTM projection. Intermediate known latitude-longitude points on these maps are accurately picked in GMDE within the ±2 pixel uncertainty of the mouse cursor. Satellite images saved from Google Earth may have to be reprojected in a UTM coordinate system using a third-party program prior to reading them into GMDE.

If the absolute location of points on the map is not of concern, the user can bypass the process of georeferencing because the program also keeps track of the local *X*-*Y* values based on pixels in the image, which have their origin at the lower left corner of the image. These are the numbers that appear in the datum information area (Fig. 1, upper right corner) for east and north coordinates, rather than the longer and more cumbersome UTM coordinates. In such cases, the scale can be set by dragging a line along the scale bar of the map, or between any two points of known distance, and then entering the distance and the units. This works well if the user is entering elevations by reading topographic contours on a scanned map.

#### MBTiles

^{level}pixels. Thus, the ground resolution is:

*radius*is the radius of the Earth in meters at the equator (6,378,137 m). The scale of the tile is then:

*screen DPI*is the resolution of the mobile device screen in dots (pixels) per inch (DPI).

Internet satellite images use a pseudo-Mercator (i.e., Spherical Mercator) projection where, for reasons of calculation speed, the Earth is assumed to be perfectly spherical rather than ellipsoidal, resulting in a 0.33% scale distortion in the north direction (Schwartz, 2018).

To prepare a tiled base map from existing raster images, a common solution is via MBTiles, an open format originally proposed by MapBox.com (Fischer et al., 2018) that uses the same math and projection as just described. An MBTiles file is actually an SQLite database in which the tiles at different resolutions are actually SQLite “blobs”. For a typical 7.5-minute quadrangle area, an MBTiles file with zoom level between 10 and 17 may contain >10,000 individual tiles. At level 17, the ground resolution is ∼1.2 m/pixel and at 96 pixels/inch display resolution, the scale of the image is ∼1:4500.

### Topography

Many mapping systems provide *X* and *Y* coordinates, but to calculate things of interest to the geologist, one also needs *Z*, the elevation. GMDE desktop can retrieve point and profile elevations from internet elevation servers in real time, or read DEMs downloaded from the internet, which is necessary for offline use and for a few intensive calculations.

#### Internet Elevation Servers

If the user is connected to the internet and the image has been georeferenced, GMDE can automatically obtain the elevation of any clicked point from elevation servers, speeding up the entry of an elevation coordinate, especially for novice users. At present, four servers are implemented: for elevations in the United States and Canada, servers provided by the USGS and the GSC (both ∼10 m resolution) can be used; alternatively, for elevations anywhere in the world (including North America), the MapQuest or GeoNames servers, based on Shuttle Radar Topography Mission (SRTM) mapping (generally with a 30 m resolution), can be used. All four servers return the elevation value at any given longitude and latitude interpolated on a grid using the WGS84 datum; errors are typically not reported. Each server returns data a bit differently, and some can return a complete topographic profile whereas others, only single point elevations.

It is instructive for the user to see the variation in elevation values returned by different servers. The given elevation at a single point can vary by >10 m depending on the server (or topographic contours) used. This uncertainty can be important when estimating errors in elevation.

#### Offline Digital Elevation Models

Digital elevation models come in a bewildering array of formats; GMDE accepts only a single but widely popular format: binary GridFloat and its text equivalent, ASCII Grid. This type of DEM is one of those offered by the USGS and other government agencies for direct download and can be written by most popular GIS programs such as ArcGIS and Global Mapper. A GridFloat DEM consists of three files: (1) a small header file (.hdr) that has information about the number of rows and columns in the DEM, the grid spacing, and other file attribute data; (2) a datum and projection file (.prj); and (3) a large binary file (.flt) containing the actual elevation data. The grid comes in two basic flavors: (1) a grid spacing and bottom left corner specified in decimal degrees, or (2) a grid spacing and bottom left corner specified in meters. GMDE can read both formats, although for the latter type, the grid must use a standard UTM central meridian. Finally, DEMs can be produced using a variety of map datums. GMDE can recognize several different types of map datums specified in the .prj file and accounts for the difference with the base-map datum. If the application does not recognize the datum, it will default to WGS84, which might not be correct.

Digital elevation model files tend to be quite large. For example, the USGS releases its ^{1}⁄_{3}-arc-second DEMs as 1° × 1° tiles, ∼0.5 Gb in size, that generally cover an area much larger than the user needs. GMDE has no problem with oversized files because it does not read the DEM into memory; instead, the program does a direct read from the position in the binary file containing the needed elevation. The desktop version of GMDE can trim the DEM to the map area and save a new smaller GridFloat DEM to disk if needed.

## GEOLOGICAL CALCULATIONS

Many basic structural geology calculations can be carried out computationally (Allmendinger et al., 2012) by converting orientations to direction cosines in a Cartesian north-east-down (NED) coordinate system and then using vector operations such as dot, cross, and dyad products (Table 1). Map data are usually referenced to an east-north-up coordinate system. Several calculations involve the transformation into a coordinate system fixed to the structure, as described below. In these cases, the transformation matrix is composed of the direction cosines of the new coordinate axes in the old coordinate system.

*X*-

*Y*-

*Z*coordinates) has uncertainties, the calculated quantities also have uncertainties. GMDE uses the standard Gaussian error propagation equation to determine the errors on the derived parameters (strike and dip, thickness, trend and plunge, etc.):

*d*is the propagated error on

*d*, which is calculated from measured parameters α, β, and γ (plus any others), and δα, δβ, and δγ are the errors on those parameters.

### The Three (or More)–Point Problem

The orientation of a plane crossing topography (or pierced in multiple boreholes) can be determined from three non-collinear points on the plane. The pole to a plane can be calculated using a vector cross product (Allmendinger et al., 2012) or with the approach outlined in the next section. Calculation has the advantage that the errors can be propagated (Allmendinger and Judge, 2013), which is especially useful for identifying degenerate cases where the points are nearly collinear.

*v*is the vector formed by subtracting the centroid from the position vector (Equation 6). Finally, the eigenvalues and eigenvectors of the covariance matrix, σ, are determined. The third eigenvector (corresponding to the smallest eigenvalue) represents the pole to the best-fitting plane, and the square root of the smallest eigenvalue gives the standard deviation of the distance of each point from the best-fitting plane (Fig. 3).

### Map Stratigraphic Thickness

*t*) and base (

*b*) of the formation. This formula is derived by performing two coordinate transformations from data in a map (ENU) coordinate system to NED then to a coordinate system defined by the bedding where axis

**X**″ is parallel to the strike,

**Y**″ parallel to the true dip direction, and

**Z**″ in the direction of the downward-pointing pole to bedding (Fig. 4). In this coordinate system, indicated by double primes in Equation 8, the thickness measured perpendicular to bedding is given by the difference in the third coordinate of position vectors (parallel to the pole) of the top (

*t*″) and base (

_{z}*b*″). The detailed derivation of this formula and propagation of errors are given in Allmendinger and Judge (2013). In areas of high relief, difficult access, and few orientation measurements made in the field, measuring three points on one formation surface and a fourth point on the other bounding surface can provide substantial quantitative thickness and orientation data in a short period of time.

_{z}### Piercing Point Calculations

Fault slip is generally determined from the offset of a linear feature displaced across a fault. The intersection of that linear feature with the planar fault surface results in two piercing points. This problem is geometrically identical to the question of an exploration geologist who wants to drill a hole into an inclined plane, e.g., a mineralized fault zone, at depth.

**p2**, and the orientation of the fault plane (represented by the pole to the plane,), and (2) two points, one on the hanging wall (HW) and one on the footwall (FW) of the fault plane, and the orientation(s) of the offset line at those two points (

**p1**and

_{HW}**p1**; Fig. 5). The unknown piercing points,

_{FW}**p**and

_{HW}**p**, must lie on the line with known trend and plunge at the surface,

_{FW}**p1**:

_{HW}and p

_{FW}) and

*u*is the distance from

**p1**to the plate.

**p1**must likewise lie on a vector in the fault plane, which is perpendicular to the pole to the plane at the known point

**p2**. Because the dot product of two perpendicular vectors is zero, we can write:

**p**and

*u*. The latter is a scalar distance which is given by:

*u*into Equation 9, we can calculate the scalar components (i.e., coordinates) of the position vector,

**p**:

We repeat the process for the offset linear feature on the other side of the fault and then calculate the distance between the two offset points to get the fault displacement.

### Fold Axis and Down-Plunge Projection

*Operations*>

*Calculate Cylindrical Fold Axis*. That fold axis corresponds to the least principal axis of the orientation tensor given by:

^{th}pole in the data set.

*a*

_{ij}, of the down-plunge projection (from the ENU coordinate system of the digitized bedding to the fold coordinate system) is a function of the fold axis trend and plunge (Allmendinger et al., 2012):

The down-plunge projection plot window in GMDE (Fig. 6) shows all of the checked contact polygons projected into the profile plane of the fold. The user can adjust the fold axis trend and plunge incrementally and see the fold reprojected in real time. A vertex selected in either the map or the down-plunge view causes the same vertex to be selected in the other view, and thus anomalies identified in the down-plunge view can be related to the map. A button in the down-plunge window allows the user to save the profile view as a high-resolution scaleable vector graphics (SVG) file for entry into any modern vector graphics program or web browser.

The profile plane of a cylindrical fold is most commonly not a vertical plane but inclined at an angle equal to 90°–*plunge*. GMDE can also calculate the elevation of a cylindrically folded contact at depth, vertically beneath a point at the surface (*Operations* > *Depth to Folded Surface*). Likewise, it can also project the digitized vertices of contacts onto vertical topographic cross sections.

### Projecting a Plane across a Digital Elevation Model

*A*,

*B*,

*C*, and I are the coefficients in the standard equations for a plane; the first three define a vector which is the pole to the plane. If one knows any three, non-collinear points on the plane, the values of

*A*,

*B*, and

*C*can be calculated from the following determinants:

In the case of our map, the three points, **Q**, **R**, and **S**, are specified in UTM coordinates and represent three known points along an exposed formation contact. Alternatively, one can specify a strike and dip at one location on the plane and, from that, calculate the position of any two other points on the plane but not necessarily at the surface. GMDE arbitrarily locates these two additional points at a user-specifiable distance (by default, 500 m) away from the point with the orientation, at rakes of 45° and 135°, either up dip or down dip from the measured point. This is done via a coordinate transformation to the plane coordinate system (strike-dip-pole) and then a reverse transformation back to geographic coordinates.

Once all three points are specified or calculated for the plane, Equation 16 is calculated with the actual elevation values at each node in the DEM grid. If a point on the surface is on the same side as the pole vector [*A*, *B*, *C*], it yields a positive number when substituted into Equation 15, and if on the opposite side of the plane, it yields a negative value (Fig. 7). GMDE then contours the zero value using a custom implementation of the contouring algorithm CONREC (Bourke, 1987), producing the actual surface trace through the entire portion of the DEM grid. A subroutine by Allmendinger joins the individual contour segments into continuous polygons and thus allows as many intersections of the plane with the land surface as necessary.

## APPLICATIONS

GMDE is equally useful to the professional geologist in the field or office, and in a classroom setting. A few usage case scenarios are outlined below, but the program is capable of considerably more than the applications described here. All data collected with GMDE Mobile can be uploaded directly from the application to the user’s account at the StraboSpot website, https://strabospot.org.

### Structural Geology and Map Reading Classroom Activities

Many of the calculations performed by GMDE are taught in standard introductory structural geology laboratory exercises. GMDE makes it easy for students to capture the actual strike and dip values on a published geologic map, which can be digitized simply by dragging a line parallel to the strike line with the dip either entered manually or using the length of the dragged line. Those values can then be used in exercises such as calculating mean vectors or cylindrical best fits to poles to bedding. Because dragging any line on the map displays both the distance and the azimuth of the line, students can make their measurements easily and precisely directly on the screen rather than fussing with protractors and rulers. Topographic profiles with apparent dips projected parallel to the fold axis trend can be constructed automatically so that students can focus on interpreting the geology rather than the tedium of profile construction.

The professor who wants to challenge their students a bit more can ask them to calculate the values using the same algorithms that GMDE uses. The program makes it easy to capture the coordinates of the position vectors, which can then be pasted into a spreadsheet or more sophisticated computing environment to do the calculations. Special menu options such as *Edit* > *Copy Special* > *All Coordinates (NED)* will copy to the system clipboard, with a single menu command, the three components of the three position vectors that go into a three-point problem (nine numbers in all), formatted so that they can be incorporated as a matrix into a spreadsheet with a single paste command.

In my own class, several of these methods constitute an end-of-semester exercise (provided as Supplemental Material^{2} to this paper) using map data from the Canadian Rocky Mountains (Price and Mountjoy, 1978). In a small area (Fig. 8), the student must determine a number of strikes and dips using three-point problems, determine the fold axis, extract map thicknesses, calculate piercing points on a reverse fault, interpret map patterns to see how many folds there are and whether the thrusting is in or out of sequence, and draw a cross section. The students must calculate all of the values from the position vectors that GMDE provides, but the exercise is also useful for those professors who want to treat the program as more of a black box so their students can focus solely on the geological interpretation.

### Evaluation of Published Maps

A professional geologist in academia, commercial enterprise, or government is far more likely to be using and evaluating published geologic maps than making their own maps. GMDE can be of substantial assistance in extracting quantitative information from the map. Strikes and dips can be digitized quickly and cylindrical best fits determined in the program without having to go to a stereonet program. GMDE is a useful starting point for cross-section construction, especially given that GMDE can use a cylindrical fold model to project surface contacts up or down plunge onto the vertical plane of a cross section of any orientation. Depths to a folded surface from any point on the surface can quickly be calculated. Piercing points can be used not only to determine fault slip but also to determine how far one must drill on a particular azimuth and plunge in order to intersect a plane of interest.

GMDE can also assist in identifying problem areas in the mapping itself. One such example is shown in Figure 9. The published strike and dip on the map (Albee and Cullins, 1975) is inconsistent with the orientation from a three-point problem, and projection of the contact using the geologist-measured strike and dip further emphasizes the problematic nature of this part of the map. Contacts constructed based on an orientation determined from a three-point problem and projection across the DEM are substantially more accurate than those drawn by the mapper, and three-point orientations describe map-scale structures more accurately than outcrop measurements (Fig. 10).

### Real-Time Hypothesis Testing in the Field

The ability to test predictions and hypotheses in the field is now enabled by a generation of smart mobile devices with computing power similar to that of their desktop counterparts. Projecting planar contacts across topography based on orientations can be carried out in the field with the iPad version of GMDE or using the Windows version on a Microsoft Surface tablet.

In the area of southeastern Idaho described previously (Albee and Cullins, 1975), we have been able to show that some stratigraphic contacts are almost perfectly planar across >3 km of strike distance. Where the contact deviates from that planar geometry, it can help the field geologist identify previously unknown faults, subtle fold hinges, or even abrupt stratigraphic variations. One such example comes from projecting the contacts of two distinctive sublithographic limestone units, the Cretaceous Peterson and Draney Limestones in the Poker Peak 1:24,000-scale quadrangle (Albee and Cullins, 1975) (Fig. 11). Projection of the contacts of these two units, based on four three-point problems, and the ensuing systematic offset identifies the existence of an east-northeast–striking fault that does not crop out anywhere. In addition, the fault displayed immediately to the south on the map does not exist based on the continuity of the contacts.

## CONCLUSIONS

Geologic maps are quantitative scientific documents, but their utility has been hampered by traditional graphical methods of extracting that information. GMDE is a versatile program for extracting quantitative information from geologic maps quickly and easily, providing substantial utility for both professional geologists and students. In the office, the program can be used to analyze published maps and satellite images, and in the field, the mobile version of the app can make predictions that are testable while the geologist is still in the field. The ability to retrieve orientations from three-point problems and project planar contacts across DEMs will result in more accurate geologic maps. The feature that makes these calculations possible is immediate access to online or offline DEMs to provide complete *x*-*y*-*z* coordinates of position vectors used in the calculations. Now that structural geologists have easy access to this information, the discipline must become more quantitative if it is to enter the era of computational field geology.

## ACKNOWLEDGMENTS

I am grateful to Nathan Niemi, Sara Carena, Paul Karabinos, and Néstor Cardozo for beta testing various versions of GMDE and suffering through many crashes and unexpected behaviors. Néstor Cardozo and Adam Cawood provided expert reviews of the manuscript, and I thank Associate Editor Claire Bond for guidance with the revisions. I was taught field geology and reverence for geologic maps by three of the best: USGS geologists Steve Oriel, Gene Boudette, and Max Crittenden, to whom I am indebted. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. The author has no competing or financial interests.

^{1}GMDE and GMDE Mobile are completely free and come without any restrictions for all users. The desktop versions of the program are available for Microsoft Windows, Linux, and Apple Mac OS for 64-bit systems. The Windows and Linux versions are also available in 32-bit versions. All can be downloaded freely and anonymously at: http://www.geo.cornell.edu/geology/faculty/RWA/programs/strikedipthickness.html. GMDE Mobile is available for iPads (Apple iOS operating system) and can be downloaded for free from the iOS App Store.

^{2}Supplemental Material. Contains a sample classroom exercise using a map from the Canadian Rocky Mountains that can be done with GMDE. Please visit https://doi.org/10.1130/GEOS.S.12861050 to access the supplemental material, and contact editing@geosociety.org with any questions.