Welcome to pydEXP’s documentation!

pyDEXP is a open-source python package aiming at processing potential field data using the dEXP theory formulated by Fedi et al., 2012. The package largely benefits from the imaging methods module of the Fatiando a terra open-source python library.

References

  • Uieda, L., V. C. Oliveira Jr, and V. C. F. Barbosa (2013), Modeling the Earth with Fatiando a Terra, Proceedings of the 12th Python in Science Conference, pp. 91 - 98.

  • Fedi, M., and M. Pilkington (2012), Understanding imaging methods for potential field data, Geophysics, 77(1), G13, doi:10.1190/geo2011-0078.1

Quick use

Clone the gitlab repository:

git clone https://github.com/BenjMy/dEXP_imaging.git

From the same dEXP_imaging directory you can import the module from source using python.

Contribute and support

Citing pydEXP

If you use the pydEXP code for you work, please cite this paper as:

Paper in preparation

Getting Started

pyDEXP aims to process Mise-à-la-masse (MALM) datasets for a variety of applications. pyDEXP has been initially developed for plant root imaging.

The simpliest processing can be achieved with the python API. You’ll first need to import the pyDEXP package:

import lib.dEXP as dEXP
import lib.plot_dEXP as pEXP

Then after loading you data make basics operation such as upward continuation (Uieda et al., 2013, 2018):

mesh, label_prop = dEXP.upwc(xp, yp, zp, U, shape,
                 zmin=0, zmax=max_elevation, nlayers=nlay,
                 qorder=qorder)

Searching for ridges (Fedi et al., 2012):

dfI,dfII, dfIII = dEXP.ridges_minmax(xp, yp, mesh, p1, p2,
                                      label=label_prop,
                                      fix_peak_nb=2,

and finally plot the results:

fig = plt.figure()
ax = plt.gca()
pEXP.plot_xy(mesh, label=label_prop, ax=ax)
pEXP.plot_ridges_harmonic(dfI,dfII,dfIII,ax=ax)

More examples are available in the Example gallery section.

References

Uieda, L., V. C. Oliveira Jr, and V. C. F. Barbosa (2013), Modeling the Earth with Fatiando a Terra, Proceedings of the 12th Python in Science Conference, pp. 91 - 98.

Uieda, L. (2018). Verde: Processing and gridding spatial data using Green’s functions. Journal of Open Source Software, 3(29), 957. doi:10.21105/joss.00957

Fedi, M., and M. Pilkington (2012), Understanding imaging methods for potential field data, Geophysics, 77(1), G13, doi:10.1190/geo2011-0078.1

Which Python?

For the moment, pydEXP is only tested on Python 3.7. First, make sure you have all dependencies (see below) installed.

Clone the gitlab repository:

git clone https://github.com/BenjMy/dEXP_imaging.git

From the same dEXP_imaging directory you can import the module from source using python.

Dependencies

pydEXP requires the following dependencies for running:

  • numpy

  • scipy

  • matplotlib

  • pandas

  • mpl_axes_aligner

In order to make upward continuation and derivate the field, pyDEXP uses: * fatiando

Optionnal Third party packages

In order to forward model the geolectrical data, pyDEXP uses:

  • Pygimli

  • Resipy

You can test the installation running one of the exemple.

API documentation

API reference: The pydEXP package

Imaging methods for potential fields.

Implements the DEXP method described in Fedi and Pilkington (2012).

Note

This is part of a larger project aiming at inverting current sources density (see more at: https://icsd-dev.readthedocs.io/en/latest/)

References

Fedi, M., and M. Pilkington (2012), Understanding imaging methods for potential field data, Geophysics, 77(1), G13, doi:10.1190/geo2011-0078.1

dEXP.dEXP(x, y, z, data, shape, zmin, zmax, nlayers, qorder=0, SI=1)[source]

DEXP model (Fedi, 2012). (NOT YET TESTED)

Calculates a physical property distribution given potential field data on a regular grid. Uses depth weights. Parameters:

  • x, y1D-arrays

    The x and y coordinates of the grid points

  • zfloat or 1D-array

    The z coordinate of the grid points

  • data1D-array

    The potential field at the grid points

  • shapetuple = (ny, nx)

    The shape of the grid

  • zmin, zmaxfloat

    The top and bottom, respectively, of the region where the physical property distribution is calculated

  • nlayersint

    The number of layers used to divide the region where the physical property distribution is calculated

  • qorderfloat

    The order of the vertical derivative

  • SIfloat

    The structural index

Returns:

  • meshfatiando.mesher.PrismMesh

    The estimated physical property distribution set in a prism mesh (for easy 3D plotting)

dEXP.dEXP_ratio(x, y, z, data, shape, zmin, zmax, nlayers, qorders=[1, 0])[source]

DEXP ratio model (NOT YET validated) Abbas, M. A., and Fedi, M. (2014). Automatic DEXP imaging of potential fields independent of the structural index. Geophysical Journal International, 199 (3), 1625-1632.

Calculates a physical property distribution given potential field data on a regular grid. Uses depth weights. Parameters:

  • x, y1D-arrays

    The x and y coordinates of the grid points

  • zfloat or 1D-array

    The z coordinate of the grid points

  • data1D-array

    The potential field at the grid points

  • shapetuple = (ny, nx)

    The shape of the grid

  • zmin, zmaxfloat

    The top and bottom, respectively, of the region where the physical property distribution is calculated

  • nlayersint

    The number of layers used to divide the region where the physical property distribution is calculated

  • qorders1D-array

    The order of the derivatives for the ratio calculation

Returns:

  • meshfatiando.mesher.PrismMesh

    The estimated physical property distribution set in a prism mesh (for easy 3D plotting)

dEXP.filter_ridges(dfIf, dfIIf, dfIIIf, minDepth, maxDepth, minlength=3, rmvNaN=False, **kwargs)[source]

Filter non-rectiligne ridges (denoising)

Parameters: dfI: dataframe

contains ridges of type I

dfII: dataframe

contains ridges of type II

dfIII: dataframe

contains ridges of type III

  • minDepth

    Text here

  • maxDepth

  • kwargs

heights: dataframe

contains heights of ridges of type I

Returns:

  • BB :

    Text here

dEXP.fit_ridges(df, rmvOutliers=False)[source]

Fit ridges and return points and fit equations to plot

Parameters:

  • df

    dataframe including all tree types of ridges

Returns:

  • BB :

    points and fit equations to plot

dEXP.peakdet(v, delta, x=None)[source]

Converted from MATLAB script at http://billauer.co.il/peakdet.html

Returns two arrays

function [maxtab, mintab]=peakdet(v, delta, x) %PEAKDET Detect peaks in a vector % [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local % maxima and minima (“peaks”) in the vector V. % MAXTAB and MINTAB consists of two columns. Column 1 % contains indices in V, and column 2 the found values. % % With [MAXTAB, MINTAB] = PEAKDET(V, DELTA, X) the indices % in MAXTAB and MINTAB are replaced with the corresponding % X-values. % % A point is considered a maximum peak if it has the maximal % value, and was preceded (to the left) by a value lower by % DELTA.

% Eli Billauer, 3.4.05 (Explicitly not copyrighted). % This function is released to the public domain; Any use is allowed.

dEXP.profile_extra(x, y, v, point1, point2, size, algorithm='cubic')[source]

Extract a profile between 2 points from spacial data.

Uses interpolation to calculate the data values at the profile points.

Parameters:

  • x, y1D arrays

    Arrays with the x and y coordinates of the data points.

  • v1D array

    Array with the scalar value assigned to the data points.

  • point1, point2lists = [x, y]

    Lists the x, y coordinates of the 2 points between which the profile will be extracted.

  • sizeint

    Number of points along the profile.

  • algorithmstring

    Interpolation algorithm. Either 'cubic', 'nearest', 'linear' (see scipy.interpolate.griddata).

Returns:

  • [xp, yp, distances, vp]1d arrays

    xp and yp are the x, y coordinates of the points along the profile. distances are the distances of the profile points from point1. vp are the data points along the profile.

dEXP.profile_noInter(x, y, v, point1, point2, size=None, **kwargs)[source]

Extract a profile between 2 points from spacial data.

NO interpolation to calculate the data values at the profile points (find nearest point).

Parameters:

  • x, y1D arrays

    Arrays with the x and y coordinates of the data points.

  • v1D array

    Array with the scalar value assigned to the data points.

  • point1, point2lists = [x, y]

    Lists the x, y coordinates of the 2 points between which the profile will be extracted.

  • sizeint

    Number of points along the profile.

Returns:

  • [xp, yp, distances, vp]1d arrays

    xp and yp are the x, y coordinates of the points along the profile. distances are the distances of the profile points from point1. vp are the data points along the profile.

dEXP.ridges_intersection_Z0(fit, ax=None, ridge_nb=None)[source]

Find intersection of ridges (NOT YET IMPLEMENTED)

Parameters:

  • a

    Text here

Returns:

  • BB :

    return intersection by pairs

dEXP.ridges_minmax(x, y, mesh, p1, p2, qorder=0, z=0, label='upwc', fix_peak_nb=None, interp=True, smooth=False, showfig=False, **kwargs)[source]

Form a multiridge set RI and RII : zeros of the first horizontal and first vertical derivatives of the potential field RIII :zeros of the potential field itself

Note

ridges generated by isolated simple sources point, line, sheet, and contact are straight lines defined by the zeros of a potential

field and its horizontal and vertical derivatives at all measured or computed levels

Parameters:

  • x, y1D-arrays

    The x and y coordinates of the grid points

  • meshfatiando.mesher.PrismMesh

    The estimated physical property distribution set in a prism mesh (for easy 3D plotting). The upward continuated field mesh (of order-q derivative). Read the property label ‘upwc’ of the mesh by default

  • p1, p21D-arrays

    The p1 and p2 coordinates of the extracted profile end points

  • qorderint

    The derivative order

  • zfloat or 1D-array

    The z coordinate of the grid points

  • labelstring

    Label of the estimated physical property distribution

  • fix_peak_nbint

    The maximum number of peak to identify

  • interpTrue or False

    If True, will interpolate values between the data points.

  • smoothTrue or False

    If True, will apply a low-pass filter values.

  • showfigTrue or False

    If True, will display the figure.

**kwargs

prominence for peak detection reverse: start peak analysis from bottom to top

Returns:

  • MinMax_peaks :

    Panda dataframe containing ridges RI, RII and RII

dEXP.ridges_minmax_plot(x, y, mesh, p1, p2, qorder=0, z=0, fix_peak_nb=None, label='upwc', interp=True, smooth=False, showfig=False, **kwargs)[source]
  • x, y1D-arrays

    The x and y coordinates of the grid points

  • meshfatiando.mesher.PrismMesh

    The estimated physical property distribution set in a prism mesh (for easy 3D plotting)

  • p1, p21D-arrays

    The p1 and p2 coordinates of the extracted profile end points

  • qorderint

    The derivative order

  • zfloat or 1D-array

    The z coordinate of the grid points

  • fix_peak_nbint

    The maximum number of peak to identify

  • labelstring

    Label of the estimated physical property distribution

  • interpTrue or False

    If True, will interpolate values between the data points.

  • smoothTrue or False

    If True, will apply a low-pass filter values.

  • showfigTrue or False

    If True, will display the figure.

Returns:

dEXP.scalEULER(df, EXTnb=[1], z0=0)[source]

Analysis (Euler deconvolution of ridges, Fedi 2019) (NOT YET IMPLEMENTED)

Parameters:

  • a

    Text here

Returns:

  • BB :

    Text here

dEXP.scalFUN(df, EXTnb=[1], z0=0, **kwargs)[source]

Analysis of ridges (NOT YET IMPLEMENTED)

Parameters:

  • a

    Text here

Returns:

  • BB :

    Text here

dEXP.upwc(x, y, z, data, shape, zmin, zmax, nlayers, qorder=0)[source]

Upward continuation model (Fedi, 2012).

Calculates the upward continuation for given potential field data on a regular grid. Parameters:

  • x, y1D-arrays

    The x and y coordinates of the grid points

  • zfloat or 1D-array

    The z coordinate of the grid points

  • data1D-array

    The potential field at the grid points

  • shapetuple = (ny, nx)

    The shape of the grid

  • zmin, zmaxfloat

    The top and bottom, respectively, of the region where the physical property distribution is calculated

  • nlayersint

    The number of layers used to divide the region where the physical property distribution is calculated

  • qorderfloat

    The order of the vertical derivative

Returns:

  • meshfatiando.mesher.PrismMesh

    The estimated physical property distribution set in a prism mesh (for easy 3D plotting)

plot_dEXP: plotter

Some functions to plot results from DEXP transformation

plot_dEXP.label2unit(label)[source]

Convert label name to unit

Parameters

label : str

Returns

unit : str

plot_dEXP.plot_line(x, y, data, p1, p2, ax=None, interp=True, **kwargs)[source]

Parameters

xTYPE

x coords of the 2d map.

yTYPE

y coords of the 2d map..

dataTYPE

2d map dataset.

p1list

initial point coordinate (x,y) of the extracted profile.

p2list

final point coordinate (x,y) of the extracted profile.

**kwargs

Returns

xxnp.array

interpolated x.

yynp.array

interpolated y.

distancenp.array

profile distance.

profileTYPE

data values along the profile.

plot_dEXP.plot_ridges_harmonic(RI=None, RII=None, RIII=None, ax=None, label=False, legend=True, **kwargs)[source]

Plot ridges in the harmonic domain

Parameters:

  • a

    Text here

Returns:

  • BB

    Text here

plot_dEXP.plot_ridges_sources(df_fit, ax=None, ridge_type=None, ridge_nb=None, z_max_source=None, **kwargs)[source]

Plot ridges in the source domain and observe intersection point

Parameters: * df_fit

dataframe fit

Returns:

  • BB

    Text here

plot_dEXP.plot_scalFUN(points, fit, ax=None, z0=None, label=None)[source]

Plot scalfun function analysis

Parameters:

  • a

    Text here

Returns:

  • BB

    Text here

plot_dEXP.plot_xy(mesh, scaled=0, label=None, ax=None, markerMax=False, **kwargs)[source]

Get vertical xy section of the mesh at max/2 and plot it

Parameters:

  • mesh

    Fatiando mesh type

  • scaled

    Txt here

  • label

    Txt here

  • ax

    Txt here

  • markerMax

    Txt here

plot_dEXP.plot_z(mesh)[source]

Get horizontal sections at z levels and plot it

Parameters:

  • mesh

    Fatiando mesh type

plot_dEXP.slice_mesh(x, y, mesh, label, p1, p2, ax=None, interp=True, **kwargs)[source]

Get vertical xy section of the mesh at max/2 and plot it

Parameters:

  • mesh

    Fatiando mesh type

  • scaled

    Txt here

  • label

    Txt here Txt here

  • markerMax

    Txt here

utils_dexp: utils

Miscellaneous utility functions.

@author: Benjamin

utils_dEXP.cor_field_B(x, y, z, u, B, rho=100, **kwargs)[source]

Calculates the potential field (electric) produced by a current injection in B (return electrode) for a given homogeneous electrical resistivity rho.

Parameters

  • x, y1D-arrays

    The x and y coordinates of the grid points

  • zfloat or 1D-array

    The z coordinate of the grid points

  • u1D-array

    The potential field at the grid points

  • Bu1D-array

    The position of the return electrode B

  • rhoint

    The estimated electrical resistivity of the medium

Returns

u_cor1D-arrays

The corrected from the B return electrode influence potential field at the grid points

utils_dEXP.load_surfer(fname, dtype='float64')[source]

Read data from a Surfer ASCII grid file.

Surfer is a contouring, griding and surface mapping software from GoldenSoftware. The names and logos for Surfer and Golden Software are registered trademarks of Golden Software, Inc.

http://www.goldensoftware.com/products/surfer

Parameters:

  • fnamestr

    Name of the Surfer grid file

  • dtypenumpy dtype object or string

    The type of variable used for the data. Default is numpy.float64. Use numpy.float32 if the data are large and precision is not an issue.

Returns:

  • datadict

    The data in a dictionary with some metadata:

    • 'file'string

      The name of the original data file

    • 'shape'tuple

      (nx, ny), the number of grid points in x (North) and y (East)

    • 'area'tuple

      (x1, x2, y1, y2), the spacial limits of the grid.

    • 'x'1d-array

      Value of the North-South coordinate of each grid point.

    • 'y'1d-array

      Value of the East-West coordinate of each grid point.

    • 'data'1d-array

      Values of the field in each grid point. Field can be for example topography, gravity anomaly, etc. If any field values are >= 1.70141e+38 (Surfers way of marking NaN values), the array will be masked at those values (i.e., 'data' will be a numpy masked array).

Indices and tables