API reference¶

The public API of PyCurious. Everything below is importable directly from the top-level pycurious namespace.

Grid and spectra¶

The base class shared by both optimisers — subgrid decomposition, the radial and window spectra, and the covariance machinery behind the uncertainties.

class pycurious.CurieGrid(grid, xmin, xmax, ymin, ymax, **kwargs)[source]¶

Bases: CurieParallel

Accepts a 2D array and Cartesian coordinates specifying the bounding box of the array

Grid must be projected in metres.

Parameters:

grid – 2D numpy array 2D array of magnetic data
xmin – float minimum x bound in metres
xmax – float maximum x bound in metres
ymin – float minimum y bound in metres
ymax – float maximum y bound in metres

grid¶: 2D numpy array 2D array of magnetic data

xmin¶: float minimum x bound in metres

xmax¶: float maximum x bound in metres

ymin¶: float minimum y bound in metres

ymax¶: float maximum y bound in metres

dx¶: float grid spacing in the x-direction in metres

dy¶: float grid spacing in the y-direction in metres

nx¶: int number of nodes in the x-direction

ny¶: int number of nodes in the y-direction

xcoords¶: 1D numpy array 1D numpy array of coordinates in the x-direction

ycoords¶: 1D numpy array 1D numpy array of coordinates in the y-direction

Notes

In all instances x indicates eastings in metres and y indicates northings. Using a grid of longitude / latitudinal coordinates (degrees) will result in incorrect Curie depth calculations.

subgrid(window, xc, yc)[source]¶

Extract a subgrid from the data at a window around the point (xc,yc)

Parameters:

xc – float x coordinate
yc – float y coordinate
window – float size of window in metres

Returns:

data –

2D array: subgrid encompassing window size

create_centroid_list(window, spacingX=None, spacingY=None)[source]¶

Create a list of xc,yc values to extract subgrids.

Parameters:

window – float size of the windows in metres
spacingX – float (optional) specify spacing in metres in the X direction will default to maximum X resolution
spacingY – float (optional) specify spacing in metres in the Y direction will default to maximum Y resolution

Returns:

xc_list –

1D array: array of x coordinates
yc_list1D array: array of y coordinates

remove_trend_linear(data)[source]¶

Remove the best-fitting linear trend from the data

This may come in handy if the magnetic data has not been reduced to the pole.

The trend is the least-squares plane. Over a regular grid the centred row and column indices are mutually orthogonal and both orthogonal to the constant, so the normal equations decouple and the plane’s three coefficients are one mean and two 1-D inner products – no design matrix and no SVD. This is an order of magnitude cheaper than the equivalent lstsq fit (the trend is subtracted once per window when computing a spectrum), and unlike an (nr, nc) vs (nc, nr) design matrix it stays correct when the grid is not square.

Parameters:: data – 2D numpy array
Returns:: data – 2D numpy array

radial_spectrum(subgrid, taper=<function hanning>, power=2.0, return_counts=False, **kwargs)[source]¶

Compute the radial spectrum for a square grid.

Wavenumber is returned in units of rad/km.

Parameters:

subgrid – 2D array window of the original data (see subgrid method)
taper – function (default=np.hanning) taper function, set to None for no taper function
power –
float raise the FFT of the magnetic anomaly to the power:
- 2.0 for Bouligand et al. (2009) use cases, which gives the log power spectrum \(\ln \Phi_{\Delta T}\)
- 1.0 for Tanaka et al. (1999) use cases, which gives the log amplitude spectrum \(\ln \Phi_{\Delta T}^{1/2}\)
return_counts – bool (default=False) also return the number of FFT cells averaged into each bin
kwargs – keyword arguments keyword arguments to pass to taper

Returns:

k –

1D array shape (n,): wavenumber in rad/km
Phi1D array shape (n,): Radial power spectrum
sigma_Phi1D array shape (n,): Standard deviation of Phi within each radial bin
counts1D array shape (n,): number of FFT cells in each radial bin. Only returned if return_counts=True.

Notes

Phi is the mean of \(\ln |FFT|\) over each annulus, so sigma_Phi describes the scatter of the individual cells, not the uncertainty of that mean. Dividing by the square root of counts gives the standard error, though note the cells are not independent – a real field has Hermitian symmetry, so roughly half of them are redundant, and tapering correlates neighbours.

While subgrid is projected in eastings / northings (in metres), the wavenumber, \(k\), is returned in units of rad/km. This is because both Bouligand et al. (2009) and Tanaka et al. (1999) require the computation of Curie depth in these units.

References

Bouligand, C., J. M. G. Glen, and R. J. Blakely (2009), Mapping Curie temperature depth in the western United States with a fractal model for crustal magnetization, J. Geophys. Res., 114, B11104, doi:10.1029/2009JB006494

Tanaka, A., Okubo, Y., & Matsubayashi, O. (1999). Curie point depth based on spectrum analysis of the magnetic anomaly data in East and Southeast Asia. Tectonophysics, 306(3–4), 461–470. doi:10.1016/S0040-1951(99)00072-4

window_spectrum(window, xc, yc, taper=<function hanning>, power=2.0, process_subgrid=None, dof_factor=None, **kwargs)[source]¶

Radial spectrum of one window, weighted ready for fitting.

Extracts the subgrid, computes its radial spectrum, and converts the within-annulus scatter into the uncertainty of the annulus mean, which is what a fit needs. Both pycurious.optimise_bouligand and pycurious.optimise_tanaka build on this.

Parameters:

window – float size of the window in metres
xc – float centroid of the window
yc – float centroid of the window
taper – function (default=np.hanning) taper function, or None for no taper
power – float raise the FFT of the anomaly to this power – 2.0 for the log power spectrum that Bouligand et al. (2009) fit, 1.0 for the log amplitude spectrum of Tanaka et al. (1999)
process_subgrid – function, optional applied to the subgrid before the spectrum is computed
dof_factor – float, optional override the effective-degrees-of-freedom deflation (see Notes)
kwargs – keyword arguments passed to taper

Returns:

k –

1D array: wavenumber in rad/km
Phi1D array: log spectrum, raised to power
sigma1D array: uncertainty of the binned mean

Usage:

>>> k, Phi, sigma = grid.window_spectrum(200e3, xc, yc, power=2)

Notes

radial_spectrum returns the scatter of the FFT cells within each annulus, whereas a fit needs the uncertainty of the annulus mean. That is the standard error, except that the cells are not independent: Hermitian symmetry makes about half of them redundant, and tapering correlates neighbours. The correction is calibrated per taper and varies with the number of cells in the bin, since a fixed number of them is lost to correlation however few there are. dof_factor overrides it with a constant.

The cells of neighbouring annuli are correlated too, which this does not address – it inflates the uncertainty of a fitted parameter rather than of any individual bin. See pycurious.optimise_bouligand.CurieOptimiseBouligand.optimise.

reduce_to_pole(data, inc, dec, sinc=None, sdec=None)[source]¶

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 (Blakely, 1996).

Parameters:

data – 1D array the total field anomaly data at each point.
inc – float / 1D array inclination of the inducing Geomagnetic field
dec – float / 1D array declination of the inducing Geomagnetic field
sinc – float / 1D array (optional) inclination of the total magnetization of the anomaly source
sdec – float / 1D array (optional) 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 –

2D array: the data reduced to the pole.

References

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

Notes

This functions performs the reduction in the frequency domain (using the FFT). The transform filter is (in the freq 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, k_y\) are the wave-numbers in the x and y directions and

\(|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\)

\(\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.

upward_continuation(data, height)[source]¶

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, \(|k|\) is the wavenumber modulus.

Parameters:

data – 2D array potential field at the grid points
height – float height increase (delta z) in meters.

Returns:

cont –

array: upward continued data

References

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

parallelise_routine(window, xc_list, yc_list, func, *args, **kwargs)¶

Implements shared memory multiprocessing to split multiple evaluations of a function centroids across processors.

Supply the window size and lists of x,y coordinates to a function along with any additional arguments or keyword arguments.

Parameters:

window – float size of window in metres
xc_list – array shape (l,) centroid x values
yc_list – array shape (l,) centroid y values
func – function Python function to evaluate in parallel
args – arguments additional arguments to pass to func
kwargs –
keyword arguments additional keyword arguments to pass to func. Two keys are reserved and consumed here rather than forwarded:
- on_error : {“raise”, “ignore”} (default=”raise”) what to do when func fails for a centroid. "ignore" fills that centroid with NaN and warns.
- seed : int or None (default=None) seeds any stochastic routine reproducibly. Each centroid receives an independent child seed derived from numpy.random.SeedSequence, passed to func as a seed keyword, so results do not depend on how many processors were used. func must accept a seed keyword if this is supplied.

Returns:

out –

list of lists: (depends on output of func - see notes)

Usage:

An obvious use case is to compute the Curie depth for many centroids in parallel.

>>> self.parallelise_routine(window, xc_list, yc_list, self.optimise)

Each centroid is assigned a new process and sent to a free processor to compute. In this case, the output is separate lists of shape(l,) for \(\beta, z_t, \Delta z, C\). If len(xc_list)=2 then,

>>> self.parallelise_routine(window, [x1,x2], [y1, y2], self.optimise)
[[beta1  beta2], [zt1  zt2], [dz1  dz2], [C1  C2]]

Another example is to parallelise the sensitivity analysis:

>>> self.parallelise_routine(window, xc_list, yc_list, self.sensitivity, nsim)

This time the output will be a list of lists for \(\beta, z_t, \Delta z, C\) i.e. if len(xc_list)=2 is the number of centroids and nsim=4 is the number of simulations then separate lists will be returned for \(\beta, z_t, \Delta z, C\).

>>> self.parallelise_routine(window, [x1,x2], [y1,y2], self.sensitivity, 4)

which would return:

[[[ beta1a , beta1b , beta1c , beta1d ],   # centroid 1 (x1,y1)
  [ beta2a , beta2b , beta2c , beta2d ]],  # centroid 2 (x2,y2)
 [[   zt1a ,   zt1b ,   zt1c ,   zt1d ],   # centroid 1 (x1,y1)
  [   zt2a ,   zt2b ,   zt2c ,   zt2d ]],  # centroid 2 (x2,y2)
 [[   dz1a ,   dz1b ,   dz1c ,   dz1d ],   # centroid 1 (x1,y1)
  [   dz2a ,   dz2b ,   dz2c ,   dz2d ]]   # centroid 2 (x2,y2)
 [[    C1a ,    C1b ,    C1c ,    C1d ],   # centroid 1 (x1,y1)
  [    C2a ,    C2b ,    C2c ,    C2d ]]]  # centroid 2 (x2,y2)

Notes

See the module docstring for the if __name__ == "__main__": requirement when calling this from a script.

Bouligand optimiser¶

class pycurious.CurieOptimiseBouligand(grid, xmin, xmax, ymin, ymax, **kwargs)[source]¶

Bases: CurieGrid

Extends the pycurious.grid.CurieGrid class to include optimisation routines see scipy.optimize.minimize for a description of the algorithm.

Parameters:

grid – 2D numpy array 2D array of magnetic data
xmin – float minimum x bound in metres
xmax – float maximum x bound in metres
ymin – float minimum y bound in metres
ymax – float maximum y bound in metres

bounds¶: list of tuples lower and upper bounds for \(\beta, z_t, \Delta z, C\). \(\Delta z\) is capped where the forward model stops evaluating, which depends on the grid spacing and is hundreds of km – far beyond any Curie depth on Earth, so it never binds on data that constrain the base. Reassign this attribute to impose a tighter one, bearing in mind that a bound near the physical range will truncate the upper tail of a skewed posterior rather than report it.

prior¶: dict dictionary of priors for \(\beta, z_t, \Delta z, C\)

grid¶: 2D numpy array 2D array of magnetic data

xmin¶: float minimum x bound in metres

xmax¶: float maximum x bound in metres

ymin¶: float minimum y bound in metres

ymax¶: float maximum y bound in metres

dx¶: float grid spacing in the x-direction in metres

dy¶: float grid spacing in the y-direction in metres

nx¶: int number of nodes in the x-direction

ny¶: int number of nodes in the y-direction

xcoords¶: 1D numpy array 1D numpy array of coordinates in the x-direction

ycoords¶: 1D numpy array 1D numpy array of coordinates in the y-direction

Notes

In all instances x indicates eastings in metres and y indicates northings. Using a grid of longitude / latitudinal coordinates (degrees) will result in incorrect Curie depth calculations.

add_prior(**kwargs)[source]¶

Add a prior to the dictionary (tuple) Available priors are \(\beta, z_t, \Delta z, C\)

Assumes a normal distribution or define another distribution from scipy.stats

Usage:

>>> add_prior(beta=(p, sigma_p))

>>> add_prior(beta=scipy.stats.norm(p, sigma_p))

reset_priors()[source]¶: Reset priors to uniform distribution

objective_routine(**kwargs)[source]¶

Evaluate the objective routine to find the misfit with priors Only keys carrying a prior will be added to the total misfit

Parameters:: kwargs – parameter values to test against their priors

Usage:

>>> objective_routine(beta=2.5)

Returns:

misfit –

float: misfit integrated over all observations and priors

objective_function(x, x0, sigma_x0, *args)[source]¶

Objective function used in objective_routine Evaluates the l2-norm misfit

Parameters:

x – float, ndarray
x0 – float, ndarray
sigma_x0 – float, ndarray

Returns:

misfit – float

residuals(x, kh, Phi, sigma_Phi, prior=None)[source]¶

Whitened residuals of the fit: the spectrum first, then one entry per prior.

min_func is the half sum of squares of this vector, and the fit covariance comes from its Jacobian, so the two cannot drift apart. A Gaussian prior \(N(p, \sigma_p)\) on a parameter \(m\) is just another observation, contributing a residual \((m - p)/\sigma_p\).

Parameters:

x – array shape (4,) \(\beta, z_t, \Delta z, C\)
kh – array shape (n,) wavenumbers (rad/km)
Phi – array shape (n,) radial power spectrum \(\Phi\)
sigma_Phi – array shape (n,) uncertainty of \(\Phi\), as returned by pycurious.grid.CurieGrid.window_spectrum
prior – dict, optional priors to use in place of self.prior

Returns:

residuals – array shape (n + number of priors,)

Notes

Warnings from pycurious.grid.bouligand2009 are suppressed because some combinations of parameters overflow, which would otherwise crash the minimiser. Any residual that comes back non-finite is replaced by a large finite value, so that one unusable bin costs the fit a fixed penalty rather than poisoning the whole vector.

min_func(x, kh, Phi, sigma_Phi, prior=None)[source]¶

Function to minimise: the negative log posterior, up to a constant.

Parameters:

x – array shape (4,) array of variables \(\beta, z_t, \Delta z, C\)
kh – array shape (n,) wavenumbers (rad/km)
Phi – array shape (n,) radial power spectrum \(\Phi\)
sigma_Phi – array shape (n,) uncertainty of \(\Phi\), as returned by pycurious.grid.CurieGrid.window_spectrum
prior – dict, optional priors to use in place of self.prior

Returns:

misfit –

float: sum of misfit (scalar)

Notes

sigma_Phi should be the uncertainty of the binned mean, not the scatter of the cells within each annulus. Weighting by the latter recovers parameters no better than not weighting at all, because it is nearly flat across the spectrum and so carries almost no information – see pycurious.grid.CurieGrid.window_spectrum.

optimise(window, xc, yc, beta=3.0, zt=1.0, dz=10.0, C=5.0, taper=<function hanning>, process_subgrid=None, dof_factor=None, return_cov=False, **kwargs)[source]¶

Find the optimal parameters of \(\beta, z_t, \Delta z, C\) for a given centroid (xc,yc) and window size, with their uncertainties.

Parameters:

window – float size of window in metres
xc – float centroid x values
yc – float centroid y values
beta – float fractal parameter (starting value)
zt – float top of magnetic layer (starting value)
dz – float thickness of magnetic layer (starting value)
C – float field constant (starting value)
taper – taper (default=`numpy.hanning`) taper function, set to None for no taper function
process_subgrid – function a custom function to process the subgrid
dof_factor – float, optional override the effective-degrees-of-freedom deflation applied to the spectral uncertainties, see pycurious.grid.CurieGrid.window_spectrum
return_cov – bool (default=False) also return the 4x4 parameter covariance matrix
kwargs – keyword arguments to pass to radial_spectrum.

Returns:

beta –

float: fractal parameter
ztfloat: top of magnetic layer
dzfloat: thickness of magnetic layer
Cfloat: field constant
sigma_beta, sigma_zt, sigma_dz, sigma_Cfloat: standard deviation of each of the above
cov2D array shape (4,4): parameter covariance. Only returned if return_cov=True.

Usage:

>>> beta, zt, dz, C, s_beta, s_zt, s_dz, s_C = grid.optimise(
...     200e3, xc, yc)
>>> CPD, sigma_CPD = grid.calculate_CPD(zt, dz, s_zt, s_dz)

Notes

The uncertainties come from the curvature of the misfit at the solution, corrected for the correlation between neighbouring spectral bins – see _covariance. They describe the scatter of the spectrum at this window and this model. They do not include the systematic error from the choice of window size or centroid, which on a small grid is larger: sweeping those over tests/test_mag_data.txt moves \(\Delta z\) by 5.5 km, where the fit reports about 3.3 km.

sigma_dz in particular should be read as a lower bound. The likelihood in \(\Delta z\) has a long upper tail (Mather & Fullea, 2019), so a symmetric interval is the wrong shape for it. Measured over 200 independent synthetics, sigma_beta, sigma_zt and sigma_C reproduce the true spread to within 1%, while sigma_dz understates it by about 40%. Use profile for an honest interval on \(\Delta z\) and on the Curie depth.

optimise_routine(window, xc_list, yc_list, beta=3.0, zt=1.0, dz=10.0, C=5.0, taper=<function hanning>, process_subgrid=None, dof_factor=None, **kwargs)[source]¶

Iterate through a list of centroids to compute the optimal values of \(\beta, z_t, \Delta z, C\) for a given window size.

Parameters:

window – float size of window in metres
xc_list – ndarray shape (l,) centroid x values
yc_list – ndarray shape (l,) centroid y values
beta – float fractal parameter
zt – float top of magnetic layer
dz – float thickness of magnetic layer
C – float field constant
taper – function taper function (default=`numpy.hanning`) set to None for no taper function
process_subgrid – func a custom function to process the subgrid
dof_factor – float, optional override the effective-degrees-of-freedom deflation applied to the spectral uncertainties, see pycurious.grid.CurieGrid.window_spectrum
kwargs – keyword arguments to pass to radial_spectrum.

Returns:

beta –

ndarray shape (l,): fractal parameters
ztndarray shape (l,): top of magnetic layer
dzndarray shape (l,): thickness of magnetic layer
Cndarray shape (l,): field constant
sigma_beta, sigma_zt, sigma_dz, sigma_Cndarray shape (l,): standard deviation of each of the above, so a map of the uncertainty comes out alongside the map of the parameter

Notes

The covariance matrix is deliberately not available here. pycurious.parallel.CurieParallel.parallelise_routine collects one array per returned quantity, which a 4x4 matrix per centroid does not fit. Call optimise directly with return_cov=True for that.

profile(window, xc, yc, target, level=0.95, npoints=21, bracket=None, beta=3.0, zt=1.0, dz=10.0, C=5.0, taper=<function hanning>, process_subgrid=None, dof_factor=None, **kwargs)[source]¶

Confidence interval for one parameter, or for the Curie depth, without assuming the posterior is symmetric.

Each point of the scan holds target fixed and re-optimises everything else, tracing the deviance \(2(F - F_{min})\). The interval is where that crosses \(\chi^2_1\) at the requested level, which is the usual likelihood-ratio construction.

This matters most for \(\Delta z\) and hence the Curie depth. Both have a long upper tail (Mather & Fullea, 2019), so the symmetric \(\pm \sigma\) that optimise reports understates how far the parameter can plausibly reach – by about 40% on synthetics.

Parameters:

window – float size of window in metres
xc – float centroid of the window
yc – float centroid of the window
target – str one of "beta", "zt", "dz", "C" or "CPD"
level – float (default=0.95) confidence level
npoints – int (default=21) nodes in the scan. The endpoints are then refined by root finding, so this sets the resolution of the returned curve rather than of the interval.
bracket – tuple, optional (min, max) of the scan. Defaults to a range either side of the fitted value, wider above than below because of the tail.
beta – float starting values for the underlying fit
zt – float starting values for the underlying fit
dz – float starting values for the underlying fit
C – float starting values for the underlying fit
taper – function (default=np.hanning) taper function, or None for no taper
process_subgrid – function, optional applied to the subgrid before the spectrum is computed
dof_factor – float, optional see pycurious.grid.CurieGrid.window_spectrum
kwargs – keyword arguments passed to radial_spectrum

Returns:

values –

1D array shape (npoints,): where the target was held
deviance1D array shape (npoints,): \(2(F - F_{min})\) at each of those
lowerfloat: lower end of the interval, -inf if the scan never crossed
upperfloat: upper end of the interval, inf if the scan never crossed

Usage:

>>> values, deviance, lo, hi = grid.profile(200e3, xc, yc, "CPD")
>>> print("Curie depth {:.1f} ({:.1f} to {:.1f}) km".format(cpd, lo, hi))

Notes

This is a profile posterior deviance rather than a profile likelihood: any priors added with add_prior contribute to F. A Gaussian prior is one more observation, so the calibration still holds, but it is not the marginal an MCMC would report – profiling takes the ridge of the posterior rather than integrating over it, and so is a little narrower for a skewed one.

Like the covariance from optimise, the interval describes the scatter of the spectrum at a fixed window, centroid and model. It does not cover the systematic error from choosing those: on tests/test_mag_data.txt the interval for \(\Delta z\) is about 3.3 km wide, where sweeping the window size and centroid moves \(\Delta z\) over 5.5 km.

metropolis_hastings(window, xc, yc, nsim, burnin, x_scale=None, beta=3.0, zt=1.0, dz=10.0, C=5.0, taper=<function hanning>, process_subgrid=None, dof_factor=None, adapt=True, seed=None, return_diagnostics=False, **kwargs)[source]¶

MCMC algorithm using a Metropolis-Hastings sampler.

Evaluates a Markov chain for starting values of \(\beta, z_t, \Delta z, C\) and returns the ensemble of model realisations.

Parameters:

window – float size of window in metres
xc – float centroid x values
yc – float centroid y values
nsim – int number of simulations
burnin – int number of burn-in simulations discarded before the nsim recorded samples
x_scale – float(4), optional initial width of the proposal in each parameter (default=`[1,1,1,1]` for [beta, zt, dz, C]). With adapt=True this is only a starting point.
beta – float fractal parameter (starting value for the search)
zt – float top of magnetic layer (starting value for the search)
dz – float thickness of magnetic layer (starting value for the search)
C – float field constant (starting value for the search)
taper – function (default=np.hanning) taper function, or None for no taper
process_subgrid – function, optional applied to the subgrid before the spectrum is computed
dof_factor – float, optional see pycurious.grid.CurieGrid.window_spectrum
adapt – bool (default=True) tune the proposal during burn-in – see Notes. Turning this off is only sensible if you have a good x_scale already.
seed – int, optional seed for reproducibility
return_diagnostics – bool (default=False) also return a dict of acceptance, burnin_acceptance, x_scale

Returns:

beta –

ndarray shape (nsim,): fractal parameter
ztndarray shape (nsim,): top of magnetic layer
dzndarray shape (nsim,): thickness of magnetic layer
Cndarray shape (nsim,): field constant
diagnosticsdict: only if return_diagnostics=True

Usage:

>>> posterior, info = grid.metropolis_hastings(
...     200e3, xc, yc, 10000, 2000, seed=1, return_diagnostics=True)
>>> print("acceptance {:.2f}".format(info["acceptance"]))

Notes

Acceptance is decided in log space. Comparing \(e^{-F}\) directly underflows to zero for any real spectrum – \(F\) runs to hundreds – at which point every proposal is rejected and the chain returns a handful of distinct states dressed up as a posterior.

The chain starts at the mode, found with the same minimiser optimise uses and from the same starting values. That costs a fraction of a second and removes the job the burn-in is worst at.

There is no tempering. It was tried – annealing the burn-in after Sambridge (2013), doi:10.1093/gji/ggt342 – and made every case worse. What motivated it was that large parts of the posterior evaluated to zero, and that was the \(e^{-F}\) underflow rather than a property of the problem, so log-space acceptance removes the reason for it. It also fights the proposal tuning below: a high temperature makes almost everything acceptable, driving the scale up, and the scale then collapses as the temperature falls, freezing the chain wherever the hot phase left it.

The shape of the proposal matters more than any of the above. The four parameters are strongly correlated – \(\beta\) with \(z_t\) at about -0.92, \(z_t\) with \(C\) at about 0.87 – and their marginal widths differ by a factor of thirty, so a proposal with one width per parameter cannot move along the ridge they lie on, and the chain sits still. The proposal is therefore drawn along the fit covariance, the same one optimise reports, with the burn-in tuning only a scalar multiplier on it towards an acceptance rate of 0.234 (Robbins-Monro). x_scale sets where that multiplier starts, and adapt=False fixes it there.

The chain respects self.bounds, which the optimiser has always done but the sampler previously did not.

Both this and sensitivity treat the spectral bins as independent, which they are not, so the posterior is narrower than the spread over independent realisations of the field. See optimise.

sensitivity(window, xc, yc, nsim, beta=3.0, zt=1.0, dz=10.0, C=5.0, taper=<function hanning>, process_subgrid=None, dof_factor=None, seed=None, **kwargs)[source]¶

Sample the uncertainty of \(\beta, z_t, \Delta z, C\) by resampling the spectrum, and the centre of each prior distribution (if provided by the user - see add_prior).

Parameters:

window – float size of window in metres
xc – float centroid x values
yc – float centroid y values
nsim – int number of Monte Carlo simulations
beta – float starting fractal parameter
zt – float starting top of magnetic layer
dz – float starting thickness of magnetic layer
C – float starting field constant
taper – function (default=`numpy.hanning`) taper function, set to None for no taper function
process_subgrid – function, optional a custom function to process the subgrid
dof_factor – float, optional override the effective-degrees-of-freedom deflation, see pycurious.grid.CurieGrid.window_spectrum
seed – int, optional seed for reproducibility

Returns:

beta –

ndarray shape (nsim,): fractal parameters
ztndarray shape (nsim,): top of magnetic layer
dzndarray shape (nsim,): thickness of magnetic layer
Cndarray shape (nsim,): field constant

Notes

Each bin of the spectrum is resampled independently, so this shares the assumption behind the fit covariance that the bins are independent. They are not – a taper correlates neighbouring annuli – so agreement between the two is not evidence that either is right. Only an ensemble over independent realisations of the field calibrates that.

calculate_CPD(zt, dz, sigma_zt=0.0, sigma_dz=0.0)[source]¶

Compute the Curie depth from the results of optimise.

Parameters:

zt – float / 1D array depth to the top of the magnetic source
dz – float / 1D array thickness of the magnetic source
sigma_zt – float / 1D array standard deviation of zt
sigma_dz – float / 1D array standard deviation of dz

Returns:

CPD –

float / 1D array: estimated Curie point depth at the base of the magnetic source
CPD_stdevfloat / 1D array: standard deviation of CPD

Usage:

>>> beta, zt, dz, C, s_beta, s_zt, s_dz, s_C = grid.optimise(
...     200e3, xc, yc)
>>> CPD, sigma_CPD = grid.calculate_CPD(zt, dz, s_zt, s_dz)

Notes

\(Z_b = z_t + \Delta z\), so the uncertainties combine as \(\sqrt{\sigma_{z_t}^2 + \sigma_{\Delta z}^2}\). The two are correlated – about 0.6 – but \(\sigma_{\Delta z}\) exceeds \(\sigma_{z_t}\) by four orders of magnitude, so including the covariance changes the answer by around 1%. optimise will hand over the full matrix with return_cov=True for anyone who wants it.

The far larger effect is that this is symmetric and the Curie depth is not: \(\Delta z\) has a long upper tail, so CPD_stdev understates how deep the base can plausibly lie. Use profile with target="CPD" for an interval that does not assume symmetry.

Matches the signature and return of pycurious.optimise_tanaka.CurieOptimiseTanaka.calculate_CPD, so code written against one behaves the same against the other.

create_centroid_list(window, spacingX=None, spacingY=None)¶

Create a list of xc,yc values to extract subgrids.

Parameters:

window – float size of the windows in metres
spacingX – float (optional) specify spacing in metres in the X direction will default to maximum X resolution
spacingY – float (optional) specify spacing in metres in the Y direction will default to maximum Y resolution

Returns:

xc_list –

1D array: array of x coordinates
yc_list1D array: array of y coordinates

parallelise_routine(window, xc_list, yc_list, func, *args, **kwargs)¶

Implements shared memory multiprocessing to split multiple evaluations of a function centroids across processors.

Supply the window size and lists of x,y coordinates to a function along with any additional arguments or keyword arguments.

Parameters:

window – float size of window in metres
xc_list – array shape (l,) centroid x values
yc_list – array shape (l,) centroid y values
func – function Python function to evaluate in parallel
args – arguments additional arguments to pass to func
kwargs –
keyword arguments additional keyword arguments to pass to func. Two keys are reserved and consumed here rather than forwarded:
- on_error : {“raise”, “ignore”} (default=”raise”) what to do when func fails for a centroid. "ignore" fills that centroid with NaN and warns.
- seed : int or None (default=None) seeds any stochastic routine reproducibly. Each centroid receives an independent child seed derived from numpy.random.SeedSequence, passed to func as a seed keyword, so results do not depend on how many processors were used. func must accept a seed keyword if this is supplied.

Returns:

out –

list of lists: (depends on output of func - see notes)

Usage:

An obvious use case is to compute the Curie depth for many centroids in parallel.

>>> self.parallelise_routine(window, xc_list, yc_list, self.optimise)

Each centroid is assigned a new process and sent to a free processor to compute. In this case, the output is separate lists of shape(l,) for \(\beta, z_t, \Delta z, C\). If len(xc_list)=2 then,

>>> self.parallelise_routine(window, [x1,x2], [y1, y2], self.optimise)
[[beta1  beta2], [zt1  zt2], [dz1  dz2], [C1  C2]]

Another example is to parallelise the sensitivity analysis:

>>> self.parallelise_routine(window, xc_list, yc_list, self.sensitivity, nsim)

>>> self.parallelise_routine(window, [x1,x2], [y1,y2], self.sensitivity, 4)

which would return:

[[[ beta1a , beta1b , beta1c , beta1d ],   # centroid 1 (x1,y1)
  [ beta2a , beta2b , beta2c , beta2d ]],  # centroid 2 (x2,y2)
 [[   zt1a ,   zt1b ,   zt1c ,   zt1d ],   # centroid 1 (x1,y1)
  [   zt2a ,   zt2b ,   zt2c ,   zt2d ]],  # centroid 2 (x2,y2)
 [[   dz1a ,   dz1b ,   dz1c ,   dz1d ],   # centroid 1 (x1,y1)
  [   dz2a ,   dz2b ,   dz2c ,   dz2d ]]   # centroid 2 (x2,y2)
 [[    C1a ,    C1b ,    C1c ,    C1d ],   # centroid 1 (x1,y1)
  [    C2a ,    C2b ,    C2c ,    C2d ]]]  # centroid 2 (x2,y2)

Notes

See the module docstring for the if __name__ == "__main__": requirement when calling this from a script.

radial_spectrum(subgrid, taper=<function hanning>, power=2.0, return_counts=False, **kwargs)¶

Compute the radial spectrum for a square grid.

Wavenumber is returned in units of rad/km.

Parameters:

subgrid – 2D array window of the original data (see subgrid method)
taper – function (default=np.hanning) taper function, set to None for no taper function
power –
float raise the FFT of the magnetic anomaly to the power:
- 2.0 for Bouligand et al. (2009) use cases, which gives the log power spectrum \(\ln \Phi_{\Delta T}\)
- 1.0 for Tanaka et al. (1999) use cases, which gives the log amplitude spectrum \(\ln \Phi_{\Delta T}^{1/2}\)
return_counts – bool (default=False) also return the number of FFT cells averaged into each bin
kwargs – keyword arguments keyword arguments to pass to taper

Returns:

k –

1D array shape (n,): wavenumber in rad/km
Phi1D array shape (n,): Radial power spectrum
sigma_Phi1D array shape (n,): Standard deviation of Phi within each radial bin
counts1D array shape (n,): number of FFT cells in each radial bin. Only returned if return_counts=True.

Notes

References

reduce_to_pole(data, inc, dec, sinc=None, sdec=None)¶

Reduce total field magnetic anomaly data to the pole.

Parameters:

data – 1D array the total field anomaly data at each point.
inc – float / 1D array inclination of the inducing Geomagnetic field
dec – float / 1D array declination of the inducing Geomagnetic field
sinc – float / 1D array (optional) inclination of the total magnetization of the anomaly source
sdec – float / 1D array (optional) 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 –

2D array: the data reduced to the pole.

References

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

Notes

This functions performs the reduction in the frequency domain (using the FFT). The transform filter is (in the freq 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, k_y\) are the wave-numbers in the x and y directions and

\(|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\)

\(\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.

remove_trend_linear(data)¶

Remove the best-fitting linear trend from the data

This may come in handy if the magnetic data has not been reduced to the pole.

Parameters:: data – 2D numpy array
Returns:: data – 2D numpy array

subgrid(window, xc, yc)¶

Extract a subgrid from the data at a window around the point (xc,yc)

Parameters:

xc – float x coordinate
yc – float y coordinate
window – float size of window in metres

Returns:

data –

2D array: subgrid encompassing window size

upward_continuation(data, 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, \(|k|\) is the wavenumber modulus.

Parameters:

data – 2D array potential field at the grid points
height – float height increase (delta z) in meters.

Returns:

cont –

array: upward continued data

References

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

window_spectrum(window, xc, yc, taper=<function hanning>, power=2.0, process_subgrid=None, dof_factor=None, **kwargs)¶

Radial spectrum of one window, weighted ready for fitting.

Parameters:

window – float size of the window in metres
xc – float centroid of the window
yc – float centroid of the window
taper – function (default=np.hanning) taper function, or None for no taper
power – float raise the FFT of the anomaly to this power – 2.0 for the log power spectrum that Bouligand et al. (2009) fit, 1.0 for the log amplitude spectrum of Tanaka et al. (1999)
process_subgrid – function, optional applied to the subgrid before the spectrum is computed
dof_factor – float, optional override the effective-degrees-of-freedom deflation (see Notes)
kwargs – keyword arguments passed to taper

Returns:

k –

1D array: wavenumber in rad/km
Phi1D array: log spectrum, raised to power
sigma1D array: uncertainty of the binned mean

Usage:

>>> k, Phi, sigma = grid.window_spectrum(200e3, xc, yc, power=2)

Notes

Tanaka optimiser¶

class pycurious.CurieOptimiseTanaka(grid, xmin, xmax, ymin, ymax, **kwargs)[source]¶

Bases: CurieGrid

Extends pycurious.grid.CurieGrid with the centroid method of Tanaka et al. (1999).

Two straight lines are fitted to separate bands of the radial amplitude spectrum, giving the top and centroid depths of the magnetic source and, from those, the Curie point depth with an uncertainty.

Parameters:

grid – 2D array 2D array of magnetic data
xmin – float minimum/maximum x bounds of the grid
xmax – float minimum/maximum x bounds of the grid
ymin – float minimum/maximum y bounds of the grid
ymax – float minimum/maximum y bounds of the grid
max_processors – int, optional processors to use in optimise_routine (default=all)

max_processors¶: int processors used by the parallel routines

Notes

All attributes of pycurious.grid.CurieGrid are inherited as well.

The grid must be projected in eastings/northings (metres), not degrees, since depths are computed in km.

check_bands(k, zt_range, z0_range, thickness=None, verbose=True)[source]¶

Report whether two fitting bands are usable, before relying on them.

Both of Tanaka’s straight-line approximations hold only over part of the spectrum. A band outside that range still produces a confident looking fit, so this is worth checking explicitly.

Parameters:

k – 1D array wavenumbers from radial_spectrum, in rad/km
zt_range – tuple (min, max) wavenumber of the \(Z_t\) band, rad/km
z0_range – tuple (min, max) wavenumber of the \(Z_0\) band, rad/km
thickness – float, optional estimated source thickness in km. Without it the validity of each approximation cannot be assessed, only the point counts.
verbose – bool (default=True) print a summary. Set False to receive any problems as UserWarning instead, for use in a script.

Returns:

diagnostics –

dict: dk, n_zt, n_z0, lambda_zt, lambda_z0 and, if thickness was given, kd_max, CPD_bias and lambda_zt_min_required. CPD_bias is an estimate in km of how far the \(|k|d\) approximation drags the Curie depth, and is negative.

Usage:

>>> k, Phi, sigma_Phi = grid.radial_spectrum(subgrid, power=1)
>>> grid.check_bands(k, (0.2, 0.6), (0.0, 0.05), thickness=20.0)

optimise(window, xc, yc, zt_range, z0_range, taper=<function hanning>, beta=None, process_subgrid=None, absolute_sigma=True, dof_factor=None, **kwargs)[source]¶

Estimate the top and centroid depths of the magnetic source for one centroid, with their uncertainties.

Parameters:

window – float size of the window in metres
xc – float centroid of the window
yc – float centroid of the window
zt_range – tuple (min, max) wavenumber in rad/km over which to fit \(Z_t\). Should cover wavelengths shorter than twice the source thickness.
z0_range – tuple (min, max) wavenumber in rad/km over which to fit \(Z_0\). Must satisfy \(|k| d \ll 1\) – check with check_bands.
taper – function (default=np.hanning) taper function, or None for no taper
beta – float, optional fractal parameter of the magnetisation. If given, its contribution is removed before fitting. Leave unset for the method exactly as Tanaka published it; beta=1 is equivalent.
process_subgrid – function, optional applied to the subgrid before the spectrum is computed
absolute_sigma – bool (default=True) treat the spectral uncertainties as absolute, so the reported errors carry their units. Set False to rescale the covariance by the reduced chi-squared instead.
dof_factor – float, optional override the effective-degrees-of-freedom deflation applied to the spectral uncertainties (see Notes)
kwargs – keyword arguments passed to radial_spectrum

Returns:

zt –

float: depth to the top of the magnetic source, km
z0float: centroid depth of the magnetic source, km
zt_interceptfloat: intercept of the \(Z_t\) fit
z0_interceptfloat: intercept of the \(Z_0\) fit
sigma_ztfloat: standard deviation of zt
sigma_z0float: standard deviation of z0

Usage:

>>> zt, z0, zt_i, z0_i, sigma_zt, sigma_z0 = grid.optimise(
...     200e3, xc, yc, (0.2, 0.6), (0.0, 0.05))
>>> CPD, sigma_CPD = grid.calculate_CPD(
...     zt, z0, sigma_zt=sigma_zt, sigma_z0=sigma_z0)

Notes

Depths are returned positive downwards, i.e. the negated gradient of each fit.

The reported uncertainties describe the scatter of the spectrum only. They do not include the systematic error from the choice of band, which is usually larger – see sensitivity.

radial_spectrum returns the scatter of the FFT cells within each annulus, whereas the fit needs the uncertainty of the annulus mean. That is the standard error, except that the cells are not independent: Hermitian symmetry makes about half of them redundant, and tapering correlates neighbours. The correction is calibrated per taper; dof_factor overrides it.

Cells in neighbouring annuli are correlated too, which _fit_band allows for. Against 80 independent synthetics the ratio of the true spread to the reported sigma_zt improves from 1.48 to 1.12 with that correction in place.

sigma_z0 remains understated, at a ratio of about 1.34, and no covariance can fix it. The centroid gradient is fitted over a handful of the longest wavelengths the window resolves, and its distribution is heavy-tailed: on a 4000 km grid with a true \(Z_0\) of 11 km, the middle 90% of estimates spanned 4.2 to 25.3 km. Treat sigma_z0, and the Curie depth that follows from it, as a lower bound.

There is no profile-likelihood alternative here, as there is on the Bouligand side. Each band is a straight-line fit, so its misfit is exactly quadratic in the gradient and the profile interval is provably the same as the covariance one – verified, a deviance of 3.841459 against a threshold of 3.841459 on both bands. It would return the number _fit_band already returns. The equivalence holds only for absolute_sigma=True; setting it False rescales the covariance by the reduced chi-squared afterwards, which the profile construction does not do.

optimise_routine(window, xc_list, yc_list, zt_range, z0_range, taper=<function hanning>, beta=None, process_subgrid=None, absolute_sigma=True, dof_factor=None, **kwargs)[source]¶

Iterate optimise over a list of centroids, in parallel.

Takes the same arguments as optimise, with lists of centroids in place of a single one. See pycurious.parallel.CurieParallel.parallelise_routine for the on_error and seed keywords.

Returns:: zt, z0, zt_intercept, z0_intercept, sigma_zt, sigma_z0 – 1D arrays, one entry per centroid

Usage:

>>> xc_list, yc_list = grid.create_centroid_list(window, 10e3, 10e3)
>>> zt, z0, zt_i, z0_i, sigma_zt, sigma_z0 = grid.optimise_routine(
...     window, xc_list, yc_list, (0.2, 0.6), (0.0, 0.05))

sensitivity(window, xc, yc, nsim, zt_range, z0_range, taper=<function hanning>, beta=None, band_scale=0.1, process_subgrid=None, absolute_sigma=True, dof_factor=None, seed=None, **kwargs)[source]¶

Sample the uncertainty of the Curie depth by perturbing both the spectrum and the fitting bands.

The fit covariance alone understates the uncertainty, because where the band edges are placed usually matters more than the scatter of the spectrum. Each simulation therefore redraws the spectrum within its uncertainty and jitters both band edges. The remaining arguments are as optimise.

Parameters:

window – float size of the window in metres
xc – float centroid of the window
yc – float centroid of the window
nsim – int number of simulations
zt_range – tuple as optimise, in rad/km. These are the centres about which the band edges are jittered.
z0_range – tuple as optimise, in rad/km. These are the centres about which the band edges are jittered.
band_scale – float (default=0.1) standard deviation of the jitter applied to each band edge, as a fraction of that band’s width. Set to 0 to perturb only the spectrum, which recovers the analytic covariance.
seed – int, optional seed for reproducibility

Returns:

zt –

1D array shape (nsim,): sampled top depths
z01D array shape (nsim,): sampled centroid depths
CPD1D array shape (nsim,): sampled Curie point depths

Usage:

>>> zt, z0, CPD = grid.sensitivity(
...     200e3, xc, yc, 500, (0.2, 0.6), (0.0, 0.05))
>>> print(CPD.mean(), CPD.std())

Notes

This samples the statistical uncertainty and the sensitivity to band placement. It does not capture the systematic error from fitting outside the range where each approximation holds, nor from unmodelled fractal magnetisation – both bias the two fits coherently rather than scattering them. Use check_bands and beta for those.

calculate_CPD(zt, z0, sigma_zt=0.0, sigma_z0=0.0)[source]¶

Compute the Curie depth from the results of optimise.

Parameters:

zt – float / 1D array depth to the top of the magnetic source
z0 – float / 1D array centroid depth of the magnetic source
sigma_zt – float / 1D array standard deviation of zt
sigma_z0 – float / 1D array standard deviation of z0

Returns:

CPD –

float / 1D array: estimated Curie point depth at the base of the magnetic source
CPD_stdevfloat / 1D array: standard deviation of CPD

Notes

\(Z_b = 2 Z_0 - Z_t\), so the uncertainties combine as \(\sqrt{\sigma_{Z_t}^2 + 4\sigma_{Z_0}^2}\). This assumes the two fits are independent, which holds well enough in practice – they use disjoint bands, and the measured correlation between them is about 0.05.

zt and z0 are expected positive downwards, as returned by optimise.

create_centroid_list(window, spacingX=None, spacingY=None)¶

Create a list of xc,yc values to extract subgrids.

Parameters:

window – float size of the windows in metres
spacingX – float (optional) specify spacing in metres in the X direction will default to maximum X resolution
spacingY – float (optional) specify spacing in metres in the Y direction will default to maximum Y resolution

Returns:

xc_list –

1D array: array of x coordinates
yc_list1D array: array of y coordinates

parallelise_routine(window, xc_list, yc_list, func, *args, **kwargs)¶

Implements shared memory multiprocessing to split multiple evaluations of a function centroids across processors.

Supply the window size and lists of x,y coordinates to a function along with any additional arguments or keyword arguments.

Parameters:

window – float size of window in metres
xc_list – array shape (l,) centroid x values
yc_list – array shape (l,) centroid y values
func – function Python function to evaluate in parallel
args – arguments additional arguments to pass to func
kwargs –
keyword arguments additional keyword arguments to pass to func. Two keys are reserved and consumed here rather than forwarded:
- on_error : {“raise”, “ignore”} (default=”raise”) what to do when func fails for a centroid. "ignore" fills that centroid with NaN and warns.
- seed : int or None (default=None) seeds any stochastic routine reproducibly. Each centroid receives an independent child seed derived from numpy.random.SeedSequence, passed to func as a seed keyword, so results do not depend on how many processors were used. func must accept a seed keyword if this is supplied.

Returns:

out –

list of lists: (depends on output of func - see notes)

Usage:

An obvious use case is to compute the Curie depth for many centroids in parallel.

>>> self.parallelise_routine(window, xc_list, yc_list, self.optimise)

Each centroid is assigned a new process and sent to a free processor to compute. In this case, the output is separate lists of shape(l,) for \(\beta, z_t, \Delta z, C\). If len(xc_list)=2 then,

>>> self.parallelise_routine(window, [x1,x2], [y1, y2], self.optimise)
[[beta1  beta2], [zt1  zt2], [dz1  dz2], [C1  C2]]

Another example is to parallelise the sensitivity analysis:

>>> self.parallelise_routine(window, xc_list, yc_list, self.sensitivity, nsim)

>>> self.parallelise_routine(window, [x1,x2], [y1,y2], self.sensitivity, 4)

which would return:

[[[ beta1a , beta1b , beta1c , beta1d ],   # centroid 1 (x1,y1)
  [ beta2a , beta2b , beta2c , beta2d ]],  # centroid 2 (x2,y2)
 [[   zt1a ,   zt1b ,   zt1c ,   zt1d ],   # centroid 1 (x1,y1)
  [   zt2a ,   zt2b ,   zt2c ,   zt2d ]],  # centroid 2 (x2,y2)
 [[   dz1a ,   dz1b ,   dz1c ,   dz1d ],   # centroid 1 (x1,y1)
  [   dz2a ,   dz2b ,   dz2c ,   dz2d ]]   # centroid 2 (x2,y2)
 [[    C1a ,    C1b ,    C1c ,    C1d ],   # centroid 1 (x1,y1)
  [    C2a ,    C2b ,    C2c ,    C2d ]]]  # centroid 2 (x2,y2)

Notes

See the module docstring for the if __name__ == "__main__": requirement when calling this from a script.

radial_spectrum(subgrid, taper=<function hanning>, power=2.0, return_counts=False, **kwargs)¶

Compute the radial spectrum for a square grid.

Wavenumber is returned in units of rad/km.

Parameters:

subgrid – 2D array window of the original data (see subgrid method)
taper – function (default=np.hanning) taper function, set to None for no taper function
power –
float raise the FFT of the magnetic anomaly to the power:
- 2.0 for Bouligand et al. (2009) use cases, which gives the log power spectrum \(\ln \Phi_{\Delta T}\)
- 1.0 for Tanaka et al. (1999) use cases, which gives the log amplitude spectrum \(\ln \Phi_{\Delta T}^{1/2}\)
return_counts – bool (default=False) also return the number of FFT cells averaged into each bin
kwargs – keyword arguments keyword arguments to pass to taper

Returns:

k –

1D array shape (n,): wavenumber in rad/km
Phi1D array shape (n,): Radial power spectrum
sigma_Phi1D array shape (n,): Standard deviation of Phi within each radial bin
counts1D array shape (n,): number of FFT cells in each radial bin. Only returned if return_counts=True.

Notes

References

reduce_to_pole(data, inc, dec, sinc=None, sdec=None)¶

Reduce total field magnetic anomaly data to the pole.

Parameters:

data – 1D array the total field anomaly data at each point.
inc – float / 1D array inclination of the inducing Geomagnetic field
dec – float / 1D array declination of the inducing Geomagnetic field
sinc – float / 1D array (optional) inclination of the total magnetization of the anomaly source
sdec – float / 1D array (optional) 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 –

2D array: the data reduced to the pole.

References

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

Notes

This functions performs the reduction in the frequency domain (using the FFT). The transform filter is (in the freq 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, k_y\) are the wave-numbers in the x and y directions and

\(|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\)

\(\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.

remove_trend_linear(data)¶

Remove the best-fitting linear trend from the data

This may come in handy if the magnetic data has not been reduced to the pole.

Parameters:: data – 2D numpy array
Returns:: data – 2D numpy array

subgrid(window, xc, yc)¶

Extract a subgrid from the data at a window around the point (xc,yc)

Parameters:

xc – float x coordinate
yc – float y coordinate
window – float size of window in metres

Returns:

data –

2D array: subgrid encompassing window size

upward_continuation(data, 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, \(|k|\) is the wavenumber modulus.

Parameters:

data – 2D array potential field at the grid points
height – float height increase (delta z) in meters.

Returns:

cont –

array: upward continued data

References

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

window_spectrum(window, xc, yc, taper=<function hanning>, power=2.0, process_subgrid=None, dof_factor=None, **kwargs)¶

Radial spectrum of one window, weighted ready for fitting.

Parameters:

window – float size of the window in metres
xc – float centroid of the window
yc – float centroid of the window
taper – function (default=np.hanning) taper function, or None for no taper
power – float raise the FFT of the anomaly to this power – 2.0 for the log power spectrum that Bouligand et al. (2009) fit, 1.0 for the log amplitude spectrum of Tanaka et al. (1999)
process_subgrid – function, optional applied to the subgrid before the spectrum is computed
dof_factor – float, optional override the effective-degrees-of-freedom deflation (see Notes)
kwargs – keyword arguments passed to taper

Returns:

k –

1D array: wavenumber in rad/km
Phi1D array: log spectrum, raised to power
sigma1D array: uncertainty of the binned mean

Usage:

>>> k, Phi, sigma = grid.window_spectrum(200e3, xc, yc, power=2)

Notes

Synthetic and analytic spectra¶

pycurious.fractal_anomaly(n=512, dx=2.0, beta=3.0, zt=1.0, dz=20.0, C=5.0, seed=0)[source]¶

Synthesise a magnetic anomaly with a prescribed radial power spectrum.

Parameters:

n – int number of points per side
dx – float grid spacing in km
beta – float fractal parameter of the magnetisation
zt – float depth to the top of the magnetic source, km
dz – float thickness of the magnetic source, km
C – float field constant
seed – int seed for the white noise realisation

Returns:

data –

2D array shape (n,n): the magnetic anomaly
extenttuple: (xmin, xmax, ymin, ymax) in metres, for pycurious.grid.CurieGrid

Usage:

>>> data, extent = pycurious.fractal_anomaly(beta=3.0, zt=1.0, dz=20.0)
>>> grid = pycurious.CurieOptimiseBouligand(data, *extent)

Notes

The Curie depth of the result is zt + dz, and its centroid depth (as sought by the Tanaka method) is zt + dz/2. Those two, and beta, are recoverable.

C is not, quite. It is a level rather than a depth, and it absorbs every constant factor between the noise and the spectrum. Normalising the transform removes the dependence on n, but two offsets remain and both depend on how the spectrum is measured rather than on how the field was made:

the radial spectrum averages \(\ln |FFT|\) rather than taking the log of the mean, which is lower by the Euler-Mascheroni constant, 0.577;
a taper removes power, by \(\ln (3/8)^2 = -1.96\) for a separable numpy.hanning.

So a fit through numpy.hanning returns C about 2.5 low, and through no taper about 0.6 low. Treat the recovered C as a nuisance parameter, not as something with a true value to check against.

pycurious.bouligand2009(kh, beta, zt, dz, C)[source]¶

Calculate the synthetic radial power spectrum of magnetic anomalies

Equation (4) of Bouligand et al. (2009)

Parameters:

kh – float / 1D array wavenumber in rad/km
beta – float / 1D array fractal parameter
zt – float / 1D array top of magnetic sources
dz – float / 1D array thickness of magnetic sources
C – float / 1D array field constant (Maus et al., 1997)

Returns:

Phi –

float / 1D array: radial power spectrum of magnetic anomalies

References

Maus, S., D. Gordon, and D. Fairhead (1997), Curie temperature depth estimation using a self-similar magnetization model, Geophys. J. Int., 129, 163-168, doi:10.1111/j.1365-246X.1997.tb00945.x

pycurious.maus1995(beta, zt, kh, C=0.0)[source]¶

Calculate the synthetic radial power spectrum of magnetic anomalies (Maus and Dimri; 1995)

This is not all that useful except when testing overflow errors which occur for the second term in Bouligand et al. (2009).

Parameters:

beta – float / 1D array fractal parameter
zt – float / 1D array top of magnetic sources
kh – float / 1D array norm of the wave number in the horizontal plane
C – float / 1D array field constant (Maus et al., 1997)

Returns:

Phi –

float / 1D array: radial power spectrum of magnetic anomalies

References

Maus, S., D. Gordon, and D. Fairhead (1997), Curie temperature depth estimation using a self-similar magnetization model, Geophys. J. Int., 129, 163-168, doi:10.1111/j.1365-246X.1997.tb00945.x

Mapping¶

Projections, netCDF, and GeoTIFF I/O (all lazily imported — see the mapping and geotiff extras).

The pycurious.mapping module of PyCurious contains various functions to help manipulate geospatial data into common formats. It handles commonly encountered operations, such as:

Gridding scattered data points
Converting between coordinate reference systems (CRS)
Importing and exporting GeoTiff files

It requires some additional dependencies:

matplotlib – for plotting
pyproj – for transforming between different CRS
cartopy – for generating maps

Beware that most global data are georeferenced in WGS84 (EPSG: 4326). The radial power spectrum must be in rad/km, which requires a transformation from longitude / latitude to a local projection in eastings / northings.

For example, EMAG2 is a global compilation of the magnetic anomaly georeferenced in WGS84 longitude / latitude. This will need to be projected in a local CRS to use with PyCurious. If, for example, we are interested in a region across Ireland we could use the IRENET95 local CRS (EPSG: 2157),

transform_coordinates(lons, lats, epsg_in=4326, epsg_out=2157)

which would return a list of eastings and northings in IRENET95 projection.

pycurious.mapping.transform_coordinates(x, y, epsg_in, epsg_out)[source]¶

Transform between any coordinate system.

Requires pyproj – install using pip.

Parameters:

x – float / 1D array x coordinates (may be in degrees or metres/eastings)
y – float / 1D array y coordinates (may be in degrees or metres/northings)
epsg_in – int CRS of x and y coordinates
epsg_out – int CRS of output

Returns:

x_out –

float / list of floats: x coordinates projected in epsg_out
y_outfloat / list of floats: y coordinates projected in epsg_out

pycurious.mapping.convert_extent(extent_in, epsg_in, epsg_out)[source]¶

Transform extent from epsg_in to epsg_out

Parameters:

extent_in – tuple bounding box [minX, maxX, minY, maxY]
epsg_in – int CRS of extent
epsg_out – int CRS of output

Returns:

extent_out –

tuple: bounding box in new CRS

pycurious.mapping.trim(coords, data, extent, buffer_amount=0.0)[source]¶

Trim a smaller section of a large dataset taking into consideration transformations into various coordinate reference systems (CRS).

Parameters:

coords – array shape (n,2) geographical / projected coordinates
data – array shape (n,) values corresponding to coordinates
extent – tuple bounding box to trim data
buffer_amount – float amount of buffer to include (default=0.0)

Returns:

coords_trim –

array shape (l,2): trimmed coordinates
data_trimarray shape (l,): trimmed data array

pycurious.mapping.grid(coords, data, extent, shape=None, epsg_in=None, epsg_out=None, **kwargs)[source]¶

Grid a smaller section of a large dataset taking into consideration transformations into various coordinate reference systems (CRS).

Requires scipy.interpolate.griddata

Parameters:

coords – array shape (n,2) geographical coordinates
data – array shape (n,) values corresponding to coordinates
extent – tuple bounding box in espg_out coordinates
shape – tuple (nrows,ncols) size of the box, if None, shape is estimated from coords spacing
epsg_in – int CRS of data (if transformation is required)
epsg_out – int CRS of grid (if transformation is required)
kwargs – keyword arguments keyword arguments to pass to griddata from scipy.interpolate.griddata

Returns:

grid –

array shape (nrows, ncols): rectangular section of data bounded by extent

pycurious.mapping.ungrid(grid, extent, coordinates, **kwargs)[source]¶

Maps the value at unstructured coordinates from a regular 2D grid. Uses the scipy.ndimage.map_coordinates function with cubic interpolation

Parameters:

grid – array shape (ny,nx) regularly spaced grid
extent – tuple bounding box of grid e.g. [xmin,xmax,ymin,ymax]
coordinates – array shape (n,2) coordinates (x,y) to interpolate grid
kwargs – keyword arguments keyword arguments to pass to map_coordinates

Returns:

grid_interp –

array shape (n,): interpolated values at coordinates

pycurious.mapping.import_geotiff(file_path)[source]¶

Import a GeoTIFF to a numpy array and prints information of the Coordinate Reference System (CRS).

Requires osgeo.

Parameters:

file_path – str path to the GeoTIFF

Returns:

data –

2D array: the raster values
extenttuple: bounding box in the projection of the GeoTIFF, e.g. [xmin, xmax, ymin, ymax]

pycurious.mapping.export_geotiff(file_path, array, extent, epsg)[source]¶

Export a GeoTIFF from a numpy array projected in a predefined Coordinate Reference System (CRS).

Requires osgeo.

Parameters:

file_path – str path to write the GeoTIFF
array – 2D array array to save to GeoTiff
extent – tuple bounding box in the projection of the GeoTIFF e.g. [xmin, xmax, ymin, ymax]
epsg – int CRS of the GeoTIFF e.g. 4326 for WGS84

pycurious.mapping.export_netcdf4(file_path, array, extent)[source]¶

Export a netCDF4 file from a numpy array over a user-defined extent.

Requires netcdf4. pip install netcdf4

Parameters:

file_path – str path to write the netCDF4 file to
array – 2D array array to save to netCDF4
extent – tuple bounding box in the projection of the netCDF4 file e.g. [xmin, xmax, ymin, ymax]

pycurious.mapping.import_netcdf4(file_path)[source]¶

Import a netCDF4 file from a numpy array over a user-defined extent.

Requires netcdf4. pip install netcdf4

Parameters:

file_path – str path to the netCDF4 file to read

Returns:

array –

ndarray: numpy array containing the gridded netCDF4 data
extenttuple: bounding box in the projection of the netCDF4 file e.g. [xmin, xmax, ymin, ymax]

Note

The layout of the netCDF4 file must follow the standard naming convention such as ‘lon’, ‘lat’, ‘z’, ‘data’, etc otherwise this import function will not work as intended.

Downloads¶

This module provides some simple functions to download or install files.

It is good practise to store a checksum for each file. A checksum is a unique signature used to identify a file, and ensure that the data contained within the file is legitimate. Here is a sample workflow for downloading EMAG2 (v3) and evaluating its checksum:

resource = {
     "local_file": "../../data/EMAG2_V3_20170530.npz",
     "md5": "c0898b6a91efb3f13783873a8b67380c",
     "url": "https://zenodo.org/record/3245551/files/EMAG2_V3_20170530.npz?download=1",
     "expected_size": "500Mb",
}

download_cached_file(resource["url"], resource["local_file"], resource["md5"], resource["expected_size"])

The file will commence downloading if it does not already exist in the local directory or if the checksum (md5) does not match. If you do not know the checksum for a file, run

md5sum(filename)

to return its unique identifier.

pycurious.download.download_file(url, local_filename, expected_size='Unknown')[source]¶

Download files from a URL to a local filename.

Parameters:

url – str URL that points to the file to be downloaded
local_filename – str download content to this filename
expected_size – float / int / str, optional expected size of the file, used only for progress reporting (default=”Unknown”)

pycurious.download.md5sum(filename)[source]¶

Compute the MD5 checksum of a file.

Parameters:

filename – str path to the file to checksum

Returns:

digest –

str: hexadecimal MD5 digest of the file contents

pycurious.download.download_cached_file(location_url, local_file, expected_md5, expected_size='Unknown')[source]¶

Download files from a URL to a local filename.

Same as download_file() but checks whether the file has already been downloaded from its checksum.

Parameters:

location_url – str URL that points to the file to be downloaded
local_file – str download content to this filename
expected_md5 – str checksum belonging to that file
expected_size – float / int / str (optional) optional size of file (default=”Unknown”)

Returns:

status –

int: 2 if the cached file was reused, 1 if the file was downloaded, 0 if the download failed

pycurious.download.report_cached_file(local_file, expected_md5)[source]¶

Report whether the local file matches its checksum

Parameters:

local_file – str string pointing to the local file
expected_md5 – str checksum assigned to local_file.

Deprecated¶

tanaka1999 and ComputeTanaka implement the centroid method without uncertainties and are superseded by CurieOptimiseTanaka. They are retained for backward compatibility.

pycurious.tanaka1999(k, lnPhi, sigma_lnPhi, kmin_range=(0.05, 0.2), kmax_range=(0.05, 0.2))[source]¶

Compute weighted linear fit of Phi over spatial frequency window kmin:kmax

Parameters:

k – float / 1D-array wavenumber in rad/km
lnPhi – float / 1D array log of the radial power spectrum, expected in ln(sqrt(S)) form (as returned by radial_spectrum with power=1)
sigma_lnPhi – float / 1D array standard deviation of lnPhi
kmin_range – tuple (default:(0.05, 0.2)) minimum and maximum range of spatial frequencies to fit for the top of magnetic sources - ideally low frequency, straight line
kmax_range – tuple (default:(0.05, 0.2)) minimum and maximum range of spatial frequencies to fit for the bottom of magnetic source - ideally low frequency, straight line

Returns:

upper_source –

tuple: (Ztr,btr,dZtr) gradient, intercept, error for the top of magnetic sources
lower_sourcetuple: (Zor,bor,dZor) gradient, intercept, error for the bottom of magnetic sources

Notes

Deprecated since version 2.0: Use pycurious.optimise_tanaka.CurieOptimiseTanaka.optimise, which fits both bands with scipy.optimize.curve_fit and returns depths positive downwards with their uncertainties.

This hand-rolled weighted least squares squares an already-squared error term, so it weights by 1/sigma**4 rather than 1/sigma**2, and it subtracts ln(k) from a standard deviation. Its uncertainties are therefore not meaningful. It is retained only so existing scripts keep running.

pycurious.ComputeTanaka(zT, dzT, z0, dz0)[source]¶

Compute the Curie depth from the results of tanaka1999.

Deprecated since version 2.0: Use pycurious.optimise_tanaka.CurieOptimiseTanaka.calculate_CPD.

Parameters:

zT – float / 1D array top of the magnetic source
dzT – float / 1D array standard deviation of zT
z0 – float / 1D array centroid depth of the magnetic source
dz0 – float / 1D array standard deviation of z0

Returns:

CPD –

float / 1D array: estimated Curie point depth at bottom of magnetic source
CPD_stdevfloat / 1D array: standard deviation

Notes

The arguments interleave the depths with their standard deviations, whereas calculate_CPD groups them. Renaming a call without also reordering the arguments computes nonsense.