).
Subsidiary windows include cross sections and spin windows. The
investigator can access attributes of the databases by selecting the
objects and requesting information.Geotouch: Software for Three and Four Dimensional GIS in the Earth Sciences
JONATHAN M. LEES
Department of Geological Sciences, UNC
PO Box 208109, New Haven CT 06520-8109
Abstract
A new program for exploratory data analysis in three and four dimensions is presented and described. Interactive communication between diverse datasets is stressed as the main gestalt of the Geotouch program. The primary kinds of data include point, vector, raster and wireframe data sets, in addition to specialized forms such as focal mechanisms and ellipsoidal information. The software includes methods for cutting cross sections at arbitrary angles, spinning objects in three space and animating time series of punctual data, such as hypocenter series and volcano eruptions. The program is written in POSIX compliant C using X-windows for Unix systems and has been ported to Linux. Free access, via the Internet, to executable binary, source code and documentation make this package an attractive alternative to expensive or unwieldy commercial options.
Key Words: Graphics, Interactive, Linux, Image, Map
Introduction
High speed computing with large memory capacity has opened new avenues for integrating large datasets of diverse origin in geology and geophysics. Difficulties still remain in the organization and database manipulation of these enormous quantities of data, and new challenges arise daily in the storage retrieval and interpretation of information that has several dimensions, including spatial and temporal variation. In this paper I illustrate and review the properties of a three- and four-dimensional GIS (Geographic Information System) system that is free, open and available over the Internet. In recent years the Geotouch program has evolved into a major program for GIS exploratory data analysis, primarily aimed at geophysical data, but also of great value to geologists. The program is written in X-windows for Unix operating systems, and works well under the popular Linux operating system.
The flexibility and broad range of tools available in Geotouch make it an attractive platform for many geological problem solving. It is easy to use and learn, the cost is unbeatable, and, since the source code is available, the program is expandable by anyone who has the skill to contribute. Other available GIS programs include the popular plotting program GMT and ESRI’s commercial package, ARCINFO, not to mention numerous image processing facilitates such as ERMAPPER. The problem with many of these alternatives is that the learning curve is especially steep and often the very implementations geoscientists prefer to apply require significant preparation of software and data. Geotouch provides a link between some of these programs for exploratory data analysis.
Geotouch is designed to be fast, so that it can
be run with the minimum number of input parameters. There are no memory
limitations, other than the limitations imposed by the local operating
system. The program accepts, either on the command line or via menu
driven dialogue boxes, names of files containing geographic objects and
other parameters, which describe the database. The main window is a map
view showing objects projected on the horizontal surface
(Figure 1).
Subsidiary windows include cross sections and spin windows. The
investigator can access attributes of the databases by selecting the
objects and requesting information.
General Software Overview
Target Region
At run time the Geotouch requires information telling it where on the globe to set up the origin. This is specified by an "origin’ file which has several different formats relating the location and dimensions of a target region on earth. Since often we do not know in advance where we need to investigate there is a facility for showing a very large target and then narrowing down the focus by visually selecting with a mouse drag. Target regions are specified as either Latitude-Longitude pairs or by specifying a point and dimensions in kilometers. Once the target is registered, all calculations in Geotouch are done in a Cartesian coordinate system using kilometers as the basic unit of reference. For this reason, the program is currently aimed primarily at local to regional scale analyses, although this is not specifically a restriction.
Projections
Map projections must also be specified at startup initialization. Several standard projections (coded from ) are provided, including: Albers Equal Area, Mercator, Miller, UTM-spherical, UTM-ellipsoidal, centered versions of these and, naturally, no projection for data that is not geographically coordinated. Installation and expansion of new projections is relatively simple, although it requires coding in C and re-compilation. All calculations in Geotouch are performed in units of kilometers after geographically coordinated data are projected onto a rectilinear space.
Geologic Maps
Geologic maps consist of several kinds of information, each coordinated to geographic coordinates (Figure 2 
and Figure 3
).
These include three basic elements: points, lines and polygons. In
Geotouch all data used in maps can be accessed on the screen by
selection. Each point in any object can be identified by its location
in the database, making editing the database relatively simple. There
is no facility in Geotouch for editing the database directly. Since all
files in Geotouch are stored in simple ASCII Text format, editing
functions are left to the user’s preferred text editor.
Point Data
Point data in Geotouch includes several types of specialized formats and one generic input. Since Geotouch was developed initially as a tool for investigating seismicity in three-dimensions, earthquake hypocenters and seismic stations are the two basic point datasets. Each of these datasets are represented as LAT-LON-Z (Depth) triplets, and they are organized in groups, such that different sets of points (e.g. hypocenters) can be associated simply. An advantage of this approach allows Geotouch to plot these in specialized ways: for example, hypocenters located with different velocity models can be linearly connected to see biases in relative re-location. This could also be applied to points changing position temporally.
The other mode of representing points in space is more generic. In this case a point has a name, a location and several parameters describing how it should be plotted, such as textual labeling, symbol style, and color. This mode of point representation is appropriate for plotting cities, geographic locations, or other textual information required in a presentation.
Finally, points may be plotted as part of the geologic maps, where sequences of points are associated by function. Examples include points plotted as part of a dashed fault trace that has no direct surface expression – lines of points that are associated geographically.
Line Data
Lines in Geotouch map files are tagged and registered as to their geologic significance. Lines or other geographic objects are given a category tag so they can be selected or treated differently depending on the needs of the investigator. Tags may include categories like roads, faults, geologic contacts, etc. (All Geotouch map objects have character string tags.) During an Geotouch session objects tagged can be selected or rejected for plotting purposes on the fly. The alphanumeric tags are generic and the user is free to determine the pattern of association as is appropriate to each situation. There is a suggested tagging scheme for those who do not wish to create their own.
Polygon Data
Polygon data are simply lines that are filled with a solid color or pattern. Since Geotouch is not a drawing program, per se, there are limitations in how patterns are used in plotting to the screen. In the Postscript output, however, there is unlimited availability of patterns and advanced users will certainly create any fill pattern that is required. Polygon filling is used mainly to plot geologic maps and a running legend is available for all polygon fill features can be plotted at any time.
Vector Data
In keeping with the three-dimensional character of Geotouch we have developed a facility for incorporating and manipulating 3D vector data. In basic form these are stored as simple LAT-LON-Z triplets organized in files for association. For example, three-dimensional well boreholes, orientation of anisotropic directions (arrows) , flow patterns and seismic raypaths can all be represented as 3D vectors in space. To add to the functionality of these data, users can specify colors for each segment of the vector individually, signifying, for example, lithology along the raypath. In the case of borehole data, dip information can be portrayed as small flags in space. When these are viewed in cross section the apparent direction of dip can easily be assessed. Another application of vector data might be displaying seismic waveform data on a map or cross-section. In this case a program would have to be written to convert the data from time to spatial coordinates, but this is a relatively simple conversion.
Color
Since Geotouch is an X-window program, color is handled via the normal X11 colormaps. By default these include mappings of 255 colors from the standard colormap to specific locations in memory. The program accesses the different colors by referring to the location of the pixel-map. When the program initiates, these colors are set and a corresponding mapping of GIS objects to specific colors is created. Geotouch allows the user to modify the default color map at run time by either changing the mapping of the colors or by re-defining which objects get colored differently.
When dealing with satellite images, the color map is handled differently. In this case one needs to determine the best colormap that fits the data. We have used an algorithm extracted from the Internet to accomplish this. In this case the satellite image is rendered by choosing the best 125 colors (Red-Green-Blue combinations) that fit the image and these are assigned on the fly to the color map. In this case, all the colors in the computer are owned by the application and flashing may occur.
Typically there is a fairly large discrepancy between how a plot on the phosphorescent screen appears versus how it looks on paper after plotting. This perennial problem exists at all levels in computer imaging. To address this issue I have made specialized Postscript color mappings to the colors chosen in the Geotouch by default. This is a fairly arbitrary choice, though, and can be modified in the Postscript files themselves for final publication.
Focal Mechanisms
A major effort has gone into the plotting of focal mechanisms in Geotouch (Figure 4
).
This is because geophysicists expend a large amount of effort in
dealing with these objects. We have restricted our plotting to simple
double couple solutions, although we may expand this to non-double
couple solutions at a later date. The basic data is stored as a series
of pairs of strike and dip, representing 2 focal planes, P and T axes
and the U and the V points. Geotouch reads in this data in a column
oriented format developed by the University of Washington, Seattle
(UW). External programs are provided to convert Strike, Dip and Rake
data to the UW format. In Geotouch focal mechanisms can be plotted as
normal beach balls, with an interpreted focal plan highlighted, or it
can be plotted with the focal plane alone. P and T axes appear as small
points on each sphere, and for some purposes, one can plot only the P-T
axes, or the strike and dip of the fault. All these different styles
can be plotted by toggling a style selector so one can quickly see the
data plotted in different ways.
The interactive facilities of Geotouch allow the user to select individual focal solutions and switch the primary and auxiliary planes. By interpreting large numbers of focal mechanisms oriented along faults where the primary focal plane lines up, one can accumulate a new data set relating the interaction of the fault planes and surface geology.
Often we have too many focal planes to deal with reasonably and a
summary is required. In this case a set of focal solutions is selected
and a summary is provided in the form of two plots: a stereonet of the
P-T axes projections and a ternary plot of all the focal mechanisms in
the selection (Figure 5
)
categorized by pure strike-slip, normal or reverse faulting . Finally,
ternary plots can be generated automatically by dividing the full
target into smaller rectangles and plotting a ternary plot for each
small box. This provides a way to see the how fault zones vary
spatially in a complex region where data sets are large and
heterogeneous (Figure 6
).
Wire Frame Objects
Wire frame
objects allow the user to import 3D surfaces for projection on
horizontal and vertical cross section. These three-dimensional objects
are stored as panel elements: a surface is broken down into small
surface pieces, each on covering a small part of the body (Figure 7
).
Panels are stored in memory and cross sections register the
intersection of the planes and projection surface. Three-dimensional
bodies might include models for a magma chamber under a volcano , or,
more simply, the three-dimensional orientation of a subsurface fault
plane. This is very useful for interpreting seismic reflection data and
surface geological observation of fault planes extrapolated to depth.
The features can be registered exactly, eliminating the potential for
misinterpretation.
Spinning facilities (rotation of objects about arbitrary axes) are currently available in many data management, statistical, and CAD programs. Geotouch has a spin module, optimized to provide specific information with which geologists are concerned: the final strike and dip of a particular view, after numerous rotations about arbitrary axes. Fault planes of aligned seismicity can be determined quickly and quantitatively. Point, vector and wire frame objects can each be selected and rotated to determine 3D relationships of imported properties.
Ellipsoids
One feature unique to Geotouch is its handling of ellipsoidal information. These were created to address the representation of error ellipsoids that arise in earthquake location routines, but its treatment in Geotouch is general applicable to any situation involving variance-covariance matrix representation of an error in 3-space. Earthquake locations are often represented as ellipsoids in computer generated hypocenter solution programs. In Geotouch these are stored as ellipsoids in space. In plotting these, they are projected onto the plane of reference, whether horizontal (map) or cross sectional. This allows the investigator to see the error bars of estimated parameters from any appropriate angle. As a practical example, the ellipsoid presentations were used extensively in displaying errors associated with multiplet earthquake clusters in the Coso geothermal field .
Raster Images
Raster images are a mainstay of GIS and remote sensing. In the geosciences these usually consist of a 2 dimensional gridded image where values in the third dimension represent spatial variation of a physical parameter such as elevation, bathymetry, gravity, aeromagnetic or seismic velocity variations to name a few. Herein lies one of the central limitations of conventional GIS software. Geophysical data is by its very nature three or even four-dimensional. Most investigators do not use a third dimension in GIS analyses, but for seismologists it is fundamental. In my own research we deal with 3D images derived from tomographic inversion . Geotouch provides a simple way to access the images visually and query parts of the image to determine precisely and interactively the relationship of anomalies to other independent datasets. Users can scroll through the model by layer so they can see how patterns change with depth or elevation. It is not the intent of this program to duplicate the image processing capabilities of popular software designated for that purpose. Rather, Geotouch provides a viewing platform for exploring the diverse relationships between the different geological parameters available. Some processing facilitates are available: for example the program will integrate a region covered by a raster grid so that volumes can be estimated. This is useful for determining volumes of mountains from topography, and for determining the relative weight of a specific property at subsurface.
When a cross section in Geotouch is requested, by dragging the mouse between two points of interest, the three-dimensional image is sliced along the vertical cross section and displayed along with other point and vector information stored in memory. All data within a specified volume are plotted onto the cross section by perpendicular projection. Wire frame bodies are sliced such that only the outlines in the cross sectional plane are shown and ellipsoids are outlined to show the covariance within a section. The three dimensional relationships of these objects are clarified and articulated in a way that can convey information to a broad audience. One features of cross sections in Geotouch is that if a contoured 2D field is in memory has units of depth, like a seismic horizon at depth, one can project this onto the cross section where they are sliced.
Fourth Dimension: Time
Finally we address the fourth dimension, time. Geotouch has a simple facility for animating earthquake hypocenters in time so apparent temporal relationships along faults or magma conduits can be illuminated. A similar analysis may be applied to volcanic eruptions or other time sequences of events. Two choices for animation are in place: in one mode events remain on the screen after they appear. In the other mode, events appear and later disappear according to their "influence" factor, calculated from event magnitudes, in the case of earthquakes. Another way of viewing temporal variations in 2D, static, raster images of changing properties is achieved by stacking 2D images into 3D raster files to create time slices through which one may scroll. While temporal analysis in Geotouch is still fairly primitive, the potential for this approach is clear.
Program Availability and Access
Geotouch is available for free on the Internet. Web page access can currently be found at ../lees.html and anonymous FTP access is accessible at ftp.unc.edu (most current) or ftp.iamg.org. There is a README file for setting up the program. The package includes executable codes for SUN Unix systems and Linux operating system, full source code and a set of demonstrations and tutorials. Over the last few years we have had at least one FTP access of Geotouch per day, on average. Suggestions for changes and improvements in the code are encouraged. Users can also add to and modify the source code as they see fit. Since official documentation is maintained on the Internet it is periodically upgraded and kept current.
Acknowledgements
The author thanks the Navy Geothermal Program (award #N68936-94-R-0139) for partial support of this project. I further acknowledge California Energy Co. Inc., Craig Nicholson, and Mike Huang for advice and program suggestions. Portions of the code were developed by Peilin Jia, Mark Lindner, and Neill Symons. Thanks to Mark Alvarez for originally suggesting I develop this code.
Bibliography
Frohlich, C., 1992. Triangle diagrams: ternary graphs to display similarity and diversity of earthquake focal mechanisms. Physics of the Earth and Planetary Interiors 75, 193-198.Frohlich, C., Apperson, K. D., 1992. Earthquake focal mechanisms, moment tensors, and the consistency of seismic activity near plate boundaries. Tectonics 11, 279-296.
Lees, J. M., 1992. The magma system of Mount St. Helens: Non-linear high resolution P-wave tomography. Journal of Volcanology and Geothermal Research 53(1-4), 103-116.
Lees, J. M., 1998. Multiplet analysis at Coso Geothermal. Bulletin of the Seismological Society of America 88(5), 1127-1143.
Lees, J. M., 1990. Tomographic P-wave velocity images of the Loma Prieta earthquake asperity. Geophysical Research Letters 17, 1433-1436.
Lees, J. M., 1995. Geotouch: Three-dimensional GIS for geology and geophysics. Seismological Research Letters 66(4), 33-37.
Lees, J. M., Crosson, R. S., 1990. Tomographic imaging of local earthquake delay times for 3-D velocity variation in western Washington. Journal of Geophysical Research 95(B4), 4763-4776.
Lees, J. M., Lindley, G. T., 1994. Three-dimensional attenuation tomography at Loma Prieta: Inverting t* for Q. Journal of Geophysical Research 99(B4), 6843-6863.
Lees, J. M., Malin, P. E., 1990. Tomographic images of P-Wave velocity variation at Parkfield, California. Journal of Geophysical Research 95(B13), 21,793-21,804.
Lees, J. M., Nicholson, C., 1993. Three-dimensional tomography of the 1992 Southern California sequence: Constraints on dynamic earthquake ruptures? Geology 21(5), 385-480.
Lees, J. M., Vandecar, J. C., 1991. Seismic tomography constrained by Bouguer gravity anomalies: Applications in Western Washington. Pure and Applied Geophysics 135(1), 31-52.
Lees, J. M., Wu, H., 1999. P-wave anisotropy, stress, and crack distribution at Coso Geothermal Field, California. Journal of Geophysical Research, 104(8), 17,955-17,973.
Pallister, J. S., R. P. Hoblitt, D. R. Crandell, D. R. Mullineaux, 1992. Mount St. Helens a decade after the 1980 eruptions: magmatic models, chemical cycles, and a revised hazards assessment. Bulletin of Volcanology: 54, 126-146.
Snyder, J. P., 1987. Map Projections- a working manual. USGS Professional Paper, 383 pp.
Wessel, P., Smith, W. H. F., 1991. Free software helps map and display data. Eos (Transactions, American Geophysical Union) 72(41), 445-446.
Figure Captions
Figure 1. Screen dump of a typical Geotouch session. The larger window in the back is the main window showing a map view of the target region. The smaller window on the left is cross section 1 as seen geographically in the Map view. The third window is a spin window showing selected elements (hypocenters and vectors) rotating in 3D.
Figure 2: A) Topographic map of Mount Rainier and Mount St. Helens regions, Washington state. Rectangular boxes (1, 2) outline volumes projected in B and C. Gray circles are earthquakes and triangles are stations from the Western Washington seismic network. B) Vertical cross section of Mount Rainier. Red and blue colors are low and high velocity regions presented as percent perturbation from the background one-dimensional models. Small circles are projected seismicity. A line drawing of the proposed magma system below Rainier's edifice is an interpretation based on tomographic analysis and similarity with Mount St. Helens inversion. C) Mount St. Helens cross section. Color codes as in B. Focal mechanisms are front projections of seismicity, primarily located in the stoping region, color coded according to rake. Line drawing is an interpretation of the plumbing system of Mount St. Helens based primarily on seismicity and geochemical analysis (Pallister, et al., 1992)
Figure 3: Magnetic anomalies (image), Bouguer gravity (contours) and seismicity in South Western Washington. The southern Washington Cascades conductor (SWCC), outlined for reference, correlates with relatively flat gravity and magnetic anomalies. Note the correlation of seismicity in the WRSZ and the SHZ with anomalous negative magnetic lineations (blue).
Figure 4: Focal Mechanisms in Geotouch can be plotted in 6 different ways: A) simple beach ball; B) Beach Ball with P and T axes marked; C) Beachball with focal plane delineated; D) just the focal plane with the slip vector marked; E) Just the P-T axes; E) the strike and the dip of the fault plane. The user can toggle between these choices interactively.
Figure 5: Ternary plots in Geotouch follow the style of Frohlich . Focal mechanism data is divided up into categories based on rake angle. This plot is useful when hundreds of focal mechanisms are being analyzed and spatial patterns are difficult to ascertain.
Figure 6: Ternary plots of focal mechanisms in the Aleutian-Kamchatka subduction zones. The plots show the trend of the mechanism changing from thrust to strike slip in the Bering fault in the far western Aleutians. When the subduction zone trends south in Kamchatka, the earthquakes are again thrust events.
Figure 7: Wire frame representation of a sinusoidal surface with hypothetical boreholes sampling the surface. The figure is shown as an example of the spinning routine in Geotouch.