geosoft.gxpy.grid_fft submodule¶

Geosoft Fast Fourier Transform processes for 2D gridded data.

Classes

GridFFT: Grid FFT

Note that ‘wavenumber’ in this module refers to cycles/unit_distance. Multiply by 2*math.pi for the angular wavenumber.

See also

geosoft.gxapi.GXFFT2

Note

Regression tests provide usage examples: Tests

class GridFFT(grid, buffer=2.0, buff_iterations=250, buffer_tolerance=None, trend_order=1, trend_edge=1, expand=10.0, fill_method=0, filter_length=0, amplitude_limit=0.0, edge_limit=- 1.0, edge_limit_cells=0, smooth=1, feather=False)[source]¶

Bases: object

Descrete Fourier Transform of a grid.

Parameters

grid – grid file name or a geosoft.gxpy.grid.Grid instance.
expand – minimum expansion percent to create a periodic function. The default is 10.
trend_order – trend order to remove, default is 1
fill_method – FILL_MAXIMUM_ENTROPY (default) or FILL_MINIMUM_CURVATURE. Maximum entropy prediction fills the expanded area in a way that preserves the character of the radially-averaged power spectrum so that spectral analysis based on the shape of the spectrum will be more reliable.

The following parameters only apply for maximum-entropy prediction. The defaults will be fine in all but exceptional situations where edge effects unduly distort the result.

Parameters

buffer – percentage buffer area, default 2%. The buffer expands the size of the grid footprint before filling internal space with a minimum-curvacture surface. This minimizes edge effects from high-amplitude features at the edge of the grid.
buff_iterations – maximum iterations to resolve the minimum curvature surface for the internal fill.
filter_length – maximum entropy filter length,
amplitude_limit – Amplitudes limiting, which starts at halve this setting. Default no limiting.
edge_limit – edge amplitude limiting, starting at half this value. Default no limiting.
edge_limit_cells – if edge limiting, start this many cells from the edge
smooth – True (default) to smooth the filled expanded area.
feather – True to feather expanded data to mean value at the expanded edges. If False the data will be periodic at the edges. Feathering may be useful should the prediction function introduce unreasonable edge effects.

New in version 9.4.

__init__(grid, buffer=2.0, buff_iterations=250, buffer_tolerance=None, trend_order=1, trend_edge=1, expand=10.0, fill_method=0, filter_length=0, amplitude_limit=0.0, edge_limit=- 1.0, edge_limit_cells=0, smooth=1, feather=False)[source]¶: Initialize self. See help(type(self)) for accurate signature.

property du¶: Wavenumber increment in the grid X direction in (cycles / grid distance uom)

property dv¶: Wavenumber increment in the grid Y direction in (cycles / grid distance uom).

property expanded_filled_grid¶: Expanded and filled grid as a geosoft.gxpy.grid.Grid instance.

filter(filters=None, trn=0, height='', mag_inclination='', mag_declination='', mag_strength='')[source]¶

Apply a pre-defined filter.

See filter reference: https://github.com/GeosoftInc/gxc/blob/master/reference/con_files/magmap.con

Parameters

filters – list of filters to apply. Each filter can be a string, or a tuple with the first item being the filter name followed by the filter parameters. See magmap.con referenced above for the full list of filters.
trn – TRN_SOURCE apply to the source transform (default) or TRN_FILTERED to apply to the current filtered transform.

The following parameter are the default for magnetic filed filters like pole/equator reduction and aparent susceptibility.

Parameters

height – survey ground clearance in grid distance units
mag_inclination – magnetic field inclination
mag_declination – magnetic field declination
mag_strength – total magnetic filed strength for converting magnetization to susceptibility.

Example upward continuation 500 grid distance units and a first vertical derivative:

import geosoft.gxpy.gx as gx
import geosoft.gxpy.grid_fft as gfft

gxc = gx.GXpy(log=print)

with gxfft.GridFFT('some_mag_grid_file.grd') as fft:

    # apply the filer
    fft.filter(['CNUP 500', 'DRVZ 1'])  # equlavalent to `fft.filter([('CNUP', 500), ('DRVZ', 1)])`

    # create an output grid of the upward-continued result
    fft.result_grid(file_name='upward_continued_500.grd')

New in version 9.4.

property filtered_transform¶: Folded descrete Fourier transform after filters applied.

log_average_spectral_density(trn=0)[source]¶

Log of the average spectral density of the transform.

Parameters: trn – TRN_SOURCE (default) source data spectrum, or TRN_FILTERED current filtered transform.

New in version 9.4.

property nu¶

Number of discrete wavenumbers in grid X direction.

The transform is folded in the x direction, will be half the transform width + 1

property nv¶: Number of discrete wavenumbers in the grid Y direction.

radially_averaged_spectrum(trn=0)[source]¶

Radially averaged spectrum as a Numpy array shaped (n_wavenumbers, 5).

Parameters: trn – TRN_SOURCE (default) return spectrum of the source data, or TRN_FILTERED return spectrum of the current filtered state.

Note

Numpy array shaped (n_wavenumbers, 5), where each row contains: [wavenumber, sample_count, log_power, 3-point_depth, 5-point_depth], wavenumber in cycles per 1000 * distance unit of measure (cycle/km for metres), and log_power is the natural log of the power.

Point depths are calculated by dividing the local slope(3 points and 5 points) of the log_power by (4 * pi) (see Spector and Grant, 1970).

For code clarity, the following index constants can be used to reference columns in the spectrum array:

I_WAVENUMBER

0

I_SAMPLE_COUNT

1

I_LOG_POWER

2

I_DEPTH_3

3

I_DEPTH_5

4

New in version 9.4.

read_uv_row(row, trn=0)[source]¶

Read a row (constant wavenumber v) from (u, v) transform.

Parameters

row – row number in (u, v) space, row 0 is minimum v
trn – TRN_SOURCE from the source transform (default) or TRN_FILTERED

Returns

(u_array, v, real_array, imaginary_array)

To calculate a wavenumber array: wavenumber = np.sqrt(u_array**2 + v**2).

Upward continuation example:

import geosoft.gxpy.gx as gx
import geosoft.gxpy.grid_fft as gfft

gxc = gx.GXpy(log=print)

with gxfft.GridFFT('some_mag_grid_file.grd') as fft:

    # for each row v in (u, v)
    for vrow in range(fft.nv):

        # read the row
        u, v, r, i = fft.read_uv_row(vrow)

        # angular wavenumber along the row
        wn = np.sqrt(u**2 + v**2) * 2. * math.pi

        # upward continue 500 grid distance units
        continuation_filter = np.exp(wn * -500.)
        r *= continuation_filter
        i *= continuation_filter

        # write the filtered result to the TRN_FILTERED transform
        fft.write_uv_row(r, i, vrow, trn=gxfft.TRN_FILTERED)

    # create an output grid of the upward-continued result
    fft.result_grid(file_name='upward_continued_500.grd')

I_WAVENUMBER	0
I_SAMPLE_COUNT	1
I_LOG_POWER	2
I_DEPTH_3	3
I_DEPTH_5	4