Allen H. Cogbill 398 Venado Street Los Alamos, NM 87544-2435 505-661-9754 E-Mail: acogbill@geopotential.com
Copyright 1997-2000, Allen H. Cogbill, all rights reserved. The program InnerTC© and DemRead© are licensed, not sold, and may not be redistributed without express written consent of the author.
Most computer users realize that only the blandest warranties accompany software. This software is no different from any other in that respect. Accordingly, there are no warranties of merchantability or fitness for a particular purpose, nor will the author be responsible for any damages (incidental, consequential, or otherwise) from the use of this program.
A program has been written to use the Digital Elevation Models produced by the U. S. Geological Survey in order to calculate inner-zone terrain corrections for gravity stations. The program provides inner-zone corrections that are many times more accurate than the corrections calculated using traditional template and chart methods. Documentation is provided herein to enable the user of the program to quickly learn how to run the program, as well as providing sufficient information to understand the assumptions involved and the inherent limitations of the procedures.
This program is designed to produce so-called "inner-zone" terrain corrections. That is, the intent of the program is to calculate the terrain correction from very close to the gravity station to a distance at which previously-developed terrain correction programs, which typically use much coarser terrain data, can be used. A typical use of the program would be to calculate the terrain correction starting a few meters from a gravity station to a distance of a few kilometers from the station. At greater distances, terrain correction programs such as that of Plouff (1966) should be used. Note that InnerTC does not provide reductions for gravity measurements; it only calculates the inner-zone terrain corrections. A future development is to incorporate the program into a more general program for complete terrain corrections and gravity reductions.
Reports of bugs in the program and suggestions for improvements should be sent to the author. The preferred method for reporting bugs or providing suggestions is by electronic mail.
The original program was written in 1988 as an outgrowth of a correction code that had been developed by me in the early 1980's. The original correction codes used irregularly sampled elevation data, often obtained by simply digitizing the contours of a topographic map. These irregularly sampled terrain data were used as input to a very accurate bivariate interpolation scheme based upon the multiquadric equation approach espoused by Hardy (1971) and later used successfully by Krohn (1976) to correct gravity measurements for terrain effects.
The advent of the U. S. Geological Survey's 30-meter Digital Elevation Models made possible the extension of the previously-developed terrain correction methods to the regularly-sampled digital elevation models. The original inner-zone correction procedures were in essence identical to the procedures used on irregularly sampled terrain models. As before, multiquadric equations were used for interpolating the elevations, at least in the immediate proximity of a gravity station.
The use of multiquadric equations to interpolate the elevation data, while highly accurate, was exceedingly expensive in both machine time and, more importantly, memory requirements. The amount of time and memory required to set up, solve, and interpolate the elevation data depends upon N×N, where N is the total number of elevation samples that form the basis of the interpolation. For this reason, I investigated other methods of interpolating the terrain data, especially in the immediate proximity of a gravity station. I settled on the triangulation and interpolation methods made available in Renka (1984) which were accurate, fairly rapid, and general enough to accomodate situations in which some of the terrain data near a station might be missing or sparse, due perhaps to a station being located near the edge of the available DEM coverage. The triangulation methods were incorporated into the correction procedures in mid-1990, not long after publication of the method in Geophysics (Cogbill, 1990). The terrain correction code was little changed from 1990 until late 1994. At that time a number of changes were made: the calculation of distances in the procedures was made more accurate by taking advantage of the fact that the DEM nodes are always of integral value, input to the program was modified to permit all parameters to be entered on the command-line, and considerable flexibility in data file structure was added.
The U. S. Geological Survey (USGS) is engaged in a program of developing Digital Elevation Models (DEM's) for 7.5' topographic quadrangles. These models consist of elevation samples every 30 m. The samples are located on the coordinate lines of the Universal Transverse Mercator (UTM) grid, rather than being located along geographic meridians and parallels. These elevation samples are produced by USGS by one of several means (for example, digital scanning of aerial photographs or simply digitizing existing topographic contour maps). On the commonly available 1:24,000 scale topographic quadrangles, the 30 by 30 m grid corresponds to 1.25 × 1.25 mm, and is thus about as fine a grid as can be accurately produced from 1:24,000 scale maps, which are produced to meet a horizontal accuracy criterion of 1/20 of an inch, or 1.27 mm. Note that the models developed by USGS are comprised of elevations samples; thus, the elevations of the model are not intended to represent the mean elevation within a 30 m square.
The USGS also distributes DEMs that are produced by the Defense Mapping Agency (DMA), now known as NIMA. Such DEM's cover 1° by 1° areas and provide elevations every 3" of arc along geodetic parallels and meridians. These elevation models were produced from 1:250,000 scale maps and are not considered accurate enough for the computation of detailed inner-zone terrain corrections; therefore, such DEM's are not used by the terrain correction program to be described, and are not discussed further.
The DEM's that are distributed by the USGS are classified into Levels 1, 2, and 3. Such classifications refer to the overall quality level of the DEM. The DEM's produced thus far by USGS are virtually all Level 1 DEM's. Such models have been produced with the following intentions:
No elevation point will have an absolute error (from mean sea level) of more than 50 m. A desired root-mean-square accuracy of 7 m is considered to be a reasonable desired accuracy for a 30 m DEM; a root-mean-square error of 15 m is the maximum permitted. The relative integrity of elevations within a 7 × 7 array should not exceed 21 m. Elevation statistics are usually produced by USGS for each DEM. Such statistics are based on a minimum of 20 test points, which are locations whose elevation is well known. An elevation is calculated from the DEM at the test point location using bilinear interpolation. The difference in elevation between the interpolated elevation and the test point elevation is used to construct the statistics for the DEM. For the DEM's used thus far, the USGS-supplied root-mean-square errors have never been greater than 7 m.
The accuracy of the DEM's is usually about the same as the topographic quadrangles that they are meant to represent. Therefore, these data can be used to calculate terrain corrections for gravity stations with at least the same accuracy as terrain corrections calculated from topographic maps using more traditional map and template method (e.g, Hammer, 1939), and are potentially more accurate, because elevations from the DEM's are probably more accurate than mean elevations estimated manually from maps. More detailed information regarding the construction of DEM's and the structure of the models themsleves is available from the guide published by USGS (1987). An on-line version of this guide is also available.
The file format for the DEMs distributed by USGS is described in the circular mentioned above. The files are basically ASCII files. The format described by the USGS is one in which each record of the DEM file is a logical record of length 1024 characters. End-of-line characters are presumed to be present at the termination of each logical record. More recently, electronic distrubution of the DEMs has been established. The electronic distribution of the files does not have any end-of-line characters within the DEM file. This format I term a tape image format. The present version of DemRead automatically detects whether a DEM file supplied to it is in the older format or in a tape image format, and processes the file appropriately.
The digital elevation data files supplied by USGS must be preprocessed in order to be used by InnerTC. A program called DemRead is supplied to perform the required preprocessing. The preprocessing is done largely to enhance the execution speed of the program, but the elevation file that is produced by DemRead for input to InnerTC is approximately 30% of the size of the ASCII data provided by USGS, and thus there is also a significant storage savings.
All the ASCII digital elevation files are provided as a simple list to DemRead. This list is provided as a simple text file (the default name is DemRead.inp, but another name can be supplied on the command line). An example of the input file is shown below:
agua_fria.dat frijoles_mesa.dat guaje_mtn.dat
The processed digital elevation data are stored in a binary file. The default name for this file is demdata.bin but the user can optionally provide his or her own name for the output file. This output file is used as the input elevation data file to InnerTC.
After constructing the list of digital elevation files with a text editor, DemRead is executed simply by typing
where inputfilename is the name of the input text file created earlier and outputfilename is the binary elevation file that will be produced by DemRead. Note that one cannot specify the output file name without also specifying the input file name (the order of the command-line arguments is important), but one can specify only an input file name, in which case the elevation data are written to demdata.bin, the default file name.
Required Files. There are three input files that are required for program execution. One of these files, the Gravity Data File, contains all the gravity measurements that are to be terrain-corrected. Another required file is the binary Digital Elevation File. This elevation file must have been created earlier using the program DemRead, which pre-processes the ASCII data files distributed by the U. S. Geological Survey. Finally, a parameter file must be supplied. The parameter file is a simple text file that supplies all the options used by InnerTC.
Optional File. In addition to the three required files, an optional auxiliary elevation file may be used. This file provides elevation data supplemental to the data provided on the DEM. A common situation, especially in a detailed gravity survey, is for a user to be able to use all the elevation data that accompany the actual gravity measurements to supplement the DEM elevations. Or, the user may supplement the DEM data with elevation data from a special survey carried out in conjunction with the gravity survey. Such supplemental elevation data can be used to improve the terrain representation in the vicinity of a gravity station, thus leading to more accurate terrain corrections. When using such auxiliary elevation data, though, the user has to be careful to insure that the auxiliary elevation data are consistent with the DEM elevation data (for example, that no bias exists between the two elevation data sets)
The data file providing the gravity data to be terrain-corrected is a text file that provides a single gravity measurement (along with its associated information, such as location and elevation) per record. Each data record in the file provides the following information, though not necessarily in the order indicated below:
station identifier (a character string of up to 8 characters in length) latitude of the gravity measurement [North latitudes are (+)]. longitude of the gravity measurement [East longitudes are (+)]. elevation of the gravity measurement, in either feet or meters. observed gravity value, in milliGals.
The latitude and longitude can be provided either in (degrees,minutes) or as (degrees,minutes,seconds), or simply as decimal degrees. The elevations can be provided in either feet of meters. Finally, the station identifier can be provided either prior to the horizontal position of the station or after supplying the horizontal position of the station. Specifically, the order in which this information is provided is specified by the value of the fileformat option, which is defined in the parameter file. The interpretation of fileformat and its possible values are shown below.
0 = station,lat,latm, long,longm,elevation,obs gravity
1 = station,latd, longd, elevation,obs gravity
2 = station,lat,latmin,lats,long,longmin,longs,elevation,obs gravity
3 = station,long,longm,lat,latm ,elevation,obs gravity
4 = station,longd,latd, ,elevation,obs gravity
5 = station,long,longmin,longs,lat,latmin,lats,elevation,obs gravity
10 = lat,latm, long,longm,station,elevation,obs gravity
11 = latd, longd,station, elevation,obs gravity
12 = lat,latmin,lats,long,longmin,longs,station,elevation,obs gravity
13 = long,longm,lat,latm ,station,elevation,obs gravity
14 = longd,latd, ,station,elevation,obs gravity
15 = long,longmin,longs,lot,latmin,lats,station,elevation,obs gravity
Note that file formats 10-15 are the same as file formats 0-5
except that the station identifier follows the (x,y) location
information, rather than preceding such data. The nomenclature
used above is interpreted as follows:
lat: degrees latitude (INTEGER), north latitudes (+)
latm: minutes latitude (REAL)
latmin: minutes latitude (INTEGER)
lats: seconds latitude (REAL)
latd: latitude in decimal degress (REAL), north latidues (+)
long: degrees longitude (INTEGER), east longitudes (+)
longm: minutes longitude (REAL)
longmin: minutes longitude (INTEGER)
longs: seconds longitude (REAL)
longd: degrees longitude (REAL), east longitudes (+)
station: station identifier (CHAR), 8 char max, enclosed in 's
elevation: station elevation in either feet or meters (REAL)
obs gravity: observed gravity in milliGal
The elevation file used for input is the file created during the preprocessing by routine DemRead. The file is a binary file that is, essentially, not human-readable. It contains all the elevation data for the DEMs that are to be used for the terrain correction calculation. The name of the elevation file is supplied to InnerTC by specifying it as the value of the variable demfile in the parameterfile.
The auxiliary elevation file is a file containing irregularly-spaced elevation samples. Such data are used to augment the elevation data from the DEM for the innermost terrain correction calculation only. Using auxiliary elevation data can permit one to utilize high-resolution terrain data that may be available in the vicinity of a local gravity survey, for example. Only that portion of the calculation that involves fitting and integrating a terrain surface near the gravity station is affected by the presence of the auxiliary elevation data. Beware of possible biases between the auxiliary elevation data and the terrain data present on the DEM, however (see "Tips on Using the Procedures"). The auxiliary elevation file is associated with the fileformat option, provided in the parameter file.
The format of the auxiliary elevation data file is identical to that used by the gravity data file. In other words, for each auxiliary elevation sample, a record containing a station id, latitude, longitude, and elevation for the sample must be present (an observed gravity value, which is required to be present on the gravity data file, may be present, but is ignored). Elevation units must be the same as used for the gravity data file. The manner in which the station name and coordinates is provided on each record is controlled by the fileformat flag in exactly the same manner as in the case of the gravity data file.
Options used by InnerTC are provided in a simple text file called a parameter file. The only options that are required for execution are the names of the Gravity Data File and the Digital Elevation File. Therefore, the parameter file can be as short as 2 lines long, for default values are taken for all other parameters. However, typically the values of all parameters are specified in the parameter file, except perhaps for the name of the auxiliary elevation file, which may not be available. The values of all parameters are provided using the simple syntax variable=value, with no more than one variable's value provided per line of the parameter file. Note that the case used for the variable's name (that is, upper- lower-, or mixed-case) is irrelevant.
The options that may be specified in the parameter file are as follows.
An example of a parameter file is provided.
InnerTC produces two output files. One is a listing file which provides, in addition to the calculated terrain corrections, a more detailed list of the correction calculation. For example, information about the estimated numerical uncertainty of the integration, and the differences between the interpolated elevation and the user-supplied elevation at a gravity station, are both noted on the listing file. A more detailed listing of the contribution of the nearby versus farther-out terrain is listed here, as well. The second file, called the output file, is a simple ASCII listing of the results of the terrain correction computations. It is intended to be used as input to other terrain correction or gravity reduction programs.
By default, the names of the two output files are InnerTC.lst (the listing file) and InnerTC.out (the output file). However, different names for these files can be specified using the listfile and outfile keywords in the parameter file.
The contents of the listing file are self-explanatory, as headings are provided for all tabular output. No such headings are supplied for the output file, however, and so an example is presented here.
--- Example of output file from InnerTC ---
TESUQ-01 35 48.7510 -105 56.2500 6848.00 125.030 0.267 0.267 0.534
TESUQ-02 35 47.2310 -105 58.2310 6577.00 121.230 0.096 0.130 0.226
TESUQ-03 35 46.4370 -105 54.7390 7061.00 119.320 0.357 0.072 0.429
TESUQ-04 35 47.8630 -105 57.7240 6440.00 122.430 0.050 0.054 0.104
Explanation of output file format:
Field 1: Station ID.
Field 2: Degrees latitude of the station location.
Field 3: Minutes latitude of the station location.
Field 4: Degrees Longitude of the station location.
Field 5: Minutes longitude of the station location.
Field 6: Elevation of the Station (in the user-supplied units).
Field 7: Observated gravity value at the station.
Field 8: Calculated outer portion of terrain correction.
Field 9: Calculated inner portion of terrain correction.
Field 10: Total inner terrain correction (sum of fields 8 and 9), in milliGals.
In versions of InnerTC prior to 5.20, all options were passed to InnerTC on the command line. As the number of options increased, such a method seemed overly complicated. Thus, starting with version 5.20, almost all options are passed in a parameter file. There are now only three options available on the command line.
First, typing InnerTC (that is, executing InnerTC without using any command-line options) will provide a brief summary of syntax to be used to execute the program. The program can be executed in one of two ways:
InnerTC parameterfile
OR
InnerTC -L
The first method supplies the name of the parameter file to InnerTC; the parameter file supplies all the parameters to the program, which reads it and executes. The second method simply displays the software license for the program; no terrain corrections are calculated. Note that only one of the above methods may be used (that is, typing InnerTC -L parameterfile is illegal).
InnerTC is normally distributed in a Win32 version (that is, a version that will execute under Win95/98/NT/ME/2000). In addition, OS/2, DOS, or Linux versions are available (the Linux version is experimental); simply request such a version from Geophysical Software. Different run-time libraries are required depending upon the operating system that is to be used. Installation instructions are available to guide ones installation of the software.
Exit Code Meaning ======== ======= 0 Normal execution (no errors) 1 No arguments provided on command-line, or help requested. 2 Command-line syntax incorrect (parameters not properly paired). 3 Gravity data file not provided or simply does not exist. 4 Elevation (DEM) file not provided or simply does not exist. 7 Neither required input file not provided or does not exist. 10 I/O error while reading Rmin (associated with -r0 flag) 11 I/O error while reading Rmed (associated with -r1 flag) 12 I/O error while reading Rmax (associated with -r2 flag) 13 I/O error while reading fileformat (associated with -ff flag) 14 I/O error while reading elevunits (associated with -eu flag) 15 I/O error while reading elevtype (associated with -et flag) -5 Unrecognized file type for digital elevation file. -4 I/O error occurs while reading the (binary) DEM file -3 Units for DEM coordinates are not feet or meters (must be one or the other) -2 Spatial resolution of the DEM differs in N-S versus E-W directions -1 Unrecognized option provided on the command line. -100-N I/O error on auxiliary elevation file at record N.
The earth reference ellipsoid may be specified, although the default value of 0, which corresponds to the 1866 Clarke Ellipsoid, is appropriate for almost all of the 30-m DEMs produced by the U. S. Geological Survey. However, newer DEMs produced by USGS may use the Geodetic Reference System 1980 (GRS80), which is equivalent to the World Geodetic System 1984 (WGS84). DEMs produced by other countries or for other countries may well use a different reference ellipsoid, in order to be consistent with local mapping practice. Therefore, an option is provided that permits the user to specify which ellipsoid should be used. This information is only used to convert the geographic coordinates of the gravity station locations to the same coordinate system as used by the DEMs. It is essential that consistent coordinate systems be used. For more detailed information on earth reference ellipsoids and related information, the reader is referred to the relevant DoD NIMA publication.
The table below briefy describes the ellipsoids available. For specific ellipsoidal radii of these ellipsoids, see the example parameter file.
Value Description
0 Clarke 1866 Ellipsoid; used by the U. S. for most maps
1 WGS'84. World-wide; the normal ellipsoid used by GPS.
2 WGS'72. World-wide, but has been superseded by WGS'84.
3 Geodetic Reference System 1980 (GRS80). Equivalent to WGS'84.
4 Clarke 1880 Ellipsoid. Used in France and most of Africa.
5 Australian Ellipsoid. Used in Australia (of course!)
6 Krasovsky Ellipsoid. Used by Russia and the former USSR.
7 International. Used in much of the world.
8 Hayford (same as International)
9 Airy Ellipsoid. Used in Great Britain.
10 Modified Airy
11 Bessel 1841 Ellipsoid. Used in Germany and Central Europe.
12 Everest Ellipsoid. Used by India, Pakistan, Burma, Thailand.
13 Modified Everest.
99 User-Specified. For cases not covered by the ellipsoids listed above.