pygmi.pfmod.grvmag3d

Gravity and magnetic field calculations.

This uses the following algorithms:

References

Singh, B., Guptasarma, D., 2001. New method for fast computation of gravity and magnetic anomalies from arbitrary polyhedral. Geophysics 66, 521-526.

Blakely, R.J., 1996. Potential Theory in Gravity and Magnetic Applications, 1st edn. Cambridge University Press, Cambridge, UK, 441 pp. 200-201

Classes

PlotTest	Plot Raster Class.
GravMag	The GravMag class holds generic magnetic and gravity modelling routines.
GeoData	Data layer class.

Functions

calc_demag(mvec, k, dxy, dz)	Calculate demagnetisation correction.
save_layer(mlist)	Routine to save the mlayer and glayer to a file.
gridmatch(lmod, ctxt, rtxt)	Match the rows and columns of the second grid to the first grid.
calc_field(lmod[, pbars, showtext, parent, ...])	Calculate magnetic and gravity field.
sum_fields(k, mgval, numx, numy, modind, aaa0, aaa1, ...)	Sum magnetic and gravity field datasets to produce final model field.
quick_model([numx, numy, numz, dxy, d_z, tlx, tly, ...])	Quick model function.
dircos(incl, decl, azim)	Compute direction cosines from inclination and declination.
dat_extent(dat, axes)	Get the extent of the dat variable.

Module Contents

class pygmi.pfmod.grvmag3d.PlotTest(data=None)

Bases: pygmi.misc.ContextModule

Plot Raster Class.

update_raster(lmod2)

Update the raster plot.

Parameters:

lmod2 (pygmi.pfmod.datatypes.LithModel) – Lithology dataset to be used.

class pygmi.pfmod.grvmag3d.GravMag(parent=None)

The GravMag class holds generic magnetic and gravity modelling routines.

Routine that will calculate the final versions of the field. Other, related code is here as well, such as the inversion routines.

Parameters:

parent – Reference to the parent routine. The default is None.

setupui()

Set up UI.

calc_field_mag()

Pre field-calculation routine.

calc_field_grav()

Pre field-calculation routine.

calc_field_mag_changes()

Calculate only magnetic field changes.

calc_field_grav_changes()

Calculate only gravity field changes.

calc_field2(showreports=False, magcalc=False)

Calculate magnetic and gravity field.

Parameters:

• showreports (bool) – Flag for showing reports. The default is False.

• magcalc (bool) – Flag for choosing the magnetic calculation. The default is False.

calc_regional()

Calculate magnetic and gravity regional.

Calculates a gravity and magnetic regional value based on a single solid lithology model. The principle is that the maximum value for a solid model with fixed extents and depth, using the most COMMON lithology, would be the MAXIMUM AVERAGE value for any model which we would do. Therefore the regional is simply:

REGIONAL = OBS GRAVITY MEAN - CALC GRAVITY MAX

This routine calculates the last term.

test_pattern()

Displays a test pattern of the data.

This is an indication of the edge of model field decay. It gives an idea about how reliable the calculated field on the edge of the model is.

update_graph(grvval, magval, modind)

Update the graph.

Parameters:

• grvval (numpy.typing.NDArray) – Array of gravity values.

• magval (numpy.typing.NDArray) – Array of magnetic values.

• modind (numpy.typing.NDArray) – Model indices.

class pygmi.pfmod.grvmag3d.GeoData(parent, ncols=10, nrows=10, numz=10, dxy=10.0, d_z=10.0, mht=80.0, ght=0.0)

Data layer class.

This class defines each geological type and calculates the field for one cube from the standard definitions.

The is a class which contains the geophysical information for a single lithology. This includes the final calculated field for that lithology only.

Parameters:

• parent – Reference to the parent routine.

• ncols (int) – Number of columns in the model.

• nrows (int) – Number of rows in the model.

• numz (int) – Number of layer in the model.

• dxy (int) – X and Y size of each voxel.

• d_z (float) – Layer thickness.

• mht (float) – Magnetic sensor height.

• ght (float) – Gravity sensor height.

calc_origin_grav(hcor=None)

Calculate the field values for the lithologies.

Parameters:

hcor (numpy.typing.NDArray) – Height corrections. The default is None.

calc_origin_mag(hcor=None, demag=False)

Calculate the field values for the lithologies.

Parameters:

• hcor (numpy.typing.NDArray) – Height corrections. The default is None.

• demag (bool) – Do demagnetization.

rho()

Return the density contrast.

Returns:

Density contrast.

Return type:

float

set_xyz(ncols, nrows, numz, g_dxy, mht, ght, d_z, dxy=None, modified=True)

Sets/updates xyz parameters.

Parameters:

• ncols (int) – Number of columns.

• nrows (int) – Number of rows.

• numz (int) – Number of layers.

• g_dxy (float) – Grid spacing in x and y direction.

• mht (float) – Magnetic sensor height.

• ght (float) – Gravity sensor height.

• d_z (float) – Model spacing in z direction.

• dxy (float, optional) – Model spacing in x and y direction. The default is None.

• modified (bool, optional) – Whether the model was modified. The default is True.

set_xyz12()

Set x12, y12, z12.

This is the limits of the cubes for the model

gboxmain(xobs, yobs, zobs, hcor)

Gbox routine by Blakely.

Note: xobs, yobs and zobs must be floats or there will be problems later.

Subroutine GBOX computes the vertical attraction of a rectangular prism. Sides of prism are parallel to x,y,z axes, and z axis is vertical down.

Input parameters: | Observation point is (x0,y0,z0). The prism extends from x1 | to x2, from y1 to y2, and from z1 to z2 in the x, y, and z | directions, respectively. Density of prism is rho. All | distance parameters in units of m;

Output parameters: | Vertical attraction of gravity, g, in mGal/rho. | Must still be multiplied by rho outside routine. | Done this way for speed.

Parameters:

• xobs (numpy.typing.NDArray) – Observation X coordinates.

• yobs (numpy.typing.NDArray) – Observation Y coordinates.

• zobs (float) – Observation Z coordinates.

• hcor (numpy.typing.NDArray) – Height corrections.

mboxmain(xobs, yobs, zobs, hcor, demag=False)

Mbox routine by Blakely.

Note: xobs, yobs and zobs must be floats or there will be problems later.

Subroutine MBOX computes the total field anomaly of an infinitely extended rectangular prism. Sides of prism are parallel to x,y,z axes, and z is vertical down. Bottom of prism extends to infinity. Two calls to mbox can provide the anomaly of a prism with finite thickness; e.g.,

call mbox(x0,y0,z0,x1,y1,z1,x2,y2,mi,md,fi,fd,m,theta,t1)

call mbox(x0,y0,z0,x1,y1,z2,x2,y2,mi,md,fi,fd,m,theta,t2)

t=t1-t2

Requires subroutine DIRCOS. Method from Bhattacharyya (1964).

Input parameters: | Observation point is (x0,y0,z0). Prism extends from x1 to | x2, y1 to y2, and z1 to infinity in x, y, and z directions, | respectively. Magnetization defined by inclination mi, | declination md, intensity m. Ambient field defined by | inclination fi and declination fd. X axis has declination | theta. Distance units are irrelevant but must be consistent. | Angles are in degrees, with inclinations positive below | horizontal and declinations positive east of true north. | Magnetization in A/m.

Output parameters: | Total field anomaly t, in nT.

Parameters:

• xobs (numpy.typing.NDArray) – Observation X coordinates.

• yobs (numpy.typing.NDArray) – Observation Y coordinates.

• zobs (float) – Observation coordinate.

• hcor (numpy.typing.NDArray) – Height corrections.

• demag (bool) – Do demagnetization.

pygmi.pfmod.grvmag3d.calc_demag(mvec, k, dxy, dz)

Calculate demagnetisation correction.

Parameters:

• mvec (numpy.typing.NDArray) – Body Magnetisation.

• k (float) – susceptibility.

• dxy (float) – cell width.

• dz (float) – cell height.

Returns:

Corrected magnetisation.

Return type:

ndarray

pygmi.pfmod.grvmag3d.save_layer(mlist)

Routine to save the mlayer and glayer to a file.

Parameters:

mlist (tuple[str, pygmi.pfmod.datatypes.LithModel]) – Tuple with 2 elements - lithology name and LithModel.

Returns:

Link to a temporary file.

Return type:

IO

pygmi.pfmod.grvmag3d.gridmatch(lmod, ctxt, rtxt)

Match the rows and columns of the second grid to the first grid.

Parameters:

• lmod (pygmi.pfmod.datatypes.LithModel) – Lithology Model.

• ctxt (str) – First grid text label.

• rtxt (str) – Second grid text label.

Returns:

Numpy array of data.

Return type:

ndarray

pygmi.pfmod.grvmag3d.calc_field(lmod, pbars=None, showtext=None, parent=None, showreports=False, magcalc=False, demag=False)

Calculate magnetic and gravity field.

This function calculates the magnetic and gravity field. It has two different modes of operation, by using the magcalc switch. If magcalc=True then magnetic fields are calculated, otherwise only gravity is calculated.

Parameters:

• lmod (pygmi.pfmod.datatypes.LithModel) – PyGMI lithological model

• pbars – progress bar routine if available. (internal use)

• showtext – showtext routine if available. (internal use)

• parent – parent routine, the default is None.

• showreports (bool) – show extra reports, the default is False.

• magcalc (bool) – if True, calculates magnetic data, otherwise only gravity.

• demag (bool) – calculate demagnetization. The default is False.

Returns:

dictionary of items of type Data.

Return type:

dict of Data

pygmi.pfmod.grvmag3d.sum_fields(k, mgval, numx, numy, modind, aaa0, aaa1, mlayers, hcorflat, mijk)

Sum magnetic and gravity field datasets to produce final model field.

Parameters:

• k (int) – k index.

• mgval (numpy.typing.NDArray) – Magnetic or gravity data being summed.

• numx (int) – Number of x elements.

• numy (int) – Number of y elements.

• modind (numpy.typing.NDArray) – model with indices representing lithologies.

• aaa0 (numpy.typing.NDArray) – x indices for offsets.

• aaa1 (numpy.typing.NDArray) – y indices for offsets.

• mlayers (numpy.typing.NDArray) – Layer fields for summation.

• hcorflat (numpy.typing.NDArray) – Height correction.

• mijk (int) – Current lithology index.

Returns:

Output summed data.

Return type:

ndarray

pygmi.pfmod.grvmag3d.quick_model(numx=50, numy=40, numz=5, dxy=100.0, d_z=100.0, tlx=0.0, tly=0.0, tlz=0.0, mht=100.0, ght=0.0, finc=-67, fdec=-17, inputliths=None, susc=None, dens=None, minc=None, mdec=None, mstrength=None, hintn=30000.0)

Quick model function.

Parameters:

• numx (int) – Number of x elements. The default is 50.

• numy (int) – Number of y elements. The default is 40.

• numz (int) – number of z elements (layers). The default is 5.

• dxy (float) – Cell size in x and y direction. The default is 100..

• d_z (float) – Layer thickness. The default is 100..

• tlx (float) – Top left x coordinate. The default is 0..

• tly (float) – Top left y coordinate. The default is 0..

• tlz (float) – Top left z coordinate. The default is 0..

• mht (float) – Magnetic sensor height. The default is 100..

• ght (float) – Gravity sensor height. The default is 0..

• finc (float) – Magnetic field inclination (degrees). The default is -67.

• fdec (float) – Magnetic field declination (degrees). The default is -17.

• inputliths (list[str] | None) – List of input lithologies. The default is None.

• susc (list[float] | None) – List of susceptibilities. The default is None.

• dens (list[float] | None) – List of densities. The default is None.

• minc (list[float] | None) – List of remanent inclinations (degrees). The default is None.

• mdec (list[float] | None) – List of remanent declinations (degrees). The default is None.

• mstrength (list[float] | None) – List of remanent magnetisations (A/m). The default is None.

• hintn (float) – Magnetic field strength (nT). The default is 30000.

Returns:

Output model.

Return type:

LithModel

pygmi.pfmod.grvmag3d.dircos(incl, decl, azim)

Compute direction cosines from inclination and declination.

Parameters:

• incl (float) – inclination in degrees positive below horizontal.

• decl (float) – declination in degrees positive east of true north.

• azim (float) – azimuth of x axis in degrees positive east of north.

Returns:

• aaa (float) – First direction cosine.

• bbb (float) – Second direction cosine.

• ccc (float) – Third direction cosine.

Return type:

tuple[float, float, float]

pygmi.pfmod.grvmag3d.dat_extent(dat, axes)

Get the extent of the dat variable.

Parameters:

• dat (pygmi.raster.datatypes.Data) – PyGMI raster dataset.

• axes (matplotlib.axes.Axes) – Matplotlib axes.

Returns:

• left – Left coordinate.

• right – Right coordinate.

• bottom – Bottom coordinate.

• top – Top coordinate.

Return type:

tuple[float, float, float, float]