Total gradient amplitude transformation

[leomiquelutti]

Description of the desired feature:

A nice transformation to add to harmonica is the total gradient amplitude (aka analytic signal), defined in Roest et al (1992) as
$$|A(x, y)| = \sqrt{\left( \frac{\partial M}{\partial x} \right)^2 + \left( \frac{\partial M}{\partial y} \right)^2 + \left( \frac{\partial M}{\partial z} \right)^2}$$

Roest, W. & Verhoef, Joyce & Pilkington, Mark. (1992). Magnetic interpretation using 3-D analytic signal. GEOPHYSICS. 57. 116-125. 10.1190/1.1443174.

Are you willing to help implement and maintain this feature?

I will try to implement it but no promises about maintenance 😅

[leomiquelutti]

Since #441 is not ready, the fft seems to be the quickest option to do so. To implement it in a way that only one fft is performed onto the grid, one option is to also implement a total_gradient_amplitude_kernel in _filters.py, following the standard of the others transformations.

It should look like this:

def total_gradient_amplitude(grid):
    return apply_filter(grid, total_gradient_amplitude_kernel)

And its kernel like this:

def total_gradient_amplitude_kernel(fft_grid):
    # do something
    return apply_filter(da_filter)

If this seems ok, I need help with the da_filter calculation. The formula above is in the space domain, but how do I do it in the wavenumber domain?

[leomiquelutti]

So far I implemented this way:

def total_gradient_amplitude_kernel(fft_grid):
    
    # Catch the dims of the Fourier transformed grid
    dims = fft_grid.dims    
    # Grab the coordinates of the Fourier transformed grid
    freq_easting = fft_grid.coords[dims[1]]
    freq_northing = fft_grid.coords[dims[0]]    
    # Convert frequencies to wavenumbers
    k_easting = 2 * np.pi * freq_easting
    k_northing = 2 * np.pi * freq_northing    
    # Compute the filter for derivatives in frequency domain
    da_filter_up = np.sqrt(k_easting**2 + k_northing**2)
    da_filter_easting = (k_easting * 1j)
    da_filter_northing = (k_northing * 1j)

    return np.sqrt(da_filter_up**2 + da_filter_easting**2 + da_filter_northing**2)

Although this can generate data and hence I can create a map of it, I don't think this return is correct. The map looks weird also.

[santisoler]

Hi @leomiquelutti! Thanks for opening this issue and also to start designing an implementation to tackle it.

The kernel functions we have in harmonica/_filters are functions that return the kernel for a given FFT filter as xarray.DataArray objects. A filter $g(\mathbf{k})$ is evaluated on the wavenumber vectors $\mathbf{k}$ and returned.

These are useful when a given filter can be applied in the frequency domain and can be written as:

$$ \mathcal{F}[G(\mathbf{T})] = \mathcal{F}[\mathbf{T}] \ g(\mathbf{k}) $$

where $T$ is the 2D grid in the space domain, $\mathcal{F}$ is the Fourier transform, and $G(\mathbf{T})$ is the transformation function in the space domain (upward continuation, spatial derivative, etc).

So the transformations that make use of these filters carry out this by:

1. Compute the FFT on the data to get to the frequency domain.
2. Evaluate the filter and multiply it by the previous FFT.
3. Apply the inverse FFT to get back to the space domain.

In the particular case of the analytic signal, I'm not sure you can express it as a single FFT filter $g(\mathcal{k})$. So I don't think you can create a filter function for it. Instead, you could use the derivative_easting, derivative_northing, derivative_upward functions to compute the first derivatives along each direction, and then get the magnitude of the gradient (as the equation you pasted does).

Something like:

def total_amplitude_gradient(grid):
    gradient = (
        derivative_easting(grid, order=1),
        derivative_northing(grid, order=1),
        derivative_upward(grid, order=1),
    )
    tga = np.sqrt(np.sum([derivative**2 for derivative in gradient]))
    return tga

Hope it helps!

[leouieda]

I've just realized I forgot about the squares in the sum. I still think it should be possible to write a filter for the TGA but the maths will be harder than I had expected. Doing as @santisoler suggests could be good enough for now and we can then worry about optimizing it once we have the function.

[leomiquelutti]

I confess that I tried to solve the Fourier Transform of the TGA, but SymPy only performs 1D FT (as far as I've dug it).

Ok, I'll try going standard on the TGA though :)

[leomiquelutti]

So far I:
[x] implemented the tga function
[x] created the tga.py example
[x] added the section Total gradient amplitude (aka Analytic Signal) to transformations.rst
[x] added one test to test_transformations.py

Am I forgetting something before the PR?

[leouieda]

@leomiquelutti please open the pull request as soon as you can so the tests run and we can review. You don't have to wait until things are done and it's best to open the PR as soon as you start working on something.