Compute upward derivative of a grid in the frequency domain (#238)

Define a new derivative_upward function for computing the spatial upward
derivative of a 2D grid in the frequency domain. The function makes use of
xrft for handling Fourier transformations of xarray objects. Add a new
filters subpackage that includes FFT filters: functions that take grids in
frequency domain and return the desired filter also in frequency domain. Add
fft and ifft wrapper functions of the xrft.fft and xrft.ifft ones. Add
a new apply_filter function that takes a grid in the spatial domain, applies
fft, the filter and ifft and returns the filtered grid also in spatial domain.
Add tests for the new features and a gallery example for the upward derivative.
Add netcdf4 as requirement for testing.

doc/api/index.rst

@@ -17,6 +17,29 @@ Gravity Corrections
 
 For the Normal Earth correction, see package :mod:`boule`.
 
+Grid Transformations
+--------------------
+
+Apply well known transformations regular gridded potential fields data.
+
+.. autosummary::
+ :toctree: generated/
+
+ derivative_upward
+
+Frequency domain filters
+------------------------
+
+Define filters in the frequency domain.
+
+.. autosummary::
+ :toctree: generated/
+
+ filters.derivative_upward_kernel
+
+Use :func:`xrft.xrft.fft` and :func:`xrft.xrft.ifft` to apply Fast-Fourier
+Transforms and its inverse on :class:`xarray.DataArray`.
+
 Equivalent Sources
 ------------------
 

doc/conf.py

@@ -45,6 +45,7 @@
  "scipy": ("https://docs.scipy.org/doc/scipy/reference", None),
  "pandas": ("http://pandas.pydata.org/pandas-docs/stable/", None),
  "xarray": ("http://xarray.pydata.org/en/stable/", None),
+ "xrft": ("https://xrft.readthedocs.io/en/stable/", None),
  "cartopy": ("https://scitools.org.uk/cartopy/docs/latest/", None),
  "pooch": ("https://www.fatiando.org/pooch/latest/", None),
  "verde": ("https://www.fatiando.org/verde/latest/", None),

doc/install.rst

@@ -33,6 +33,7 @@ Required:
 * `scikit-learn <https://scikit-learn.org>`__
 * `pooch <http://www.fatiando.org/pooch/>`__
 * `verde <http://www.fatiando.org/verde/>`__
+* `xrft <https://xrft.readthedocs.io/>`__
 
 Optional:
 

env/requirements-tests.txt

@@ -5,3 +5,4 @@ pytest-cov
 coverage
 pyvista
 vtk>=9
+netcdf4

environment.yml

@@ -19,10 +19,12 @@ dependencies:
  - pooch>=0.7.0
  - verde>=1.5.0
  - xarray
+ - xrft
  # Optional requirements
  - pyvista
  - vtk>=9
  # Testing requirements
+ - netcdf4
  - pytest
  - pytest-cov
  - coverage

examples/upward_derivative.py

@@ -0,0 +1,95 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Upward derivative of a regular grid
+===================================
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+import pyproj
+import verde as vd
+import xrft
+
+import harmonica as hm
+
+# Fetch the sample total-field magnetic anomaly data from Great Britain
+data = hm.datasets.fetch_britain_magnetic()
+
+# Slice a smaller portion of the survey data to speed-up calculations for this
+# example
+region = [-5.5, -4.7, 57.8, 58.5]
+inside = vd.inside((data.longitude, data.latitude), region)
+data = data[inside]
+
+# Since this is a small area, we'll project our data and use Cartesian
+# coordinates
+projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean())
+easting, northing = projection(data.longitude.values, data.latitude.values)
+coordinates = (easting, northing, data.altitude_m)
+
+# Grid the scatter data using an equivalent layer
+eql = hm.EquivalentSources(depth=1000, damping=1).fit(
+ coordinates, data.total_field_anomaly_nt
+)
+
+grid_coords = vd.grid_coordinates(
+ region=vd.get_region(coordinates), spacing=500, extra_coords=1500
+)
+grid = eql.grid(grid_coords, data_names=["magnetic_anomaly"])
+grid = grid.magnetic_anomaly
+
+# Pad the grid to increase accuracy of the FFT filter
+pad_width = {
+ "easting": grid.easting.size // 3,
+ "northing": grid.northing.size // 3,
+}
+grid = grid.drop_vars("upward") # drop extra coordinates due to bug in xrft.pad
+grid_padded = xrft.pad(grid, pad_width)
+
+# Compute the upward derivative of the grid
+deriv_upward = hm.derivative_upward(grid_padded)
+
+# Unpad the derivative grid
+deriv_upward = xrft.unpad(deriv_upward, pad_width)
+
+# Show the upward derivative
+print("\nUpward derivative:\n", deriv_upward)
+
+# Plot original magnetic anomaly and the upward derivative
+fig, (ax1, ax2) = plt.subplots(
+ nrows=1, ncols=2, figsize=(9, 8), sharex=True, sharey=True
+)
+
+# Plot the magnetic anomaly grid
+grid.plot(
+ ax=ax1,
+ cmap="seismic",
+ cbar_kwargs={"label": "nT", "location": "bottom", "shrink": 0.8, "pad": 0.08},
+)
+ax1.set_title("Magnetic anomaly")
+
+# Plot the upward derivative
+scale = np.quantile(np.abs(deriv_upward), q=0.98) # scale the colorbar
+deriv_upward.plot(
+ ax=ax2,
+ vmin=-scale,
+ vmax=scale,
+ cmap="seismic",
+ cbar_kwargs={"label": "nT/m", "location": "bottom", "shrink": 0.8, "pad": 0.08},
+)
+ax2.set_title("Upward derivative")
+
+# Scale the axes
+for ax in (ax1, ax2):
+ ax.set_aspect("equal")
+
+# Set ticklabels with scientific notation
+ax1.ticklabel_format(axis="x", style="sci", scilimits=(-2, 2))
+ax1.ticklabel_format(axis="y", style="sci", scilimits=(-2, 2))
+
+plt.tight_layout()
+plt.show()

harmonica/__init__.py

@@ -18,6 +18,7 @@
 from .gravity_corrections import bouguer_correction
 from .io import load_icgem_gdf
 from .isostasy import isostasy_airy
+from .transformations import derivative_upward
 
 
 def test(doctest=True, verbose=True, coverage=False, figures=False):

harmonica/filters/__init__.py

@@ -0,0 +1,10 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Frequency domain filters meant to be applied on regular grids
+"""
+from ._filters import derivative_upward_kernel

harmonica/filters/_fft.py

@@ -0,0 +1,83 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Wrap xrft functions to compute FFTs and inverse FFTs
+"""
+
+from xrft.xrft import fft as _fft
+from xrft.xrft import ifft as _ifft
+
+
+def fft(grid, true_phase=True, true_amplitude=True, drop_bad_coords=True, **kwargs):
+ """
+ Compute Fast Fourier Transform of a 2D regular grid
+
+ Parameters
+ ----------
+ grid : :class:`xarray.DataArray`
+ A two dimensional :class:`xarray.DataArray` whose coordinates are
+ evenly spaced (regular grid). Its dimensions should be in the following
+ order: *northing*, *easting*. Its coordinates should be defined in the
+ same units.
+ true_phase : bool (optional)
+ Take the coordinates into consideration, keeping the original phase of
+ the coordinates in the spatial domain (``direct_lag``) and multiplies
+ the FFT with an exponential function corresponding to this phase.
+ Defaults to True.
+ true_amplitude : bool (optional)
+ If True, the FFT is multiplied by the spacing of the transformed
+ variables to match theoretical FT amplitude.
+ Defaults to True.
+ drop_bad_coords : bool (optional)
+ If True, only the indexes of the array will be kept before passing it
+ to :func:`xrft.fft`. Any extra coordinate should be drooped, otherwise
+ :func:`xrft.fft` raises an error after finding *bad coordinates*.
+ Defaults to True.
+
+ Returns
+ -------
+ fourier_transform : :class:`xarray.DataArray`
+ Array with the Fourier transform of the original grid.
+ """
+ if drop_bad_coords:
+ bad_coords = tuple(c for c in grid.coords if c not in grid.indexes)
+ grid = grid.drop(bad_coords)
+ return _fft(grid, true_phase=true_phase, true_amplitude=true_amplitude, **kwargs)
+
+
+def ifft(fourier_transform, true_phase=True, true_amplitude=True, **kwargs):
+ """
+ Compute Inverse Fast Fourier Transform of a 2D regular grid
+
+ Parameters
+ ----------
+ fourier_transform : :class:`xarray.DataArray`
+ Array with a regular grid defined in the frequency domain.
+ Its dimensions should be in the following order:
+ *freq_northing*, *freq_easting*.
+ true_phase : bool (optional)
+ Take the coordinates into consideration, recovering the original
+ coordinates in the spatial domain returning to the the original phase
+ (``direct_lag``), and multiplies the iFFT with an exponential function
+ corresponding to this phase.
+ Defaults to True.
+ true_amplitude : bool (optional)
+ If True, output is divided by the spacing of the transformed variables
+ to match theoretical IFT amplitude.
+ Defaults to True.
+
+ Returns
+ -------
+ grid : :class:`xarray.DataArray`
+ Array with the inverse Fourier transform of the passed grid.
+ """
+ return _ifft(
+ fourier_transform,
+ true_phase=true_phase,
+ true_amplitude=true_amplitude,
+ **kwargs
+ )

harmonica/filters/_filters.py

@@ -0,0 +1,64 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Frequency domain filters meant to be applied on regular grids
+"""
+import numpy as np
+
+
+def derivative_upward_kernel(fft_grid, order=1):
+ r"""
+ Filter for upward derivative in frequency domain
+
+ Return a :class:`xarray.DataArray` with the values of the frequency domain
+ filter for computing the upward derivative. The filter is built upon the
+ frequency coordinates of the passed ``fft_grid`` and is defined as follows:
+
+ .. math::
+
+ g(\mathbf{k}) = |\mathbf{k}| ^ n
+
+ where :math:`\mathbf{k}` is the wavenumber vector
+ (:math:`\mathbf{k} = 2\pi \mathbf{f}` where :math:`\mathbf{f}` is the
+ frequency vector) and :math:`n` is the order of the derivative.
+
+ Parameters
+ ----------
+ fft_grid : :class:`xarray.DataArray`
+ Array with the Fourier transform of the original grid.
+ Its dimensions should be in the following order:
+ *freq_northing*, *freq_easting*.
+ Use :func:`xrft.xrft.fft` and :func:`xrft.xrft.ifft` functions to
+ compute the Fourier Transform and its inverse, respectively.
+ order : int
+ The order of the derivative. Default to 1.
+
+ Returns
+ -------
+ da_filter : :class:`xarray.DataArray`
+ Array with the kernel for the upward derivative filter in frequency
+ domain.
+
+ References
+ ----------
+ [Blakely1995]_
+
+ See also
+ --------
+ harmonica.derivative_upward
+ """
+ # 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 upward derivative in frequency domain
+ da_filter = np.sqrt(k_easting**2 + k_northing**2) ** order
+ return da_filter

harmonica/filters/_utils.py

@@ -0,0 +1,64 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Utility functions for FFT filters
+"""
+import numpy as np
+
+from ._fft import fft, ifft
+
+
+def apply_filter(grid, fft_filter, **kwargs):
+ """
+ Apply a filter to a grid and return the transformed grid in spatial domain
+
+ Computes the Fourier transform of the given grid, builds the filter,
+ applies it and returns the inverse Fourier transform of the filtered grid.
+
+ Parameters
+ ----------
+ grid : :class:`xarray.DataArray`
+ A two dimensional :class:`xarray.DataArray` whose coordinates are
+ evenly spaced (regular grid). Its dimensions should be in the following
+ order: *northing*, *easting*. Its coordinates should be defined in the
+ same units.
+ fft_filter : func
+ Callable that builds the filter in the frequency domain.
+ kwargs :
+ Any additional keyword argument that should be passed to the
+ ``fft_filter``.
+
+ Returns
+ -------
+ filtered_grid : :class:`xarray.DataArray`
+ A :class:`xarray.DataArray` with the filtered version of the passed
+ ``grid``. Defined are in the spatial domain.
+ """
+ # Catch the dims of the grid
+ dims = grid.dims
+ # Check if the array has two dimensions
+ if len(dims) != 2:
+ raise ValueError(
+ f"Invalid grid with {len(dims)} dimensions. "
+ + "The passed grid must be a 2 dimensional array."
+ )
+ # Check if the grid has nans
+ if np.isnan(grid).any():
+ raise ValueError(
+ "Found nan(s) on the passed grid. "
+ + "The grid must not have missing values before computing the "
+ + "Fast Fourier Transform."
+ )
+ # Compute Fourier Transform of the grid
+ fft_grid = fft(grid)
+ # Build the filter
+ da_filter = fft_filter(fft_grid, **kwargs)
+ # Apply the filter
+ filtered_fft_grid = fft_grid * da_filter
+ # Compute inverse FFT
+ filtered_grid = ifft(filtered_fft_grid).real
+ return filtered_grid