Introduction
Forty years ago, we spent a few years developing software for directional analysis of gravity and magnetic field data. The project grew out of an interest in determining how strongly linear structures observed in landscape features were reflected in the geophysical data. Gravity data, and to some extent airborne magnetic data, often have deeper sources than most surface features and can therefore provide a different perspective on tectonic structure.
The software was written entirely in Fortran. At the time, personal computers were quite limited for this kind of work, so the main version ran on a Unix mainframe, while a scaled-down version was compiled for the PC. Although this arrangement served its purpose, the software was not easily accessible to others.
Recently, we decided to revisit the old files, translate them into Python programs with the help of AI tools, and make them available to anyone with an interest in this work. Although directional analysis was our original objective, a number of additional tools had to be developed along the way and they are included here.
The programs are described in the sections listed in the sidebar. They can also be downloaded as ZIP files, along with some data files. The sidebar also includes a section containing published papers and public talks based on this work.
Freyr Þórarinsson and Stefan G. Magnusson
Directional filtering
Maps of topography, gravity fields, and other such phenomena are characterized by various components: a regional slope in one direction, hills and valleys, and linear structures that reflect fracture systems in the Earth’s crust. Interest in studying the last of these was the motivation for the program STRIKE, and others that accompanied it.
Frequency analysis transforms amplitude into frequency and phase. Two-dimensional frequency analysis of a map converts it into two components, frequency and phase. The frequency plot has the lowest frequency at the center and the highest at the edges. The direction of a ray from the center combines all features with the same orientation in the original map into a single ray. Instead of [x, y] coordinates, one obtains [direction, frequency]. This makes it possible to examine the strength of amplitudes in different directions, isolate strong directions with directional filters, and transform them back into [x, y] coordinates in order to locate them on the original map.
The program Strike.py includes routines for reading a data file, gridding the data, and selecting outlines in the landscape to examine together with the data. It also contains routines for preparing the transformation into the frequency domain by removing the background trend and tapering the data down to zero at the edges. The FFT routines handle the forward and inverse transformation of the data between [x, y] and [direction, frequency]. In the frequency domain, it is possible to select directions and frequencies of interest and filter out the rest. In addition, there are routines for examining the frequency domain in various ways.
A more detailed account of these routines is given in the manual for Strike. It may also be helpful to read the papers we wrote, which can be found here in the sidebar.
Gridding data
The gridding of data in Strike.py is highly restricted. The grid must be square and the length of the sides must be a power of 2: 4, 8, 16, 32, ... This is necessary to be able to transform maps into the frequence domain using FFT.
The program Grids.py is not restricted by this and the grid only need to be rectangular. The program has routines for reading a data file, gridding the data, and selecting outlines in the landscape to examine together with the data. There is also a routine to remove a background trend. All this is similar to what is offered in Strike.py
What is added here is the possibility to work with two equally sized grids, add them together or subtract one from the other. Also offered is the the multiplication of a grid by a constant. An example of this being used is in the article Rætur Ísland (The roots of Iceland), see Papers and talks in the sidebar.
Finally there is a routine to read a 6-column file (X, Y, Elevation, Station name, Free-air value, Bouguer value) and reduce to a three column XYZ file for further processing. The Z value optional: Elevation, Free-air value or Bouguer value.
A more detailed account of these routines is given in the manual for Grids. It may also be helpful to read the papers we wrote, which can be found in the sidebar.
Bouguer density
The program Bouguer.py is intended to examine the effect of different Bouguer densities on the results of gravity measurements, that is, the Bouguer values. In the gravity data supplied with these programs, the rock is assigned a density of 2.6 g/cm3 when correcting for topography, but this Bouguer density can easily be changed and the effect examined.
The X and Y coordinates in gravity data are usually longitude and latitude, but for map making it is more convenient to use kilometers. The first routine therefore allows the data to be transformed between these two coordinate systems.
Surfaces that appear equally rough at all scales have a constant fractal dimension, which is a measure of roughness. This applies to topography supported by a rigid crust, but not when distances become so great that the topography is in isostatic equilibrium. The fractal dimension routine estimates the fractal dimension over selected distance ranges. If the input file contains three values, XYZ, then the fractal dimension of Z is calculated, but if it is a gravity data file with 6 values, there is a choice of which value to use.
The bouguer density routine creates new gravity data with a new Bouguer density and saves it. The last two routines evaluate the effect of different Bouguer densities by minimizing
• the correlation between the gravity field and the topography (Nettleton’s method),
• the roughness of the gravity field over the constant fractal distance range.
Choosing the best Bouguer density on the basis of the smallest correlation with topography is based on the assumption that topography and the gravity field are independent factors. That assumption does not hold where the crust is thin and relatively flexible. We therefore developed the method of looking for the lowest fractal dimension of the gravity field, because the topography has a higher fractal dimension, and overestimating its effect increases the fractal dimension of the gravity field.
A more detailed account of these routines is given in the manual for Bouguer. It may also be helpful to read the papers we wrote, which can be found in the sidebar.
Programs & data
The Python code: Strike.py Grids.py Bouguer.py
Download the three programs as a ZIP file
How to set up and run the programs
Download data examples as a ZIP file
Download files with physical map features as a ZIP file
References
- [1]
- [2]
- [3]
- [4]