r.gravity.terrain
A GRASS tool to calculate gravity terrain corrections
r.gravity.terrain elevation=name points=name output=name minimum_distance=float maximum_distance=float crustal_density=float [nprocs=integer] [--overwrite] [--verbose] [--quiet] [--qq] [--ui]
Example:
r.gravity.terrain elevation=name points=name output=name minimum_distance=0.0 maximum_distance=0.0 crustal_density=2670
grass.script.run_command("r.gravity.terrain", elevation, points, output, minimum_distance, maximum_distance, crustal_density=2670, nprocs=0, overwrite=None, verbose=None, quiet=None, superquiet=None)
Example:
gs.run_command("r.gravity.terrain", elevation="name", points="name", output="name", minimum_distance=0.0, maximum_distance=0.0, crustal_density=2670)
grass.tools.Tools.r_gravity_terrain(elevation, points, output, minimum_distance, maximum_distance, crustal_density=2670, nprocs=0, overwrite=None, verbose=None, quiet=None, superquiet=None)
Example:
tools = Tools()
tools.r_gravity_terrain(elevation="name", points="name", output="name", minimum_distance=0.0, maximum_distance=0.0, crustal_density=2670)
This grass.tools API is experimental in version 8.5 and expected to be stable in version 8.6.
Parameters
elevation=name [required]
Name of input elevation raster map
points=name [required]
Name of input vector map
Input vector map of points containing the gravity stations
output=name [required]
Name for output file
minimum_distance=float [required]
Minimum distance from gravity station to calculate terrain correction
maximum_distance=float [required]
Maximum distance from gravity station to calculate terrain correction
crustal_density=float [required]
Crustal density in kg/m^3 used in the terrain correction (e.g.2670 kg/m^3)
Default: 2670
nprocs=integer
Number of threads for parallel computing
0: use OpenMP default; >0: use nprocs; <0: use MAX-nprocs
Default: 0
--overwrite
Allow output files to overwrite existing files
--help
Print usage summary
--verbose
Verbose module output
--quiet
Quiet module output
--qq
Very quiet module output
--ui
Force launching GUI dialog
elevation : str, required
Name of input elevation raster map
Used as: input, raster, name
points : str, required
Name of input vector map
Input vector map of points containing the gravity stations
Used as: input, vector, name
output : str, required
Name for output file
Used as: output, file, name
minimum_distance : float, required
Minimum distance from gravity station to calculate terrain correction
maximum_distance : float, required
Maximum distance from gravity station to calculate terrain correction
crustal_density : float, required
Crustal density in kg/m^3 used in the terrain correction (e.g.2670 kg/m^3)
Default: 2670
nprocs : int, optional
Number of threads for parallel computing
0: use OpenMP default; >0: use nprocs; <0: use MAX-nprocs
Default: 0
overwrite: bool, optional
Allow output files to overwrite existing files
Default: None
verbose: bool, optional
Verbose module output
Default: None
quiet: bool, optional
Quiet module output
Default: None
superquiet: bool, optional
Very quiet module output
Default: None
elevation : str, required
Name of input elevation raster map
Used as: input, raster, name
points : str, required
Name of input vector map
Input vector map of points containing the gravity stations
Used as: input, vector, name
output : str, required
Name for output file
Used as: output, file, name
minimum_distance : float, required
Minimum distance from gravity station to calculate terrain correction
maximum_distance : float, required
Maximum distance from gravity station to calculate terrain correction
crustal_density : float, required
Crustal density in kg/m^3 used in the terrain correction (e.g.2670 kg/m^3)
Default: 2670
nprocs : int, optional
Number of threads for parallel computing
0: use OpenMP default; >0: use nprocs; <0: use MAX-nprocs
Default: 0
overwrite: bool, optional
Allow output files to overwrite existing files
Default: None
verbose: bool, optional
Verbose module output
Default: None
quiet: bool, optional
Quiet module output
Default: None
superquiet: bool, optional
Very quiet module output
Default: None
DESCRIPTION
The purpose of r.gravity.terrain is to compute gravity terrain corrections used in Complete Bouguer anomaly calculations from digital elevation models (DEM). Parameter elevation is the DEM used in the terrain correction. The region resolution should be set to match the DEM. Calculations are done in meters, so an appropriate coordinate system (such as UTM) should be used. Parameter points is the vector points file with the measured gravity station locations. Parameter output is the output file with the calculated terrain corrections in mgal. This tool also allows parallel processing to speed calculations.
NOTES
This tool calculates the final term in the below equation (g(Terrain correction)).
g(Bouguer complete) = g(Observed) - g(Latitude) + g(Free Air) - g(Simple Bouguer) + g(Terrain correction)
A simple Bouguer correction approximates the gravity of rock from sea level to the measurement location as a flat slab. This works well in areas of low topographic relief. The gravity terrain correction here accounts for the 3-D variability of crustal gravity due to spatial variations in elevation. In mountainous regions, this can be larger than the geological signature the gravity measurements are used to interpret.
This GRASS tool uses the Nowell (1999) formula 7 (after Blais and Ferland, 1984 Eq 4) to calculate terrain corrections for each DEM point with respect to each gravity measurement location. The tool must be used in a projected coordinate system as the calculations depend on Cartesian geometry.
Formula used in terrain correction
g = G * rho * Acs * ((1/d) – (1/d’))
This (g) is the gravitational effect difference from a uniform slab correction for a single DEM pixel with respect to the gravity station measurement location. It is a calculation for a flat-topped rectangular prism. There are other methods that could be implemented, but this is computationally efficient. Other methods include using prisms with slanted tops to better "fit" topography or fitting polynomial curves to topographic data. The thought is that with higher resolution DEMs there will be less difference between the various methods, and the one implemented here is computationally fast so that higher resolution DEMs can be used more readily. All calculations are done in meters, kg, and seconds, but final terrain correction values are reported in mgal (m/s^2 *100,000).
Formula variables
G = 0.000000000066743:
This is the universal gravitational constant ((m^3/(kg*s^2)) and is fixed.
Rho = 2670 kg/m^2:
This is the rock density at a DEM pixel. Typically, in Bouguer terrain corrections
this is a fixed value across all points. This GRASS tool lets the user determine
the density used in the terrain correction. Future versions could spatially vary
this according to a raster map for more realism. This is an entered variable.
Acs = 8100 m^2: (for a 90 meter grid spacing)
This is the cross-sectional area of a square DEM pixel in meters. It matches
the GRASS region grid spacing and is automatically computed.
d = ((X2 - X1)^2 + (Y2 - Y1)^(1/2)):
Distance from Gravity station to the center of the DEM pixel at constant elevation
(calculated using the distance formula)
X1 = gravity station easting,
Y1 = gravity station northing,
X2 = DEM pixel center point easting,
Y2 = DEM pixel center point northing
d’ = (((d^2) + (hdif^2))^(1/2))):
Distance from the Gravity station to the top of the DEM prism point in meters
(calculated using the Pythagorean theorem)
hdif = | Z2 - Z1 | (absolute value)
Z1 = Gravity station elevation in meters
Z2 = DEM point elevation in meters
Use of minimum and maximum distance parameters
These parameters set a minimum and maximum distance (d) from the Gravity station measurement (X1, Y1) for DEM terrain correction values to be calculated and summed.
The terrain correction formula is an approximation and gives incorrect results if the closest DEM center is too close to the Gravity measurement point. The threshold value is empirically between 1/3 and 1/2 the DEM grid spacing value. For example, when using a 90m DEM, a minimum value of 45-90 m is suggested.
The maximum value is also subjective. Calculating a terrain correction for all DEM points out to 167 km from the station measurement point is considered complete as that is the horizon that affects an individual point on a spherical Earth. However, historically, full calculations out to 167 km at a single grid spacing are not done due to computational intensity. This may or may not be an issue anymore. To address this, an inner terrain correction at a smaller grid spacing and an outer terrain correction at larger grid spacing have been used. The inner and outer corrections are summed to create a complete gravity terrain correction. It might still be useful to do this if high-resolution DEM data (e.g., 3 meter LIDAR) is used for the inner terrain correction. This would involve running the terrain correction routine twice with different grid (region) spacing and minimum/maximum distances. The outer and inner terrain corrections are then summed together.
EXAMPLE
g.region raster=90mSRTM_DEM
r.gravity.terrain input=90mSRTM_DEM maximum_distance=30000 minimum_distance=90 points=stations output=corrections.txt
Figure 1: This is an example of a gravity terrain correction calculated out to
30km from the measurement site (located in the Sangre de Cristo Mountains /
Spanish Peaks of Colorado). This tool does not output the above raster, but
instead it is a visualization of an intermediate calculation product. The terrain
correction value at each DEM pixel is summed and then output in this tool. The
color scale shows gravity terrain correction values of individual pixels in mgal.
For the above gravity point, the total terrain correction is 8.56 mgal for
an n of 3,0163,432 cells.
Figure 2: Above is a comparison of the terrain correction output of this tool with
values from the El Paso PACES gravity database downloaded in 2015 using 90m SRTM
DEM data out to 167 km. The dataset is 131 gravity points from the Sangre de Cristo
Mountains and Raton basin of Colorado. There is a strong linear correlation between
the two correction techniques (R^2 of 0.989), but the one used here outputs somewhat
smaller values at high terrain correction magnitudes. However, this could also be
due to the small number of points with terrain correction values greater than
20 mgal.
REFERENCES
Nowell, D.A.G., 1999. Gravity terrain corrections—an overview. Journal of Applied Geophysics, 42(2), pp.117-134.
Blais, J.A.R. and Ferland, R., 1984. Optimization in gravimetric terrain corrections. Canadian Journal of Earth Sciences, 21(5), pp.505-515.
SEE ALSO
AUTHORS
David W. Farris - East Carolina University
e-mail: farrisd19@ecu.edu
https://geology.ecu.edu/people/davidfarris/
SOURCE CODE
Available at: r.gravity.terrain source code
(history)
Latest change: Thursday Jun 12 13:09:39 2025 in commit 671c05d