geosoft.gxpy.grid_fft submodule¶
Geosoft Fast Fourier Transform processes for 2D gridded data.
Classes: |
|
---|
Note that ‘wavenumber’ in this module refers to cycles/unit_distance. Multiply by 2*math.pi for the angular wavenumber.
See also
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)¶ 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. IfFalse
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.
-
du
¶ Wavenumber increment in the grid X direction in (cycles / grid distance uom)
-
dv
¶ Wavenumber increment in the grid Y direction in (cycles / grid distance uom).
-
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='')¶ 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) orTRN_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.
- 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
-
filtered_transform
¶ Folded descrete Fourier transform after filters applied.
-
log_average_spectral_density
(trn=0)¶ Log of the average spectral density of the transform.
Parameters: trn – TRN_SOURCE
(default) source data spectrum, orTRN_FILTERED
current filtered transform.New in version 9.4.
-
nu
¶ Number of discrete wavenumbers in grid X direction.
The transform is folded in the x direction, will be half the transform width + 1
-
nv
¶ Number of discrete wavenumbers in the grid Y direction.
-
radially_averaged_spectrum
(trn=0)¶ Radially averaged spectrum as a Numpy array shaped (n_wavenumbers, 5).
Parameters: trn – TRN_SOURCE
(default) return spectrum of the source data, orTRN_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)¶ 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) orTRN_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')
See also
New in version 9.4.
-
result_grid
(file_name=None, overwrite=False)¶ Produce a filter result grid.
Parameters: - file_name – result grid file, default greates a temporary grid
- overwrite –
True
to overwrite an existing grid
Returns: geosoft.gxpy.grid.Grid
instanceNew in version 9.4.
-
source_grid
¶ Source grid as a
geosoft.gxpy.grid.Grid
instance.
-
source_transform
¶ Folded descrete Fourier transform as a
geosoft.gxpy.grid.Grid
instance.
-
spectrum_grid
(trn=0, file_name=None, overwrite=False)¶ Return the 2D log(power) amplitude as a grid in wavenumber domain (u, v).
Amplitude = log(real**2 + imaginary**2)
Parameters: - trn –
TRN_SOURCE
source spectrum (default) orTRN_FILTERED
filtered spectrum - file_name – name for the grid file, default is a temporary grid.
- overwrite –
True
to overwrite existing grid.
Returns: geosoft.gxpy.grid.Grid
instanceNew in version 9.4.
- trn –
-
tr_row_from_uv
(i)¶ Returns transform row index from (u, v) space row index.
New in version 9.4.
-
u0
¶ First u (X-direction) wavenumber, always 0.
-
uv_row_from_tr
(i)¶ Returns (u, v) space row index of a transform row.
New in version 9.4.
-
v0
¶ First v (Y-direction) wavenumber.
-
write_uv_row
(r, i, row, trn=0)¶ Write a row (constant wavenumber v) to the (u, v) transform.
Parameters: - r – reals as a numpy array length half the width of the transform (as returned from
read_row
). - i – imaginary as a numpy array, matches r.
- row – row number in (u, v) space, row 0 is minimum v
- trn –
TRN_SOURCE
from the source transform (default) orTRN_FILTERED
See also
New in version 9.4.
- r – reals as a numpy array length half the width of the transform (as returned from
- grid – grid file name or a
-
exception
GridFFTException
(message)¶ Bases:
geosoft.GXRuntimeError
Exceptions from this module.