Allen H. Cogbill 398 Venado Street Los Alamos, NM 87544-2435 505-216-1611 505-412-9328 [cell] 720-367-7833 [FAX] E-Mail: acogbill@geopotential.com SkypeID: acogbill
Copyright 1998-2010, Geophysical Software, all rights reserved. The program RasterTC© is licensed, not sold, and may not be redistributed without express written consent of Geophysical Software.
DISCLAIMER
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.
RasterTC is a program for calculating terrain corrections for land gravity stations. It is similar both conceptually and procedurally to its sister program InnerTC. A user provides gravity station locations as well a grid of terrain data; from this information, extremely accurate terrain corrections are calculated by. Currently, terrain grids may be provided in either Geosoft Version 1 format or in the so-called GS binary format of Golden Software; the latter format is the standard grid format used by Golden Software's SURFER program (version 6.04 or earlier).
Reports of bugs in the program and suggestions for improvements should be sent to the author. In particular, users needing support for other formats should not hesitate to send their needs. The preferred method for reporting bugs or providing suggestions is by electronic mail.
RasterTC is designed to calculate extremely accurate terrain corrections out to several km from a gravity station. It uses a series of terrain grids that are provided in rectanguar coordinates. Each successive grid is coarser in discretization than the grid used before it; the use of a series of grids is done strictly for computation speed, and is not absolutely necessary: one could use a single grid, for example. Also, please note that RasterTC does not provide gravity reductions, but only calculates terrain corrections.
The original terrain correction methodology was developed 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 DTM 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 use of 30-m, USGS DEMs was convenient in the United States (apart from Alaska and Hawaii), but the USGS DEM format is unique to USGS; thus, elsewhere in the world, the method could not be used. RasterTC is an attempt to make available the high accuracy and overall methodology used by InnerTC to other localities: users who can provide their own terrain grids can use RasterTC to calculate high-accuracy terrain corrections for their gravity stations.
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. Other required files are the binary digital terrain grid(s) (hereafter referred to as a DTM). Up to five terrain grid files may be used. Such grids should be provided in decreasing horizontal resolution (that is, grids having finer discretization should be provided first. Finally, a parameter file must be supplied. The parameter file is a simple text file that supplies all the options used by RasterTC.
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 DTM. 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 DTM elevations. Or, the user may supplement the DTM 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 DTM 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 16 characters in length) Easting (x-coordinate) of the gravity measurement Northing (y-coordinate) of the gravity measurement Elevation of the gravity measurement, in either feet or meters. Observed gravity value, in milliGals.
The horizontal coordinates (Eastings and Northings) should be provided in meters.. The elevations can be provided in either feet or 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 = stationID,easting,northing,elevation,obs gravity
1 = stationID,northing,easting,elevation,obs gravity
2 = easting,northing,stationID,elevation,obs gravity
3 = northing,easting,stationID,elevation,obs gravity
4 = easting,northing,elevation,obs gravity
5 = northing,easting,elevation,obs gravity
6 = easting,northing,elevation
7 = northing,easting,elevation
easting: east (x-) coordinate of the gravity station, in meters
northing: noth (y-) coordinate of the gravity station, in meters
stationID: station identifier, 8-characters max.
elevation: station elevation in either feet or meters.
obs gravity: observed gravity (mGal)
The auxiliary elevation file is a file containing irregularly-spaced elevation samples. Such data are used to augment the elevation data from the DTM 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 DTM, 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, easting, northing, and elevation for the sample must be present (although a station ID and observed gravity value may be provided, their values are ignored). The elevation units must be the same as used by 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 RasterTC 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 Terrain 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.
debug A Boolean switch. If true, cosiderable output is printed on stdout (the console); this output can be redirected to a file. Also, if true, the elevations associated with each inner-zone gravity calculation are written to files named NODExxxx.DAT, where xxxx = the sequence number of the gravity station. Note that considerable output can be sent to stdout when debug is on; therefore, it is usually advisable to redirect stdout to a file when debug is on.
verbose A Boolean switch. If true, more information than usual is printed to the listing file.
gravfile The name of the Gravity Data File (REQUIRED, no default).
auxelevfile The name of the Auxiliary Elevation File (OPTIONAL, no default)
listfile The name of the Listing File to be created (default: RasterTC.lst)
outfile The name of the Output File to be created (default: RasterTC.out)
csvfile The name of the comma-separated-value (csv) File to be created (default: RasterTC.csv)
rmin The minimum radius for the terrain correction calculation (default: 30 meters)
rmed The intermediate radius for the terrain correction calculation (default: 250 m)
angle The angular increment used for the azimuthal portion of the innermost terrain correction calculation (that is, for the numerical integration). This value must be provided in degrees and must divide evenly in 360. The default is 6 degrees.
rcurve The distance, in meters, beyond which the effect of the earth's curvature is taken into account. The default value is 15,000 meters. If rcurve is set to 0, the earth's curvature will not be taken into account.
gridfile The name of one of the DTM files to be used.
rmaxgrid The maximum radius for the terrain correction calculation (no default). Note that the parameters gridfile and rmaxgrid are provided in pairs: up to five pairs of such parameters may be provided.
gridtype A switch indicating which format is used for the DTM. Currently, gridtype should be either 1 (for Golden Software's GS binary format) or 2 (for Geosoft version 2 format). The default type is GS binary.
density The terrain density to be used, in grams/cc (default: 2.67)
elevunits A switch indicating whether the elevation units supplied with the gravity observations are feet (0) or meters (1). The default is 0 [feet].
elevbias A constant that is added to all the DTM elevations prior to using such elevations for the calculation of the terrain corrections. Typically, elevbias is used when the user provides auxiliary elevations. Such auxiliary elevations frequently differ by a constant (or nearly a constant) from the DTM elevations; using a non-zero elevbias value can force the two elevation sets to be consistent with one another.
fileformat A switch indicating the type file format that used for the Gravity Data File. Legal values are 0,..,7; the default value is 0. The uses of the fileformat flag are explained in the section describing the Gravity Data File.
An example of a parameter file is provided.
RasterTC produces three 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 text 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. The third file is a comma-separated-variable (csv) file, which is very useful for importing the results of the calculations into a spreadsheet application such Microsoft Excel or OpenOffice's scalc.
By default, the names of the three output files are RasterTC.lst (the listing file), RasterTC.out (the output file), and RasterTC.csv (the csv file). However, different names for these files can be specified using the listfile, outfile, and csvfile keywords in the parameter file. Examples of both the listing file and the output file are provided.
Almost all of the options needed for RasterTC's processing are provided in the parameter file. There are now only three options available on the command line.
First, typing RasterTC (that is, executing RasterTC 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:
RasterTC parameterfile
OR
RasterTC -L
The first method supplies the name of the parameter file to RasterTC; 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 RasterTC -L parameterfile is illegal).
RasterTC is normally distributed in a Win32 version (that is, a version that will execute under /NT/2000/XP/Vista/Win7).
Indicators of quality are provided for both the innermost terrain correction and the outer terrain correction. On the listing and the output files, such indicators are listed under the column headings "QF-Inner" and "QF-Outer", respectively. If the TC calculation is OK, the quality factors should be 0 for both the inner- and the outer-zone portion of the calculation. Otherwise, the quality factors provide an indication of what might be wrong.
A rough indicator of how well the terrain in the immediate vicinity of a gravity station is represented by the available elevation samples is obtained by examining the spatial distribution of the elevation samples. In the radial interval Rmin to Rmed, RasterTC counts the number of samples falling within the 8 octants surrounding the station. If any of these octants are missing elevation samples, that fact is noted, and the tabulated quality factor simply notes how many of octants are missing samples.
|
QF-Inner |
Meaning of Error Code |
|---|---|
|
0 |
Inner-zone terrain calculation OK |
|
1 |
No elevation samples occur in 1 octant surrounding the gravity station |
|
2 |
No elevation samples occur in 2 octants surrounding the gravity station |
|
3 |
No elevation samples occur in 3 octants surrounding the gravity station |
|
4 |
No elevation samples occur in 4 octants surrounding the gravity station |
|
5 |
No elevation samples occur in 5 octants surrounding the gravity station |
|
6 |
No elevation samples occur in 6 octants surrounding the gravity station |
|
7 |
No elevation samples occur in 7 octants surrounding the gravity station |
|
22 |
Duplicate elevation nodes encountered while calculating terrain gradients |
|
23 |
All elevation nodes collinear or triangulation structure corrupted |
|
24 |
Invalid parameters passed to gradient calculation routine |
|
26 |
Convergence not attanined in calculation of terrain gradients |
|
48 |
Internal logic error while attempting to delete duplicate nodes |
|
49 |
Unable to allocate memory while deleting duplicate nodes |
|
96 |
Internal error occurs while attempting to triangulate nodes near station |
|
99 |
Less than 3 elevation samples are available near station |
|
101 or greater |
Duplicate nodes occur in elevation samples near station |
|
-99 |
Unable to allocate sufficient memory for inner-zone TC calculation |
The QF-Outer codes are simple to interpret. A result of 0 means that the outer-zone calculation proceeded successfully. If a portion of the outer-zone terrain is missing from the elevation grids supplied, the value of QF-Outer will reflect the per cent of terrain that was available (rounded to the nearest per cent). For example, if QF-Outer is 91, the implication is that 9% of the terrain in the outer zones was missing for some reason, and that the terrain correction calculated for that particular station is too small by some amount.