Gravity Reduction and Terrain Correction Program

Table of Contents

Overview

A program designed to reduce measured gravity data to free-air and Bouguer gravity anomalies is described. The program also calculates terrain corrections appropriate for the Bouguer gravity anomalies. Such corrections are calculated from digital terrain data provided by the user; the digital terrain data are used to correct for the effect of terrain from an inner radius provided on the digital terrain file to an outer radius also provided on the digital terrain file. The outer radius is typically 166.7 km; the inner radius used is provided by the user, but can be no smaller than one-half the length of the diagonal of the smallest terrain compartment.

The program is a modification of one developed by Donald Plouff of the U. S. Geological Survey (Plouff, 1966). The user must supply two input files for execution. One of these files contains the measured gravity data, including station identifiers, geographic coordinates, elevations, measured gravity values, and so-called "inner-zone" terrain corrections (the latter may be 0). The other file contains the digital terrain data, which must be provided in a prescribed format. Essentially, the digital terrain data are used to correct the gravity data for the effects of terrain from the maximum distance used in the inner-zone correction to a large distance from each station, typically 166.7 km.

The current version was adapted in May 1981 from a version supplied by Don Plouff of the U. S. Geological Survey (Plouff, 1966). Modifications made at that time were largely in the area of i/o and use of a more accurate method of calculating the gravity reference field. Terrain corrections are calculated exactly as the Plouff version would calculate such corrections.

Input Files

Gravity Data Input File

The gravity data are supplied to the program on a sequential file that essentially consists of blocks of card images as follows.
  • Image 1 Title Card Image 1
  • Image 2 Title Card Image 2
  • Image 3 . . . . .
  • Image 4 . . . . .
  • Image 5 . . . . .
  • Image N Title Card Image N
  • Image N+1 last
  • Image N+2 Running Title Card Image (only one line permitted)
  • Image N+3 Parameter Card for Gravity data (described below)
  • Image N+4 Card Image containing format for data input
  • Image N+5 Gravity Data Card Image 1
  • Image N+6 Gravity Data Card Image 2
  • Image N+7 . . . . .
  • Image N+8 . . . . .
  • Image N+9 . . . . .
  • Image N+4+M Gravity Data Card Image M
    •  Any number of title card images can be placed in the first block. This block of card images is simply echoed on the listing file created by the program; its purpose is to permit the user to enter any comments or other pertinent information regarding the gravity data being processed onto the listing file, as well as into the file that is being processed. The block of title card images is terminated by entering the 4-character string "Last" into the card image immediately following the title block (the case, upper, lower, or mixed, of the 4-character string is irrelevant). Note that blank lines can be inserted into the title card block. Immediately following the title card block is a single card image that contains a running title. The running title is placed on each page of the listing file that is to be created. Only the first 80 characters of the running title are used.

    Following the running title, a card image containing parameters to be used during the reduction process is required. The variables read from this card image and their meanings are as follows.

    1. NTYPE The number of map sets that will be provided on the map file.
    2. METRES An INTEGER switch. If METRES=0, the elevations accompanying the gravity data are assumed to be in feet; else, they are assumed to be in metres.
    3. GCONST DOUBLE PRECISON scalar, a constant that is to be added to all the gravity data read from the input file. The use of an additive constant permits the data to be recorded as values with respect to a base station. The base station value can then be added as GCONST.
    4. DENSE2 REAL scalar, a 2nd reduction density for the Bouguer gravity anomalies. The first reduction density is always 2.67 grams/cc. The default value for DENSE2 is 2.45 grams/cc.
    5. SOURCE A 4-character code used to indicate the organization that supplied the gravity data to be reduced. For example, "USGS" or "LANL" are codes that might be used to indicate the origin of the data.
    6. SGNLAT A 1-character code reflecting the sign of the latitudes of the gravity stations. If SGNLAT is blank, "+" is assumed. "N" or "S" are common alternatives.
    7. SGNLON A 1-character code reflecting the sign of the longitudes of the gravity stations. IF SGNLON is blank, "+" is assumed. "E" or "W" are common alternatives.
    8. APPLYCC An INTEGER switch; if APPLYCC=0, the normal curvature correction is applied when calculating the Bouguer gravity anomaly. Otherwise, no curvature correction is applied. Note that the value of APPLYCC only affects the normal calculation of the Bouguer gravity anomaly; it does not affect whether a curvature correction is applied during the calculation of the terrain correction (the parameter RGR, which is provided on the first record of the map data file, can be used to control the application of a curvature correction during the terrain correction calculation).
    9. DEBUG An INTEGER switch; if DEBUG=0, only the total machine correction for each gravity station is listed. Otherwise, the terrain correction for each digital terrain quadrangle and for each station is listed. Setting DEBUG non-zero may produce voluminous printed output, and is thus normally used when the user suspects that incorrect results are being produced.

    Example Input File

    The list below shows an example of the gravity data input file: 
    The gravity data below are provided as an example of an input file to
the terrain correction program of Don Plouff, as modified by Allen Cogbill.
Any number of title card images may be used as part of a "header" of
information. Such information is simply printed (echoed) to the output
listing file that the program creates; the facility is designed to permit
comments regarding the data set to be printed along with the terrain
correction and reduced gravity data set.
last
Example of a running title card image; this image limited to <= 80 chars.
3 1 979000.00d0  2.20  'TEST' '+' '-'  0 0
(a8,i2,f8.4,i5,f8.4,f9.2,4x,f7.3,15x,f6.3,1x,a20)
800     31 10.5420 -105 21.7799  1403.04 979004.011  0.136  0.015  0.151 Comment # 1
200     31 28.4880 -105 13.2600  1325.40 979064.778  0.023  0.004  0.027 Comment # 2
201     31 28.2360 -105 13.2600  1330.22 979062.765  0.018  0.006  0.025 Comment # 3
202     31 28.0380 -105 13.3200  1333.55 979061.327  0.018  0.006  0.024 Comment # 4
    The map parameter card image is read using a Fortran free-format read, with variables either comma- or blank-delimited. Thus, CHARACTER variables (such as SOURCE) should be enclosed in single quotation marks.

    Immediately following the map parameter card image is a card image containing the user-supplied data input format for the gravity data to follow. This card image must contain a format that conforms to Fortran syntax; the format must be begin and terminate using a left and right parenthesis. It should be designed to permit a single gravity observation and associated data to be read by the reduction program for each execution of a Fortran READ. The gravity data are placed in the file immediately following this card image. Each record of the gravity data must contain the following information, in the order shown below.

    1. IDENT Station identifier for the gravity station. IDENT is limited to 8 characters in length.
    2. LAT INTEGER scalar, degrees latitude of the gravity station.
    3. ALATM REAL scalar, minutes latitude of the gravity station.
    4. LONG INTEGER scalar, degrees longitude of the gravity station.
    5. ALONGM REAL scalar, the minutes longitude of the gravity station.
    6. ELEV REAL scalar, the elevation with respect to mean sea level of the gravity station. Note that the units of ELEV are either metres or feet, depending upon the value supplied for the switch METRES on the parameter card.
    7. OBS REAL scalar, the observed gravity value, in milliGals, at the gravity station. Note that the constant GCONST (as supplied on the parameter card) will be added to the value of OBS prior to performing any reduction calculations.
    8. TCHAND REAL scalar, the value of the "inner-zone" terrain correction for the gravity station. TCHAND must be supplied in milligals and be appropriate for a reduction density of 2.67 grams/cc.
    9. REMARK CHARACTER string of length 40 characters or less. The string in REMARK is not used by the program, but is placed on the output files created by the program. Blank fields are acceptable.
    Any number of gravity stations can be placed in the input file, although on 500 stations are processed at a time. If more than 500 stations are provided, stations are internally processed in blocks of 500 stations (or less, for the last block). All the values read by the reduction program are echoed on the listing file prior to any processing being performed on the data.

    Return to Table of Contents

    Map Data Input File

    If terrain corrections are desired, a map data input file must be provided. The map file is an ASCII file that not only provides the digital terrain data needed to calculate the terrain corrections, but also provides parameters that dictate how the terrain correction calculation will be performed. The map data are provided as one or more map sets, each set corresponding to a particular compartment size. The first map set is the set of maps having the smallest compartments, the next set has larger compartments, etc. The compartment size of each set of maps must be an integral multiple of the compartment size of the previous map set. For example, suppose the first map set has a compartment size of 15 sec of arc. The next map set might then contain compartment sizes of 1 minute (i.e, 60 sec) of arc. The 3rd map set might be made up of compartments of, say, 2 minutes or 4 minutes of arc, etc.

    The provision of providing multiple map sets permits the program to calculate terrain corrections very rapidly, using small compartment sizes close to a gravity station, and much larger compatment sizes at greater distances. Construction of the various map data sets required to produce such efficiency can be tedious. Alternatives are to either use a single map set having a small compartment size or to automatically produce the required map sets using a computer program. The later is not generally available, as far as I know, but special instances exist for the use of 30-sec compartment sizes.

    There is no requirement that multiple map sets be provided. The entire terrain correction can be readily calculated using a single type of terrain data. The ability to combine different map types by using overlapping map sets is strictly a provision that can supply very high efficiency in calculating the terrain corrections. From the user's point of view, using a single type of map data set is much easier to set up, and may be preferable in some cases.

    Description of the Map Data File

    The map file contains almost all the information needed to process the elevations used to calculate the terrain corrections (an exception is that the number of map data sets is supplied from the gravity data input file). Map data sets consist of quadrangle-formatted elevation data; each map "set" is defined by having a quadrangles of the same size as well as having terrain compartments that all all the same size within a quadrangle. The general flow of the map file is illustrated below: 
      Header for Map Data Set
      Quadrangle Header Record
           Elevation Record 1
           Elevation Record 2
                .       .   .
                .       .   .
                .       .   .
           Last Elevation Record
      Quadrangle Header Record
           Elevation Record 1
           Elevation Record 2
                .       .   .
                .       .   .
                .       .   .
           Last Elevation Record
      Quadrangle Header Record
           Elevation Record 1
           Elevation Record 2
                .       .   .
                .       .   .
                .       .   .
           Last Elevation Record
      Quad Header with '999999999' in first 12 columns
      (the nines signal the end of this map data set)
    If more than one map data set is to be used, the additional sets are simply appended, using the same organizational structure, to the previous map data.

    Specific format of the map data file

    Record 1; Fortran format: (2i3,2i2,2i3,4f8.4,12x,2i2). This record is the header record for all the map quadrangles of a map set. 
    LATM    N-S length of the quadrangles, in minutes of arc.
LONM    E-W length of the quadrangles, in minutes of arc.
KCLAT   N-S length of the smallest compartment in the quad, in min or sec.      
LCLON   E-W length of the smallest compartment in the quad, in min or sec.
KBLOK   N-S length of composite block, in min or sec
LBLOK   E-W length of composite block, in min or sec
RMAX    Maximum distance to which terrain correction will be carried (km)
RMIN    Distance at which terrain correction will be started (km)
RMED    Intermediate distance, at which the conversion to composite blocks occurs (km)
RGR     Distance beyound which curvaure correction is incorporated
        into terrain correction calculation (km)
KCIRC   Set to 0 if piecemeal join to previous terrain correction block is
        desired; else, set to non-zero if a circular inner join to
        previous terrain corrections is desired. Typically, KCIRC is set
        to non-zero for the first (that is, the closest) map set, then set
        to 0 for succeeding sets. The circular joint permits the
        calculated correction to be correctly added to "inner-zone"
        corrections, which are normally carried out to a particular
        radius.
KSEC    Set KSEC to 1 if compartment sizes and block sizes are provided in
        minutes; else, set KSEC to 60 if compartment sizes and block sizes are
        provided in seconds. Note that quadrangle sizes (LATM, LONM)
        are always provided in minutes.
    Record 2; Fortran format: (a12,2i3,i4,i3,i2). This record is the header record for a particular map quadrangle. After reading this record and the associated elevation data for the quadrangle, execution continues by reading another quadrangle, and so forth, until a quadrangle name of all nines (i.e., '999999999999') is read. Note that the elevaton data must be terminated by a record having nines for its quad name; an end-of-file condition will terminate the program prematurely. 
    DENTQ   -  12-character name of the quadrangle.
            Note that the end of this particular map set is signaled
            by placing all 9's in columns 1 through 12.
LTDQ    -  Latitude degrees of NW corner of the quad.
LTMQ    -  Latitude minutes of NW corner of the quad.
LNDQ    -  Longitude degrees of NW corner of the quad.
LNMQ    -  Longitude minutes of NW corner of the quad.
KFATH   -  Controls the interpretation of the elevations on the quadrangle.  Specifically,
       A)  If KFATH < 31,
           A.1  All elevations are in feet unless KFATH = 6.
                If KFATH <> 6, negative elevations imply
                sea-water compartments.
           A.2  If KFATH = 6, positive elevations are in feet,
                negative elevations are in fathoms and imply
                sea-water compartments.
       B)  If KFATH = 31, positive elevations are in metres,
           negative elevations are in feet and imply sea-water
           compartments.
       C)  If KFATH > 31,
           C.1  If KFATH = 32, same as KFATH = 31 [see (b)].
           C.2  If KFATH = 33, both positive and negative
                elevations are in metres and negative elevations
                imply sea-water compartments.
           C.3  If KFATH > 33,
                C.3.1  If KFATH <> 73, positive elevations
                       are in metres and negative elevations
                       are in fathoms and imply sea-water comp.
                C.3.2  If KFATH = 73, both positive and
                       negative elevations are in metres and
                       negative elevations imply land compartments.
    Elevation data; Fortran format (20f5.0). Following the quadrangle header record, the elevation data for the quadrangle appear on subsequent records. The elevation data are arranged in scan-line format, starting with the compartment elevations in the NW corner of the quadrangle. Scan-line format implies that the data are arranged in west-east rows, as illustrated below for a 15-minute quadrangle having compartments of size 30 seconds of arc: 
            NW Corner
           |
           |
  Row 1    X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
           X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
  Row 3    X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
           X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
  Row 5    X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
           X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
  Row 7    X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
           X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
  Row 9    X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
           X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
  Row 11   X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
           X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
  Row 13   X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
           X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
  Row 15   X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
           X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
  Row 17   X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
           X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
  Row 19   X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
           X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
  Row 21   X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
           X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
  Row 23   X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
           X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
  Row 25   X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
           X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
  Row 27   X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
           X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
  Row 29   X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
  Row 30   X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
                                                                     |
                                                                     |
                                                                 SE Corner
     

    Following the elevation data for this quadrangle, control returns to "Record 2"; that is, another quadrangle header record is read. When a quadrangle name of all nines is read from the quad header, the program attempts to read "Record 1"; that is, the header record for the next map data set, if another data set has been requested.

    Example File. An example (145 Kb) of a map data file is available that may help clarify the organization of the file.

    Remark. A sometimes-annoying feature of the terrain correction program is that the number of map data sets to be used to terrain-correct the gravity data is supplied in the gravity data input file (it is the first variable read from that file), rather than in the map data input file. This can lead to easy mistakes by the user, when the number of map ata sets is entered incorrectly in the gravity data input file.

    Return to Table of Contents

    Command-line Options

    The only command-line options are the names of the input or output files that are to be used. The syntax for runnng outertc is

    outertc in=inputfile list=listfile map=mapfile out=mapfile

    where inputfile is the input file [contains the gravity data], listfile is the name of the listing file that will be created, and mapfile is the name of the map data file used for outer-zone TCs. An additional file called gravcorr.dat is created by the program. This last file contains the results of the gravity reduction and the terrain correction calculations (if requested). If no map data file is placed on the command-line, no terrain corrections will be calculated; the gravity data will have standard reductions applied, minus any terrain corrections.

    The only file required for input is the gravity data input file. If no listing file name is supplied, the default name of gravlist.lst will be used.

    Return to Table of Contents

    References Cited

    Plouff, Donald, 1996, Digital terrain corrections based on geographic coordinates, paper presented at the 36th annual meeting of the Society of Exploration Geophysicists, Houston, Texas, 8 pages + 6 figures.

    Return to Table of Contents