Functions to caluclate the radial average power spectrum

[santisoler]

Fixes #224

I've retaken the ideas of the closed #227 and #230.

This time I've written two new functions in gravmag.transform: a simple one that calculates the Power-Density Spectra (PDS) of regular gridded data and the other one radially averages this PDS (or any real 2D-array in the wavenumber domain).

As discussed in #230, this new functions has the following characteristics:

• The power_density_spectra() function makes use of the _fftfreq() function.
• The inputs for the radial_average() function are the wavenumbers array (kx and ky) and the pds (or any wavenumber-domain array). The outputs of the power_density_spectra() function are suitable to use as inputs of this one.
• The are no restrictions of squared x, y grids or equally spaced.
• All the calculations are made with the kx and ky distances instead of using indexes.
• The code is more numpy-ish, making the calculation even faster.
• It allows to optionally change the rings' width (making the averages more or less populated) and changing the radius of the biggest ring (in order to leave highest wavenumbers out of the averaging, or in case that our ky, ky grid isn't square to introduce those points that are left outside by default).

Now the schematic draw looks like this:

I leave an example with its outputs:

import numpy as np
import matplotlib.pyplot as plt
from fatiando.datasets import fetch_bouguer_alps_egm
from fatiando.gridder import load_surfer
from fatiando.gravmag import transform


# Fetch Alps Bouguer Anomaly Data
fetch_bouguer_alps_egm(fname="alps-bouguer.grd")
x, y, bg, shape = load_surfer("alps-bouguer.grd")

plt.contourf(y.reshape(shape), x.reshape(shape), bg.reshape(shape), 150)
plt.axes().set_aspect('equal')
plt.title("Bouguer Anomaly of the Alps Region")
plt.colorbar()
plt.savefig("bouguer.png")
plt.show()

# Calculate the Power Density Spectra
kx, ky, pds = transform.power_density_spectra(x, y, bg, shape)

plt.contourf(ky.reshape(shape), kx.reshape(shape), np.log(pds.reshape(shape)),
             150)
plt.axes().set_aspect('equal')
plt.title("Power Density Spectra of the Bouguer Anomaly")
plt.colorbar()
plt.savefig("pds.png")
plt.show()


# Radially Average the Power Density Spectra
default_width = max(kx[1, 0], ky[0, 1])
default_radius = min(kx.max(), ky.max())

# Default
k_radial, pds_radial = transform.radial_average(kx, ky, pds)
plt.plot(k_radial, np.log(pds_radial), 'o-')
plt.title("Radially Averaged PDS")
plt.savefig("default.png")
plt.show()

# Wider Rings
k_radial, pds_radial = transform.radial_average(kx, ky, pds,
                                                ring_width=4*default_width)

plt.plot(k_radial, np.log(pds_radial), 'o-')
plt.title("Radially Averaged PDS (wider rings)")
plt.savefig("wider_rings.png")
plt.show()

# Bigger Max Radius
k_radial, pds_radial = transform.radial_average(kx, ky, pds,
                                                max_radius=0.5*default_radius)

plt.plot(k_radial, np.log(pds_radial), 'o-')
plt.title("Radially Averaged PDS (max radius reduced to half)")
plt.savefig("max_radius.png")
plt.show()

Checklist:

• Make tests for new code (at least 80% coverage)
• Create/update docstrings
• Include relevant equations and citations in docstrings
• Docstrings follow the style conventions
• Code follows PEP8 style conventions
• Code and docs have been spellchecked
• Include new dependencies in doc/install.rst, requirements.txt, environment.yml, ci/requirements-conda.txt and ci/requirements-pip.txt.
• Documentation builds properly (run cd doc; make locally)
• Changelog entry (leave for last to avoid conflicts)
• Firt-time contributor? Add yourself to doc/contributors.rst (leave for last to avoid conflicts)

[leouieda]

@santis19 thank you for opening this again!

I've taken a quick look at the code and it's very nice and clean! There isn't much to add to the code. I'll leave some comments inline with some suggestions to improve the interface and docstrings. Great job!

A nice thing about the separation is that we can now test the radial_average function on simpler input. I'd suggest testing by creating some data that is made up of rings of integers. Then, if we specify the correct ring_width, we should get the actual integer values as output.

Another test could then be to create a data array that is the distance from the center of the grid. The radial average should be the distance returned to us.

These two test cases should cover most of the code and usage of this function.

I'm having a harder time thinking how we could test the power_density_spectra function. Maybe we could test it on a sin or cos oscillating in one dimension only? Then go and test it on the other dimension? The tests should also check is kx and ky are reasonable. This kind of thing is much harder to test (when you don't really know the output).

Sidetrack

Another point I'm considering is if it would be better to make a sane_fft2d function instead of a power spectra function (or maybe both). The sane_fft2d function would calculate the FFT and return it and the wavenumber arrays. It should also pad the data using the gridder padding functions by @drandykass (optionally), since this is done in most FFT applications anyway. We can offer the power spectra function function as well that takes the space-domain data and uses the sane_fft2d function, much like tga does for the derivatives.

I like this approach better because it gets rid of repeated code in the transform module and would help people develop other functions for FFT processing. But this should be kept for another PR 😄 I'll open an issue for this so that we can keep track.

End sidetrack

For this PR, I think you should focus now on implementing the tests discussed above. They should go in the fatiando/gravmag/tests/test_transform.py file.

Let me know if you need help with these. I recently updated the reduction to the pole tests and the tilt tests are brand new so you should use those as inspiration.

[leouieda]

@mtb-za and @drandykass would you care to weight in on this as well?

[leouieda]

This function is generic to any 2D array, so it might be better to rename the inputs to x, y, data instead of kx, ky, pds to indicate this. We can offer an example in the docstring to make it clear how to use this to calculate the power spectra.

[santisoler]

First I thought to make it explicitly generic to any 2D array, but I noticed that the kx and ky arrays (as outputs of _fftfreq) has special properties:

• they always include the origin (kx, ky) = (0, 0)
• they are ordered in special way (first the positive and then the negatives)
• if the original array has even amount of points per axes then the kx and ky arrays are not symmetric to the origin, instead they have one extra point in the positive semiaxe.

So I decided to stick to frequency domain (or wavenumber domain) 2D arrays only. Maybe we can generalize this function. In order to do that I think that we will have to overcome two possible problems:

• How the distances \Delta k_x and \Delta k_y are calculated.
• What would happen if the x and y arrays don't include the origin (that would lead to a nan in the radial_pds array)?

[leouieda]

OK, I hadn't thought of that. You're right, this has to be specific to something returned by FFT.

I don't think it's worth it trying to make this overly general. It will be a lot of work and there are many subtle ways that bugs can hide in the code. It's better to make it clear that this function only operates on Fourier transforms.

Since this is the case, then it might be better to rename the function to radial_average_spectrum to make it clear that it won't work for general arrays. However, it will work for other quantities in the frequency domain, not only the power spectrum. So removing explicit mentions of it should still be done.

[santisoler]

I've just finished the radial_average tests: one with integers and another one with distances (as @leouieda proposed).

I have a problem with the last one: the first ring (not the one related to k=0, but the next one) has more than 0.1 relative error. I think that may be related to the fact that this ring is the least populated, so the values are not so randomly spread inside it. But for the rest of the array the test doesn't fail.

I leave the codes here. I don't know if I should start another PR for this test or just push the edited test_transform.py into my branch...please help this github newbie! haha

from __future__ import division
import numpy as np
import numpy.testing as npt
from fatiando.gravmag import transform
from fatiando import gridder


def test_radial_average_distances():
    shape = (201, 201)
    area = (-100, 100, -100, 100)
    x, y = gridder.regular(area, shape)
    x, y = x.reshape(shape), y.reshape(shape)
    x, y = np.fft.ifftshift(x), np.fft.ifftshift(y)
    distances = np.sqrt(x**2 + y**2)
    k, radial_distances = transform.radial_average(x, y, distances)
    npt.assert_allclose(k, radial_distances, rtol=0.1)  # this fails
    npt.assert_allclose(k[2:], radial_distances[2:], rtol=0.1)  # this doesn't


def test_radial_average_integers():
    x = np.arange(-10, 11, 1)
    y = np.arange(-10, 11, 1)
    y, x = np.meshgrid(y, x)
    x, y = np.fft.ifftshift(x), np.fft.ifftshift(y)
    r = np.sqrt(x**2 + y**2)
    z = np.zeros(x.shape)
    integers = np.arange(0, x.max() + 1, 1)
    for i in integers:
        if i == 0:
            inside = r <= 0.5
        else:
            inside = np.logical_and(r > i - 0.5, r <= i + 0.5)
        z[inside] = i
    k, radial_z = transform.radial_average(x, y, z)
    npt.assert_allclose(integers, radial_z, rtol=0.1)

[santisoler]

I've edited the description of the radial_average function following @leouieda suggestions.

[leouieda]

@santis19 you can push the test here as well. A PR will only be merged when tests are implemented and passing, that's a way to guarantee that all new code works properly. You'll notice that one of the reasons why "All checks have failed" is that the code test coverage decreased (meaning that you added untested code). Just avoid pushing things that do have anything to do with these two functions (typo fixes elsewhere, bugs in other functions, etc). Those should be submitted as separate PRs.

[leouieda]

Once you push the tests to this branch we can take a look at the TravisCI output and I'll comment on the code.

[leouieda]

Great work on the tests by the way!

[santisoler]

I've just pushed a few more commits: the tests have been added and I've renamed the radial_average function to radial_average_spectrum as @leouieda proposed. Now all checks have passed!

[leouieda]

@santis19 perfect! Sorry for the delay, I was away on vacation 🏖️

I've sent you a pull request to add the changelog entry. You can merge that and add yourself to doc/contributors.rst and we're done! Let me know as soon as you do this and I'll merge this right away.

Thank you so much for taking the time to implement this! I'm working on ways of making the PR processes faster and less of a burden.

[leouieda]

@santis19 great work! I'm merging this now. Sorry it took so long. I'm rethinking the contribution process to try to make things faster.

Thanks for putting time into this!