geoist.pfm package

Submodules

geoist.pfm.euler module

  • EulerDeconv: The classic 3D solution to Euler’s equation for potential fields (Reid et al., 1990). Runs on the whole dataset.
  • EulerDeconvEW: Run Euler deconvolution on an expanding window over the data set and keep the best estimate.
  • EulerDeconvMW: Run Euler deconvolution on a moving window over the data set to produce a set of estimates.

References

Reid, A. B., J. M. Allsop, H. Granser, A. J. Millett, and I. W. Somerton (1990), Magnetic interpretation in three dimensions using Euler deconvolution, Geophysics, 55(1), 80-91, doi:10.1190/1.1442774.

class geoist.pfm.euler.EulerDeconv(x, y, z, field, xderiv, yderiv, zderiv, structural_index)

Bases: geoist.inversion.misfit.Misfit

Classic 3D Euler deconvolution of potential field data.

Follows the formulation of Reid et al. (1990). Performs the deconvolution on the whole data set. For windowed approaches, use EulerDeconvMW (moving window) and EulerDeconvEW (expanding window).

Works on any potential field that satisfies Euler’s homogeneity equation (both gravity and magnetic, assuming simple sources):

\[(x_i - x_0)\dfrac{\partial f_i}{\partial x} + (y_i - y_0)\dfrac{\partial f_i}{\partial y} + (z_i - z_0)\dfrac{\partial f_i}{\partial z} = \eta (b - f_i),\]

in which \(f_i\) is the given potential field observation at point \((x_i, y_i, z_i)\), \(b\) is the base level (a constant shift of the field, like a regional field), \(\eta\) is the structural index, and \((x_0, y_0, z_0)\) are the coordinates of a point on the source (for a sphere, this is the center point).

The Euler deconvolution estimates \((x_0, y_0, z_0)\) and \(b\) given a potential field and its x, y, z derivatives and the structural index. However, this assumes that the sources are ideal (see the table below). We recommend reading Reid and Thurston (2014) for a discussion on what the structural index means and what it does not mean.

Warning

Please read the paper Reid et al. (2014) to avoid doing horrible things with Euler deconvolution. Uieda et al. (2014) offer a practical tutorial using geoist code and show some common misinterpretations.

After Reid et al. (2014), values of the structural index (SI) can be:

Source type SI (Mag) SI (Grav)
Point, sphere 3 2
Line, cylinder, thin bed fault 2 1
Thin sheet edge, thin sill, thin dyke 1 0

Use the fit() method to run the deconvolution. The estimated coordinates \((x_0, y_0, z_0)\) are stored in the estimate_ attribute and the estimated base level \(b\) is stored in baselevel_.

Note

Using structural index of 0 is not supported yet.

Note

The data does not need to be gridded for this! So long as you can calculate the derivatives of non-gridded data (using an Equivalent Layer, for example).

Note

x is North, y is East, and z is down.

Note

Units of the input data (x, y, z, field, derivatives) must be in SI units! Otherwise, the results will be in strange units. Use functions in geoist.pfm.giutils to convert between units.

Parameters:

  • x, y, z : 1d-arrays
    The x, y, and z coordinates of the observation points
  • field : 1d-array
    The potential field measured at the observation points
  • xderiv, yderiv, zderiv : 1d-arrays
    The x-, y-, and z-derivatives of the potential field (measured or calculated) at the observation points
  • index : float
    The structural index of the source

References:

Reid, A. B., J. M. Allsop, H. Granser, A. J. Millett, and I. W. Somerton (1990), Magnetic interpretation in three dimensions using Euler deconvolution, Geophysics, 55(1), 80-91, doi:10.1190/1.1442774.

Reid, A. B., J. Ebbing, and S. J. Webb (2014), Avoidable Euler Errors – the use and abuse of Euler deconvolution applied to potential fields, Geophysical Prospecting, doi:10.1111/1365-2478.12119.

Reid, A., and J. Thurston (2014), The structural index in gravity and magnetic interpretation: Errors, uses, and abuses, GEOPHYSICS, 79(4), J61-J66, doi:10.1190/geo2013-0235.1.

Uieda, L., V. C. Oliveira Jr., and V. C. F. Barbosa (2014), Geophysical tutorial: Euler deconvolution of potential-field data, The Leading Edge, 33(4), 448-450, doi:10.1190/tle33040448.1.

baselevel_
fmt_estimate(p)

Separate the (x, y, z) point coordinates from the baselevel.

Coordinates are stored in estimate_ and a base level is stored in baselevel_.

jacobian(p)
predicted(p)

Calculate the predicted data for a given parameter vector.

Parameters:

  • p : 1d-array or None
    The parameter vector used to calculate the predicted data. If None, will use the current estimate stored in estimate_.

Returns:

  • predicted : 1d-array or list of 1d-arrays
    The predicted data. If this is the sum of 1 or more Misfit instances, will return the predicted data from each of the summed misfits in the order of the sum.
class geoist.pfm.euler.EulerDeconvEW(x, y, z, field, xderiv, yderiv, zderiv, structural_index, center, sizes)

Bases: geoist.pfm.euler.EulerDeconv

Euler deconvolution using an expanding window scheme.

Uses data inside a window of growing size to perform the Euler deconvolution. Keeps the best result, judged by the estimated error.

The deconvolution is performed as in EulerDeconv.

Use the fit() method to produce an estimate. The estimated point is stored in the attribute estimate_ and the base level in baselevel_.

Parameters:

  • x, y, z : 1d-arrays
    The x, y, and z coordinates of the observation points
  • field : 1d-array
    The potential field measured at the observation points
  • xderiv, yderiv, zderiv : 1d-arrays
    The x-, y-, and z-derivatives of the potential field (measured or calculated) at the observation points
  • index : float
    The structural index of the source
  • center : [x, y]
    The x, y coordinates of the center of the expanding windows.
  • sizes : list or 1d-array
    The sizes of the windows.
fit()

Perform the Euler deconvolution with expanding windows.

The estimated point is stored in estimate_, the base level in baselevel_.

class geoist.pfm.euler.EulerDeconvMW(x, y, z, field, xderiv, yderiv, zderiv, structural_index, windows, size, keep=0.2)

Bases: geoist.pfm.euler.EulerDeconv

Solve an Euler deconvolution problem using a moving window scheme.

Uses data inside a window moving to perform the Euler deconvolution. Keeps only a top percentage of the estimates from all windows.

The deconvolution is performed as in EulerDeconv.

Use the fit() method to produce an estimate. The estimated points are stored in estimate_ as a 2D numpy array. Each line in the array is an [x, y, z] coordinate of a point. The base levels are stored in baselevel_.

Parameters:

  • x, y, z : 1d-arrays
    The x, y, and z coordinates of the observation points
  • field : 1d-array
    The potential field measured at the observation points
  • xderiv, yderiv, zderiv : 1d-arrays
    The x-, y-, and z-derivatives of the potential field (measured or calculated) at the observation points
  • index : float
    The structural index of the source
  • windows : (ny, nx)
    The number of windows in the y and x directions
  • size : (dy, dx)
    The size of the windows in the y and x directions
  • keep : float
    Decimal percentage of solutions to keep. Will rank the solutions by increasing error and keep only the first keep percent.
baselevel_
fit()

Perform the Euler deconvolution on a moving window.

The estimated points are stored in estimate_, the base levels in baselevel_.

fmt_estimate(p)

Separate the (x, y, z) point coordinates from the baselevel.

Coordinates are stored in estimate_ and a base level is stored in baselevel_.

geoist.pfm.giconstants module

Author : Steve Chen <chenshi@cea-igp.ac.cn> Affiliation : Institute of Geophysics, CEA. Version : 0.1.0 Copyright : Copyright (C) 2018-2020 GEOIST Development Team. All Rights Reserved. License : Distributed under the MIT License. See LICENSE.txt for more info. Github : https://igp-gravity.github.io/ Description :

Holds all physical constants and unit conversions used in geoist. All modules should import the constants from here! All constants should be in SI, unless otherwise stated!
geoist.pfm.giconstants.CM = 1e-07

Proportionality constant used in the magnetic method in henry/m (SI)

geoist.pfm.giconstants.G = 6.673e-11

The gravitational constant in \(m^3 kg^{-1} s^{-1}\)

geoist.pfm.giconstants.MEAN_EARTH_RADIUS = 6378137.0

The mean earth radius in meters

geoist.pfm.giconstants.PERM_FREE_SPACE = 1.2566370614359173e-06

Permeability of free space in \(N A^{-2}\)

geoist.pfm.giconstants.SI2EOTVOS = 1000000000.0

\(1/s^2 = 10^9\ Eotvos\)

Type:Conversion factor from SI units to Eotvos
geoist.pfm.giconstants.SI2MGAL = 100000.0

\(1\ m/s^2 = 10^5\ mGal\)

Type:Conversion factor from SI units to mGal
geoist.pfm.giconstants.T2NT = 1000000000.0

Conversion factor from tesla to nanotesla

geoist.pfm.giconstants.THERMAL_DIFFUSIVITY = 1e-06

The default thermal diffusivity in \(m^2/s\)

geoist.pfm.giconstants.THERMAL_DIFFUSIVITY_YEAR = 31.5576

The default thermal diffusivity but in \(m^2/year\)

geoist.pfm.giutils module

Author : Steve Chen <chenshi@cea-igp.ac.cn> Affiliation : Institute of Geophysics, CEA. Version : 0.1.0 Copyright : Copyright (C) 2018-2020 GEOIST Development Team. All Rights Reserved. License : Distributed under the MIT License. See LICENSE.txt for more info. Github : https://igp-gravity.github.io/ Description : Miscellaneous utility functions. ——————————————————————————-

class geoist.pfm.giutils.SparseList(size, elements=None)

Bases: object

Store only non-zero elements on an immutable list.

Can iterate over and access elements just like if it were a list.

Parameters:

  • size : int
    Size of the list.
  • elements : dict
    Dictionary used to initialize the list. Keys are the index of the elements and values are their respective values.

Example:

>>> l = SparseList(5)
>>> l[3] = 42.0
>>> print len(l)
5
>>> print l[1], l[3]
0.0 42.0
>>> l[1] += 3.0
>>> for i in l:
...     print i,
0.0 3.0 0.0 42.0 0.0
>>> l2 = SparseList(4, elements={1:3.2, 3:2.8})
>>> for i in l2:
...     print i,
0.0 3.2 0.0 2.8
next()
geoist.pfm.giutils.ang2vec(intensity, inc, dec)

Convert intensity, inclination and declination to a 3-component vector

Note

Coordinate system is assumed to be x->North, y->East, z->Down. Inclination is positive down and declination is measured with respect to x (North).

Parameter:

  • intensity : float or array
    The intensity (norm) of the vector
  • inc : float
    The inclination of the vector (in degrees)
  • dec : float
    The declination of the vector (in degrees)

Returns:

  • vec : array = [x, y, z]
    The vector

Examples:

>>> import numpy
>>> print ang2vec(3, 45, 45)
[ 1.5         1.5         2.12132034]
>>> print ang2vec(numpy.arange(4), 45, 45)
[[ 0.          0.          0.        ]
 [ 0.5         0.5         0.70710678]
 [ 1.          1.          1.41421356]
 [ 1.5         1.5         2.12132034]]
geoist.pfm.giutils.check_hash(fname, known_hash, hash_type='sha256')

Assert that the hash of a file is equal to a known hash.

Useful for checking if a file has changed or been corrupted.

Parameters:

  • fname : string
    The name of the file.
  • known_hash : string
    The known (recorded) hash of the file.
  • hash_type : string
    What kind of hash is it. Can be anything defined in Python’s hashlib.

Raises:

  • AssertionError if the hash is different from the known hash.
geoist.pfm.giutils.contaminate(data, stddev, percent=False, return_stddev=False, seed=None)

Add pseudorandom gaussian noise to an array.

Noise added is normally distributed with zero mean.

Parameters:

  • data : array or list of arrays
    Data to contaminate
  • stddev : float or list of floats
    Standard deviation of the Gaussian noise that will be added to data
  • percent : True or False
    If True, will consider stddev as a decimal percentage and the standard deviation of the Gaussian noise will be this percentage of the maximum absolute value of data
  • return_stddev : True or False
    If True, will return also the standard deviation used to contaminate data
  • seed : None or int
    Seed used to generate the pseudo-random numbers. If None, will use a different seed every time. Use the same seed to generate the same random sequence to contaminate the data.

Returns:

if return_stddev is False:

  • contam : array or list of arrays
    The contaminated data array

else:

  • results : list = [contam, stddev]
    The contaminated data array and the standard deviation used to contaminate it.

Examples:

>>> import numpy as np
>>> data = np.ones(5)
>>> noisy = contaminate(data, 0.1, seed=0)
>>> print noisy
[ 1.03137726  0.89498775  0.95284582  1.07906135  1.04172782]
>>> noisy, std = contaminate(data, 0.05, seed=0, percent=True,
...                          return_stddev=True)
>>> print std
0.05
>>> print noisy
[ 1.01568863  0.94749387  0.97642291  1.03953067  1.02086391]
>>> data = [np.zeros(5), np.ones(3)]
>>> noisy = contaminate(data, [0.1, 0.2], seed=0)
>>> print noisy[0]
[ 0.03137726 -0.10501225 -0.04715418  0.07906135  0.04172782]
>>> print noisy[1]
[ 0.81644754  1.20192079  0.98163167]
geoist.pfm.giutils.dircos(inc, dec)

Returns the 3 coordinates of a unit vector given its inclination and declination.

Note

Coordinate system is assumed to be x->North, y->East, z->Down. Inclination is positive down and declination is measured with respect to x (North).

Parameter:

  • inc : float
    The inclination of the vector (in degrees)
  • dec : float
    The declination of the vector (in degrees)

Returns:

  • vect : list = [x, y, z]
    The unit vector
geoist.pfm.giutils.eotvos2si(value)

Convert a value from Eotvos to SI units.

Parameters:

  • value : number or array
    The value in Eotvos

Returns:

  • value : number or array
    The value in SI
geoist.pfm.giutils.gaussian(x, mean, std)

Non-normalized Gaussian function

\[G(x,\bar{x},\sigma) = \exp\left(-\frac{(x-\bar{x})^2}{\sigma^2} \right)\]

Parameters:

  • x : float or array
    Values at which to calculate the Gaussian function
  • mean : float
    The mean of the distribution \(\bar{x}\)
  • std : float
    The standard deviation of the distribution \(\sigma\)

Returns:

  • gauss : array
    Gaussian function evaluated at x
geoist.pfm.giutils.gaussian2d(x, y, sigma_x, sigma_y, x0=0, y0=0, angle=0.0)

Non-normalized 2D Gaussian function

Parameters:

  • x, y : float or arrays
    Coordinates at which to calculate the Gaussian function
  • sigma_x, sigma_y : float
    Standard deviation in the x and y directions
  • x0, y0 : float
    Coordinates of the center of the distribution
  • angle : float
    Rotation angle of the gaussian measure from the x axis (north) growing positive to the east (positive y axis)

Returns:

  • gauss : array
    Gaussian function evaluated at x, y
geoist.pfm.giutils.mgal2si(value)

Convert a value from mGal to SI units.

Parameters:

  • value : number or array
    The value in mGal

Returns:

  • value : number or array
    The value in SI
geoist.pfm.giutils.nt2si(value)

Convert a value from nanoTesla to SI units.

Parameters:

  • value : number or array
    The value in nanoTesla

Returns:

  • value : number or array
    The value in SI
geoist.pfm.giutils.safe_diagonal(matrix)

Get the diagonal of a matrix using the appropriate method.

Parameters:

  • matrix : 2d-array, matrix, sparse matrix
    The matrix…

Returns:

  • diag : 1d-array
    A numpy array with the diagonal of the matrix
geoist.pfm.giutils.safe_dot(a, b)

Make the dot product using the appropriate method.

If a and b are dense, will use numpy.dot(). If either is sparse (from scipy.sparse) will use the multiplication operator (i.e., *).

Parameters:

  • a, b : array or matrix
    The vectors/matrices to take the dot product of.

Returns:

  • prod : array or matrix
    The dot product of a and b
geoist.pfm.giutils.safe_inverse(matrix)

Calculate the inverse of a matrix using an apropriate algorithm.

Uses the standard numpy.linalg.inv() if matrix is dense. If it is sparse (from scipy.sparse) then will use scipy.sparse.linalg.inv().

Parameters:

  • matrix : 2d-array
    The matrix

Returns:

  • inverse : 2d-array
    The inverse of matrix
geoist.pfm.giutils.safe_solve(matrix, vector)

Solve a linear system using an apropriate algorithm.

Uses the standard numpy.linalg.solve() if both matrix and vector are dense.

If any of the two is sparse (from scipy.sparse) then will use the Conjugate Gradient Method (scipy.sparse.cgs()).

Parameters:

  • matrix : 2d-array
    The matrix defining the linear system
  • vector : 1d or 2d-array
    The right-side vector of the system

Returns:

  • solution : 1d or 2d-array
    The solution of the linear system
geoist.pfm.giutils.si2eotvos(value)

Convert a value from SI units to Eotvos.

Parameters:

  • value : number or array
    The value in SI

Returns:

  • value : number or array
    The value in Eotvos
geoist.pfm.giutils.si2mgal(value)

Convert a value from SI units to mGal.

Parameters:

  • value : number or array
    The value in SI

Returns:

  • value : number or array
    The value in mGal
geoist.pfm.giutils.si2nt(value)

Convert a value from SI units to nanoTesla.

Parameters:

  • value : number or array
    The value in SI

Returns:

  • value : number or array
    The value in nanoTesla
geoist.pfm.giutils.sph2cart(lon, lat, height)

Convert spherical coordinates to Cartesian geocentric coordinates.

Parameters:

  • lon, lat, height : floats
    Spherical coordinates. lon and lat in degrees, height in meters. height is the height above mean Earth radius.

Returns:

  • x, y, z : floats
    Converted Cartesian coordinates
geoist.pfm.giutils.vec2ang(vector)

Convert a 3-component vector to intensity, inclination and declination.

Note

Coordinate system is assumed to be x->North, y->East, z->Down. Inclination is positive down and declination is measured with respect to x (North).

Parameter:

  • vector : array = [x, y, z]
    The vector

Returns:

  • [intensity, inclination, declination] : floats
    The intensity, inclination and declination (in degrees)

Examples:

>>> s = vec2ang([1.5, 1.5, 2.121320343559643])
>>> print "%.3f %.3f %.3f" % tuple(s)
3.000 45.000 45.000

geoist.pfm.grdio module

Name : grdio.py Created on : 2018/11/24 08:57 Author : Steve Chen <chenshi@cea-igp.ac.cn> Affiliation : Institute of Geophysics, CEA. Version : 0.1.0 Copyright : Copyright (C) 2018-2020 GEOIST Development Team. All Rights Reserved. License : Distributed under the MIT License. See LICENSE.txt for more info. Github : https://igp-gravity.github.io/ Description : Application for processing grid data of potential dataset.

class geoist.pfm.grdio.grddata

Bases: object

Grid Data Object .. attribute:: data

array to contain raster data

type:numpy masked array
xmin

min value X coordinate of raster grid

Type:float
ymin

min value Y coordinate of raster grid

Type:float
xdim

x-dimension of grid cell

Type:float
ydim

y-dimension of grid cell

Type:float
typeofdata

number of datatype

Type:int
dataname

data name or id

Type:str
rows

number of rows for each raster grid/band

Type:int
cols

number of columns for each raster grid/band

Type:int
nullvalue

grid null or nodata value

Type:float
norm

normalized data

Type:dictionary
gtr

projection information

Type:tuple
wkt

projection information

Type:str
units

description of units to be used with color bars

Type:str
export_ascii(fname)

Export Ascii file

Parameters:data (grid Data) – dataset to export
export_surfer(fname, flag=True, file_format='binary')

Export a surfer grid

Parameters:
  • fname (filename of grid dataset to export) –
  • flag (True - Output Grid Grid) – False - Output Bak Grid Grid
  • file_format (binary/b - output binary format) – ascii/a - output ascii format
fill_nulls(method='nearest')

Fill in the NaNs or masked values on interpolated points using nearest neighbors. method=’nearest’ or ‘linear’ or ‘cubic’

grd2xyz(flag=True)

Return x,y,z 1-D array data from 2-D grid array.

Parameters:flag – True - Output Grid Grid False - Output Bak Grid Grid
Returns:x,y,z 1-D array data
load_ascii(fname, dtype='float64')

Load Ascii file

Parameters:data (grid Data) – dataset to export
load_grd(fname, *args, **kwargs)
load_surfer(fname, *args, **kwargs)

Read data from a Surfer grid file.

Parameters:

  • fname : str
    Name of the Surfer grid file
  • dtype : numpy dtype object or string
    The type of variable used for the data. Default is numpy.float64 for ascii data and is ‘=f’ for binary data. Use numpy.float32 if the data are large and precision is not an issue.
  • header_format : header format (excluding the leading ‘DSBB’) following
    the convention of the struct module. Only used for binary data.

Returns:

geoist.pfm.grdio.regular(area, shape, z=None)

Create a regular grid.

The x directions is North-South and y East-West. Imagine the grid as a matrix with x varying in the lines and y in columns.

Returned arrays will be flattened to 1D with numpy.ravel.

Parameters:

  • area
    (x1, x2, y1, y2): Borders of the grid
  • shape
    Shape of the regular grid, ie (nx, ny).
  • z
    Optional. z coordinate of the grid points. If given, will return an array with the value z.

Returns:

  • [x, y]
    Numpy arrays with the x and y coordinates of the grid points
  • [x, y, z]
    If z given. Numpy arrays with the x, y, and z coordinates of the grid points

Examples:

>>> x, y = regular((0, 10, 0, 5), (5, 3))
>>> print(x)
[  0.    0.    0.    2.5   2.5   2.5   5.    5.    5.    7.5   7.5   7.5
  10.   10.   10. ]
>>> print(x.reshape((5, 3)))
[[  0.    0.    0. ]
 [  2.5   2.5   2.5]
 [  5.    5.    5. ]
 [  7.5   7.5   7.5]
 [ 10.   10.   10. ]]
geoist.pfm.grdio.spacing(area, shape)

Returns the spacing between grid nodes

Parameters:

  • area
    (x1, x2, y1, y2): Borders of the grid
  • shape
    Shape of the regular grid, ie (nx, ny).

Returns:

  • [dx, dy]
    Spacing the y and x directions

Examples:

>>> print(spacing((0, 10, 0, 20), (11, 11)))
[1.0, 2.0]
>>> print(spacing((0, 10, 0, 20), (11, 21)))
[1.0, 1.0]
>>> print(spacing((0, 10, 0, 20), (5, 21)))
[2.5, 1.0]
>>> print(spacing((0, 10, 0, 20), (21, 21)))
[0.5, 1.0]

geoist.pfm.igrf module

Name : igrf.py Created on : 2018/09/11 17:00 Author : Steve Chen <chenshi@cea-igp.ac.cn> Affiliation : Institute of Geophysics, CEA. Version : 0.1.0 Copyright : Copyright (C) 2018-2020 GEOIST Development Team. All Rights Reserved. License : Distributed under the MIT License. See LICENSE.txt for more info. Github : https://igp-gravity.github.io/ Description : Application for ***.

class geoist.pfm.igrf.IGRF(parent=None)

Bases: object

IGRF field calculation

This produces two datasets. The first is an IGRF dataset for the area of interest, defined by some input magnetic dataset. The second is the IGRF corrected form of that input magnetic dataset.

To do this, the input dataset must be reprojected from its local projection to degrees, where the IGRF correction will take place. This is done within this class.

parent

reference to the parent routine

Type:parent
indata

dictionary of input datasets

Type:dictionary
outdata

dictionary of output datasets

Type:dictionary
Parameters:
  • altmin (Double) – Minimum height of selected model.
  • altmax (Double array) – array of MAXMOD Maximum height of model.
  • maxalt (Double) – Maximum height of selected model.
  • d (float) – Declination of the field from the geographic north (deg).
  • sdate (float) – start date inputted
  • ddot (float) – annual rate of change of decl. (arc-min/yr)
  • alt (float) – altitude above WGS84 Ellipsoid
  • epoch (list) – list of MAXMOD epoch of model.
  • latitude (float) – Latitude.
  • longitude (float) – Longitude.
  • gh (numpy array) – Schmidt quasi-normal internal spherical harmonic coeff. Schmidt quasi-normal internal spherical harmonic coeff. Coefficients of resulting model. Coefficients of rate of change model.
  • i (float) – Inclination (deg).
  • idot (float) – Rate of change of i (arc-min/yr).
  • igdgc (int) – Flag for geodetic or geocentric coordinate choice.
  • irec_pos (int array) – array of MAXMOD Record counter for header
  • fileline (int) – Current line in file (for errors)
  • max1 (list, int) – array of MAXMOD Main field coefficient.
  • max2 (list, int) – array of MAXMOD Secular variation coefficient.
  • max3 (list, int) – array of MAXMOD Acceleration coefficient.
  • minyr (float) – Min year of all models
  • maxyr (float) – Max year of all models
  • yrmax (list, float) – array of MAXMOD Max year of model.
  • yrmin (list, float) – array of MAXMOD Min year of model.
dihf(gh)

Computes the geomagnetic d, i, h, and f from x, y, and z.

FORTRAN : A. Zunde, USGS, MS 964, box 25046 Federal Center, Denver,
CO. 80225
C : C. H. Shaffer, Lockheed Missiles and Space Company, Sunnyvale CA
Parameters:
  • x (float) – northward component
  • y (float) – eastward component
  • z (float) – vertically-downward component
Returns:

  • d (float) – declination
  • i (float) – inclination
  • h (float) – horizontal intensity
  • f (float) – total intensity
extrapsh(date, dte1, nmax1, nmax2, gh)

Extrapolates linearly a spherical harmonic model with a rate-of-change model.

FORTRAN : A. Zunde, USGS, MS 964, box 25046 Federal Center, Denver,
CO. 80225
C : C. H. Shaffer, Lockheed Missiles and Space Company, Sunnyvale CA
Parameters:
  • date (float) – date of resulting model (in decimal year)
  • dte1 (float) – date of base model
  • nmax1 (int) – maximum degree and order of base model
  • gh (numpy array) – Schmidt quasi-normal internal spherical harmonic coefficients of base model and rate-of-change model
  • nmax2 (int) – maximum degree and order of rate-of-change model
Returns:

  • gh (numpy array) – Schmidt quasi-normal internal spherical harmonic coefficients
  • nmax (int) – maximum degree and order of resulting model
getshc(file, iflag, strec, nmax_of_gh, gh)

Reads spherical harmonic coefficients from the specified model into an array.

FORTRAN: Bill Flanagan, NOAA CORPS, DESDIS, NGDC, 325 Broadway,
Boulder CO. 80301
C: C. H. Shaffer, Lockheed Missiles and Space Company, Sunnyvale CA
Parameters:
  • file (file) – reference to a file object
  • iflag – Flag for SV equal to 1 or not equal to 1 for designated read statements
  • strec (int) – Starting record number to read from model
  • nmax_of_gh (int) – Maximum degree and order of model
Returns:

gh – Schmidt quasi-normal internal spherical harmonic coefficients
Return type:

list
interpsh(date, dte1, nmax1, dte2, nmax2, gh)

Interpolates linearly, in time, between two spherical harmonic models.

FORTRAN : A. Zunde, USGS, MS 964, box 25046 Federal Center, Denver,
CO. 80225
C : C. H. Shaffer, Lockheed Missiles and Space Company, Sunnyvale CA
Parameters:
  • date (float) – date of resulting model (in decimal year)
  • dte1 (float) – date of earlier model
  • nmax1 (int) – maximum degree and order of earlier model
  • gh (numpy array) – Schmidt quasi-normal internal spherical harmonic coefficients of earlier model and internal model
  • dte2 (float) – date of later model
  • nmax2 (int) – maximum degree and order of later model
Returns:

  • gh (numpy array) – coefficients of resulting model
  • nmax (int) – maximum degree and order of resulting model
pnt(latitude, longitude, alt, magdata, Model='IGRF12')

Settings Dialog. This is the main entrypoint into this routine. It also contains the main IGRF code.

shval3(igdgc, flat, flon, elev, nmax, gh)

Calculates field components from spherical harmonic (sh) models.

Based on subroutine ‘igrf’ by D. R. Barraclough and S. R. C. Malin, report no. 71/1, institute of geological sciences, U.K.

FORTRAN : Norman W. Peddie, USGS, MS 964, box 25046 Federal Center,
Denver, CO. 80225
C : C. H. Shaffer, Lockheed Missiles and Space Company, Sunnyvale CA
Parameters:
  • igdgc (int) – indicates coordinate system used set equal to 1 if geodetic, 2 if geocentric
  • latitude (float) – north latitude, in degrees
  • longitude (float) – east longitude, in degrees
  • elev (float) – WGS84 altitude above ellipsoid (igdgc=1), or radial distance from earth’s center (igdgc=2)
  • a2,b2 (float) – squares of semi-major and semi-minor axes of the reference spheroid used for transforming between geodetic and geocentric coordinates or components
  • nmax (int) – maximum degree and order of coefficients
Returns:

  • x (float) – northward component
  • y (float) – eastward component
  • z (float) – vertically-downward component

geoist.pfm.normgra module

Reference ellipsoids

This module uses instances of ReferenceEllipsoid to store the physical constants of ellipsoids. To create a new ellipsoid, just instantiate ReferenceEllipsoid and give it the semimajor axis a, the flattening f, the geocentric gravitational constant GM, and the angular velocity omega. All other quantities, like the gravity at the poles and equator, eccentricities, etc, are computed by the class from these 4 parameters.

Available ellipsoids:

  • WGS84 (values taken from Hofmann-Wellenhof and Moritz, 2005):

    >>> from geoist.pfm.normgra import WGS84
>>> print(WGS84.name)
World Geodetic System 1984
>>> print('{:.0f}'.format(WGS84.a))
6378137
>>> print('{:.17f}'.format(WGS84.f))
0.00335281066474748
>>> print('{:.10g}'.format(WGS84.GM))
3.986004418e+14
>>> print('{:.7g}'.format(WGS84.omega))
7.292115e-05
>>> print('{:.4f}'.format(WGS84.b))
6356752.3142
>>> print('{:.8f}'.format(WGS84.E)) # Linear eccentricity
521854.00842339
>>> print('{:.15f}'.format(WGS84.e_prime)) # second eccentricity
0.082094437949696
>>> print('{:.10f}'.format(WGS84.gamma_a)) # gravity at the equator
9.7803253359
>>> print('{:.11f}'.format(WGS84.gamma_b)) # gravity at the pole
9.83218493786
>>> print('{:.14f}'.format(WGS84.m))
0.00344978650684

Normal gravity

  • gamma_somigliana(): compute the normal gravity using the Somigliana formula (Hofmann-Wellenhof and Moritz, 2005). Calculated on the ellipsoid.
  • gamma_somigliana_free_air(): compute normal gravity at a height using the Somigliana formula plus the free-air correction \(-0.3086H\ mGal/m\).
  • gamma_closed_form(): compute normal gravity using the closed form expression from Li and Gotze (2001). Can compute anywhere (on, above, under the ellipsoid).

Bouguer

  • bouguer_plate(): compute the gravitational attraction of an infinite plate (Bouguer plate). Calculated on top of the plate.

You can use pfm.prism and pfm.tesseroid to calculate the terrain effect for a better correction.

References

Hofmann-Wellenhof, B. and H. Moritz, 2005, Physical Geodesy, Springer-Verlag Wien, ISBN-13: 978-3-211-23584-3

Li, X. and H. J. Gotze, 2001, Tutorial: Ellipsoid, geoid, gravity, geodesy, and geophysics, Geophysics, 66(6), p. 1660-1668, doi: 10.1190/1.1487109

class geoist.pfm.normgra.ReferenceEllipsoid(name, a, f, GM, omega)

Bases: object

A generic reference ellipsoid.

It stores the physical constants defining the ellipsoid and has properties for computing other (derived) quantities.

All quantities are expected and returned in SI units.

Parameters:

  • a : float
    The semimajor axis (the largest one, at the equator). In meters.
  • f : float
    The flattening. Adimensional.
  • GM : float
    The geocentric gravitational constant of the earth, including the atmosphere. In \(m^3 s^{-2}\).
  • omega : float
    The angular velocity of the earth. In \(rad s^{-1}\).
E

Linear eccentricity

GM

Geocentric gravitational constant (including the atmosphere)

a

Semimajor axis

b

Semiminor axis

e_prime

Second eccentricity

f

Flattening

gamma_a

Normal gravity at the equator

gamma_b

Normal gravity at the poles

m

\(\omega^2 a^2 b / (GM)\)

omega

Angular velocity

geoist.pfm.normgra.bouguer_plate(topography, density_rock=2670, density_water=1040)

Calculate the gravitational effect of an infinite Bouguer plate.

Note

The effect is calculated on top of the plate.

Uses the famous \(g_{BG} = 2 \pi G \rho t\) formula, where t is the height of the topography. On water (i.e., t < 0), uses \(g_{BG} = 2 \pi G (\rho_{water} - \rho_{rock})\times (-t)\).

Parameters:

  • topography : float or numpy array
    The height of topography (in meters).
  • density_rock : float
    The density of crustal rocks
  • density_water : float
    The density of ocean water

Returns:

  • g_bouguer : float or array
    The computed gravitational effect of the Bouguer plate
geoist.pfm.normgra.gamma_closed_form(latitude, height, ellipsoid=<geoist.pfm.normgra.ReferenceEllipsoid object>)

Calculate normal gravity at a height using the closed form expression of Li and Gotze (2001).

Parameters:

  • latitude : float or numpy array

    The latitude where the normal gravity will be computed (in degrees)

  • height : float or numpy array

    The height of computation (in meters). Should be ellipsoidal (geometric) heights for geophysical purposes.

  • ellipsoid : ReferenceEllipsoid

    The reference ellipsoid used.

Returns:

  • gamma : float or numpy array
    The computed normal gravity (in mGal).
geoist.pfm.normgra.gamma_somigliana(latitude, ellipsoid=<geoist.pfm.normgra.ReferenceEllipsoid object>)

Calculate the normal gravity by using Somigliana’s formula.

This formula computes normal gravity on the ellipsoid (height = 0).

Parameters:

  • latitude : float or numpy array

    The latitude where the normal gravity will be computed (in degrees)

  • ellipsoid : ReferenceEllipsoid

    The reference ellipsoid used.

Returns:

  • gamma : float or numpy array
    The computed normal gravity (in mGal).
geoist.pfm.normgra.gamma_somigliana_free_air(latitude, height, ellipsoid=<geoist.pfm.normgra.ReferenceEllipsoid object>)

Calculate the normal gravity at a height using Somigliana’s formula and the free-air correction.

Parameters:

  • latitude : float or numpy array

    The latitude where the normal gravity will be computed (in degrees)

  • height : float or numpy array

    The height of computation (in meters). Should be ellipsoidal (geometric) heights for geophysical purposes.

  • ellipsoid : ReferenceEllipsoid

    The reference ellipsoid used.

Returns:

  • gamma : float or numpy array
    The computed normal gravity (in mGal).

geoist.pfm.pftrans module

Potential field transformations, like upward continuation and derivatives. .. note:: Most, if not all, functions here required gridded data.

geoist.pfm.pftrans.derivx(x, y, data, shape, order=1, method='fd')

Calculate the derivative of a potential field in the x direction.

Note

Requires gridded data.

Warning

If the data is not in SI units, the derivative will be in strange units! I strongly recommend converting the data to SI before calculating the derivative (use one of the unit conversion functions of geoist.utils). This way the derivative will be in SI units and can be easily converted to what unit you want.

Parameters:

  • x, y : 1D-arrays
    The x and y coordinates of the grid points
  • data : 1D-array
    The potential field at the grid points
  • shape : tuple = (nx, ny)
    The shape of the grid
  • order : int
    The order of the derivative
  • method : string
    The method used to calculate the derivatives. Options are: 'fd' for central finite-differences (more stable) or 'fft' for the Fast Fourier Transform.

Returns:

  • deriv : 1D-array
    The derivative
geoist.pfm.pftrans.derivy(x, y, data, shape, order=1, method='fd')

Calculate the derivative of a potential field in the y direction.

Note

Requires gridded data.

Warning

If the data is not in SI units, the derivative will be in strange units! I strongly recommend converting the data to SI before calculating the derivative (use one of the unit conversion functions of geoist.utils). This way the derivative will be in SI units and can be easily converted to what unit you want.

Parameters:

  • x, y : 1D-arrays
    The x and y coordinates of the grid points
  • data : 1D-array
    The potential field at the grid points
  • shape : tuple = (nx, ny)
    The shape of the grid
  • order : int
    The order of the derivative
  • method : string
    The method used to calculate the derivatives. Options are: 'fd' for central finite-differences (more stable) or 'fft' for the Fast Fourier Transform.

Returns:

  • deriv : 1D-array
    The derivative
geoist.pfm.pftrans.derivz(x, y, data, shape, order=1, method='fft')

Calculate the derivative of a potential field in the z direction.

Note

Requires gridded data.

Warning

If the data is not in SI units, the derivative will be in strange units! I strongly recommend converting the data to SI before calculating the derivative (use one of the unit conversion functions of geoist.utils). This way the derivative will be in SI units and can be easily converted to what unit you want.

Parameters:

  • x, y : 1D-arrays
    The x and y coordinates of the grid points
  • data : 1D-array
    The potential field at the grid points
  • shape : tuple = (nx, ny)
    The shape of the grid
  • order : int
    The order of the derivative
  • method : string
    The method used to calculate the derivatives. Options are: 'fft' for the Fast Fourier Transform.

Returns:

  • deriv : 1D-array
    The derivative
geoist.pfm.pftrans.power_density_spectra(x, y, data, shape)

Calculates the Power Density Spectra of a 2D gridded potential field through the FFT:

\[\Phi_{\Delta T}(k_x, k_y) = | F\left{\Delta T \right}(k_x, k_y) |^2\]

Note

Requires gridded data.

Note

x, y, z and height should be in meters.

Parameters:

  • x, y : 1D-arrays
    The x and y coordinates of the grid points
  • data : 1D-array
    The potential field at the grid points
  • shape : tuple = (nx, ny)
    The shape of the grid

Returns:

  • kx, ky : 2D-arrays
    The wavenumbers of each Power Density Spectra point
  • pds : 2D-array
    The Power Density Spectra of the data
geoist.pfm.pftrans.radial_average_spectrum(kx, ky, pds, max_radius=None, ring_width=None)

Calculates the average of the Power Density Spectra points that falls inside concentric rings built around the origin of the wavenumber coordinate system with constant width.

The width of the rings and the inner radius of the biggest ring can be changed by setting the optional parameters ring_width and max_radius, respectively.

Note

To calculate the radially averaged power density spectra use the outputs of the function power_density_spectra as input of this one.

Parameters:

  • kx, ky : 2D-arrays
    The wavenumbers arrays in the x and y directions
  • pds : 2D-array
    The Power Density Spectra
  • max_radius : float (optional)
    Inner radius of the biggest ring. By default it’s set as the minimum of kx.max() and ky.max(). Making it smaller leaves points outside of the averaging, and making it bigger includes points nearer to the boundaries.
  • ring_width : float (optional)
    Width of the rings. By default it’s set as the largest value of \(\Delta k_x\) and \(\Delta k_y\), being them the equidistances of the kx and ky arrays. Making it bigger gives more populated averages, and making it smaller lowers the ammount of points per ring (use it carefully).

Returns:

  • k_radial : 1D-array
    Wavenumbers of each Radially Averaged Power Spectrum point. Also, the inner radius of the rings.
  • pds_radial : 1D array
    Radially Averaged Power Spectrum
geoist.pfm.pftrans.reduce_to_pole(x, y, data, shape, inc, dec, sinc, sdec)

Reduce total field magnetic anomaly data to the pole.

The reduction to the pole if a phase transformation that can be applied to total field magnetic anomaly data. It “simulates” how the data would be if both the Geomagnetic field and the magnetization of the source were vertical (\(90^\circ\) inclination) (Blakely, 1996).

This functions performs the reduction in the frequency domain (using the FFT). The transform filter is (in the frequency domain):

\[RTP(k_x, k_y) = \frac{|k|}{ a_1 k_x^2 + a_2 k_y^2 + a_3 k_x k_y + i|k|(b_1 k_x + b_2 k_y)}\]

in which \(k_x\) and \(k_y\) are the wave-numbers in the x and y directions and

\[\begin{split}|k| = \sqrt{k_x^2 + k_y^2} \\ a_1 = m_z f_z - m_x f_x \\ a_2 = m_z f_z - m_y f_y \\ a_3 = -m_y f_x - m_x f_y \\ b_1 = m_x f_z + m_z f_x \\ b_2 = m_y f_z + m_z f_y\end{split}\]

\(\mathbf{m} = (m_x, m_y, m_z)\) is the unit-vector of the total magnetization of the source and \(\mathbf{f} = (f_x, f_y, f_z)\) is the unit-vector of the Geomagnetic field.

Note

Requires gridded data.

Warning

The magnetization direction of the anomaly source is crucial to the reduction-to-the-pole. Wrong values of *sinc* and *sdec* will lead to a wrong reduction.

Parameters:

  • x, y : 1d-arrays
    The x, y, z coordinates of each data point.
  • data : 1d-array
    The total field anomaly data at each point.
  • shape : tuple = (nx, ny)
    The shape of the data grid
  • inc, dec : floats
    The inclination and declination of the inducing Geomagnetic field
  • sinc, sdec : floats
    The inclination and declination of the total magnetization of the anomaly source. The total magnetization is the vector sum of the induced and remanent magnetization. If there is only induced magnetization, use the inc and dec of the Geomagnetic field.

Returns:

  • rtp : 1d-array
    The data reduced to the pole.

References:

Blakely, R. J. (1996), Potential Theory in Gravity and Magnetic Applications, Cambridge University Press.

geoist.pfm.pftrans.tga(x, y, data, shape, method='fd')

Calculate the total gradient amplitude (TGA).

This the same as the 3D analytic signal of Roest et al. (1992), but we prefer the newer, more descriptive nomenclature suggested by Reid (2012).

The TGA is defined as the amplitude of the gradient vector of a potential field \(T\) (e.g. the magnetic total field anomaly):

\[TGA = \sqrt{ \left(\frac{\partial T}{\partial x}\right)^2 + \left(\frac{\partial T}{\partial y}\right)^2 + \left(\frac{\partial T}{\partial z}\right)^2 }\]

Note

Requires gridded data.

Warning

If the data is not in SI units, the derivatives will be in strange units and so will the total gradient amplitude! I strongly recommend converting the data to SI before calculating the TGA is you need the gradient in Eotvos (use one of the unit conversion functions of geoist.utils).

Parameters:

  • x, y : 1D-arrays
    The x and y coordinates of the grid points
  • data : 1D-array
    The potential field at the grid points
  • shape : tuple = (nx, ny)
    The shape of the grid
  • method : string
    The method used to calculate the horizontal derivatives. Options are: 'fd' for finite-difference (more stable) or 'fft' for the Fast Fourier Transform. The z derivative is always calculated by FFT.

Returns:

  • tga : 1D-array
    The amplitude of the total gradient

References:

Reid, A. (2012), Forgotten truths, myths and sacred cows of Potential Fields Geophysics - II, in SEG Technical Program Expanded Abstracts 2012, pp. 1-3, Society of Exploration Geophysicists.

Roest, W., J. Verhoef, and M. Pilkington (1992), Magnetic interpretation using the 3-D analytic signal, GEOPHYSICS, 57(1), 116-125, doi:10.1190/1.1443174.

geoist.pfm.pftrans.tilt(x, y, data, shape, xderiv=None, yderiv=None, zderiv=None)

Calculates the potential field tilt, as defined by Miller and Singh (1994)

\[tilt(f) = tan^{-1}\left( \frac{ \frac{\partial T}{\partial z}}{ \sqrt{\frac{\partial T}{\partial x}^2 + \frac{\partial T}{\partial y}^2}} \right)\]

When used on magnetic total field anomaly data, works best if the data is reduced to the pole.

It’s useful to plot the zero contour line of the tilt to represent possible outlines of the source bodies. Use matplotlib’s pyplot.contour or pyplot.tricontour for this.

Note

Requires gridded data if xderiv, yderiv and zderiv are not given.

Parameters:

  • x, y : 1D-arrays
    The x and y coordinates of the grid points
  • data : 1D-array
    The potential field at the grid points
  • shape : tuple = (nx, ny)
    The shape of the grid. Ignored if xderiv, yderiv and zderiv are given.
  • xderiv : 1D-array or None
    Optional. Values of the derivative in the x direction. If None, will calculated using the default options of derivx()
  • yderiv : 1D-array or None
    Optional. Values of the derivative in the y direction. If None, will calculated using the default options of derivy()
  • zderiv : 1D-array or None
    Optional. Values of the derivative in the z direction. If None, will calculated using the default options of derivz()

Returns:

  • tilt : 1D-array
    The tilt angle of the total field in radians.

References:

Miller, Hugh G, and Vijay Singh. 1994. “Potential Field Tilt — a New Concept for Location of Potential Field Sources.” Journal of Applied Geophysics 32 (2–3): 213-17. doi:10.1016/0926-9851(94)90022-1.

geoist.pfm.pftrans.upcontinue(x, y, data, shape, height)

Upward continuation of potential field data.

Calculates the continuation through the Fast Fourier Transform in the wavenumber domain (Blakely, 1996):

\[F\{h_{up}\} = F\{h\} e^{-\Delta z |k|}\]

and then transformed back to the space domain. \(h_{up}\) is the upward continue data, \(\Delta z\) is the height increase, \(F\) denotes the Fourier Transform, and \(|k|\) is the wavenumber modulus.

Note

Requires gridded data.

Note

x, y, z and height should be in meters.

Note

It is not possible to get the FFT of a masked grid. The default geoist.gridder.interp() call using minimum curvature will not be suitable. Use extrapolate=True or algorithm='nearest' to get an unmasked grid.

Parameters:

  • x, y : 1D-arrays
    The x and y coordinates of the grid points
  • data : 1D-array
    The potential field at the grid points
  • shape : tuple = (nx, ny)
    The shape of the grid
  • height : float
    The height increase (delta z) in meters.

Returns:

  • cont : array
    The upward continued data

References:

Blakely, R. J. (1996), Potential Theory in Gravity and Magnetic Applications, Cambridge University Press.

geoist.pfm.tide module

Name : tide.py Created on : 2018/11/03 17:00 Author : Steve Chen <chenshi@cea-igp.ac.cn> Affiliation : Institute of Geophysics, CEA. Version : 0.1.0 Copyright : Copyright (C) 2018-2020 GEOIST Development Team. All Rights Reserved. License : Distributed under the MIT License. See LICENSE.txt for more info. Github : https://igp-gravity.github.io/ Description : code for UTC time, 1.16 don’t multiply.

class geoist.pfm.tide.TideModel

Bases: object

Class to encapsulate the Longman 1959 tide model.

calculate_julian_century(timestamp)

Calculate the julian century and hour.

Take a datetime object and calculate the decimal Julian century and floating point hour. This is in reference to noon on December 31, 1899 as stated in the Longman paper.

Parameters:timestamp (datetime) – Time stamp to convert
Returns:Julian century and hour
Return type:float, float
plot()

Plot the model results.

Make a simple plot of the gravitational tide results from the model run.

run_model()

Run the model for a range of times.

Runs the tidal model beginning at start_time with time steps of increment seconds for days.

solve_longman(lat, lon, alt, time)

Solve the tide model.

Given the location and datetime object, computes the current gravitational tide and associated quantities. Latitude and longitude and in the traditional decimal notation, altitude is in meters, time is a datetime object.

Parameters:
  • lat (float) – latitude (in degrees)
  • lon (float) – longitude (in degrees)
  • alt (float) – altitude (in meters)
  • time (datetime) – time at which to solve the model
Returns:

lunar, solar, and total gravitational tides
Return type:

float, float, float
write(fname)

Write model results to file.

Write results out of a file for later analysis or reading into another method for analysis/correction of data.

Parameters:fname (string) – name of file to save

Module contents