Pad grids automatically before applying transformations

doc/user_guide/transformations.rst

@@ -42,41 +42,6 @@ And plot it:
    :func:`verde.project_grid` function and a map projection like the ones
    available in :mod:`pyproj`.
 
-Since all the grid transformations we are going to apply are based on FFT
-methods, we usually want to pad them in order their increase the accuracy.
-We can easily do it through the :func:`xrft.pad` function.
-First we need to define how much padding we want to add along each direction.
-We will add one third of the width and height of the grid to each side:
-
-.. jupyter-execute::
-
-    pad_width = {
-        "easting": magnetic_grid.easting.size // 3,
-        "northing": magnetic_grid.northing.size // 3,
-    }
-
-And then we can pad it, but dropping the ``height`` coordinate first (this is
-needed by the :func:`xrft.pad` function):
-
-.. jupyter-execute::
-
-    import xrft
-
-    magnetic_grid_no_height = magnetic_grid.drop_vars("height")
-    magnetic_grid_padded = xrft.pad(magnetic_grid_no_height, pad_width)
-    magnetic_grid_padded
-
-.. jupyter-execute::
-
-    tmp = magnetic_grid_padded.plot(cmap="seismic", center=0, add_colorbar=False)
-    plt.gca().set_aspect("equal")
-    plt.title("Padded magnetic anomaly grid")
-    plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
-    plt.colorbar(tmp, label="nT")
-    plt.show()
-
-Now that we have the padded grid, we can apply any grid transformation.
-
 
 Upward derivative
 -----------------
@@ -88,15 +53,7 @@ magnetic anomaly grid using the :func:`harmonica.derivative_upward` function:
 
     import harmonica as hm
 
-    deriv_upward = hm.derivative_upward(magnetic_grid_padded)
-    deriv_upward
-
-This grid includes all the padding we added to the original magnetic grid, so
-we better unpad it using :func:`xrft.unpad`:
-
-.. jupyter-execute::
-
-    deriv_upward = xrft.unpad(deriv_upward, pad_width)
+    deriv_upward = hm.derivative_upward(magnetic_grid)
     deriv_upward
 
 And plot it:
@@ -160,14 +117,12 @@ frequency domain:
 
 .. jupyter-execute::
 
-    deriv_easting = hm.derivative_easting(magnetic_grid_padded, method="fft")
-    deriv_easting = xrft.unpad(deriv_easting, pad_width)
+    deriv_easting = hm.derivative_easting(magnetic_grid, method="fft")
     deriv_easting
 
 .. jupyter-execute::
 
-    deriv_northing = hm.derivative_northing(magnetic_grid_padded, method="fft")
-    deriv_northing = xrft.unpad(deriv_northing, pad_width)
+    deriv_northing = hm.derivative_northing(magnetic_grid, method="fft")
     deriv_northing
 
 .. jupyter-execute::
@@ -213,17 +168,9 @@ to upward continue it a height displacement of 500m:
 .. jupyter-execute::
 
     upward_continued = hm.upward_continuation(
-        magnetic_grid_padded, height_displacement=500
+        magnetic_grid, height_displacement=500
     )
 
-This grid includes all the padding we added to the original magnetic grid, so
-we better unpad it using :func:`xrft.unpad`:
-
-.. jupyter-execute::
-
-    upward_continued = xrft.unpad(upward_continued, pad_width)
-    upward_continued
-
 And plot it:
 
 .. jupyter-execute::
@@ -269,11 +216,8 @@ remanence), then we can apply the reduction to the pole passing only the
 .. jupyter-execute::
 
     rtp_grid = hm.reduction_to_pole(
-        magnetic_grid_padded, inclination=inclination, declination=declination
+        magnetic_grid, inclination=inclination, declination=declination
     )
-
-    # Unpad the reduced to the pole grid
-    rtp_grid = xrft.unpad(rtp_grid, pad_width)
     rtp_grid
 
 And plot it:
@@ -296,15 +240,12 @@ magnetization vector of the sources, we can specify the
     mag_inclination, mag_declination = -25, 21
 
     tmp = rtp_grid = hm.reduction_to_pole(
-        magnetic_grid_padded,
+        magnetic_grid,
         inclination=inclination,
         declination=declination,
         magnetization_inclination=mag_inclination,
         magnetization_declination=mag_declination,
     )
-
-    # Unpad the reduced to the pole grid
-    rtp_grid = xrft.unpad(rtp_grid, pad_width)
     rtp_grid
 
 .. jupyter-execute::
@@ -337,24 +278,17 @@ Let's define a cutoff wavelength of 5 kilometers:
 
     cutoff_wavelength = 5e3
 
-Then apply the two filters to our padded magnetic grid:
+Then apply the two filters to our magnetic grid:
 
 .. jupyter-execute::
 
     magnetic_low_freqs = hm.gaussian_lowpass(
-        magnetic_grid_padded, wavelength=cutoff_wavelength
+        magnetic_grid, wavelength=cutoff_wavelength
     )
     magnetic_high_freqs = hm.gaussian_highpass(
-        magnetic_grid_padded, wavelength=cutoff_wavelength
+        magnetic_grid, wavelength=cutoff_wavelength
     )
 
-And unpad them:
-
-.. jupyter-execute::
-
-    magnetic_low_freqs = xrft.unpad(magnetic_low_freqs, pad_width)
-    magnetic_high_freqs = xrft.unpad(magnetic_high_freqs, pad_width)
-
 .. jupyter-execute::
 
     magnetic_low_freqs
@@ -422,11 +356,8 @@ We can apply it through the :func:`harmonica.total_gradient_amplitude` function.
 .. jupyter-execute::
 
     tga_grid = hm.total_gradient_amplitude(
-        magnetic_grid_padded
+        magnetic_grid
     )
-
-    # Unpad the total gradient amplitude grid
-    tga_grid = xrft.unpad(tga_grid, pad_width)
     tga_grid
 
 And plot it:

harmonica/_transformations.py

@@ -54,7 +54,7 @@ def derivative_upward(grid, order=1):
     --------
     harmonica.filters.derivative_upward_kernel
     """
-    return apply_filter(grid, derivative_upward_kernel, order=order)
+    return apply_filter(grid, derivative_upward_kernel, filter_kwargs={"order": order})
 
 
 def derivative_easting(grid, order=1, method="finite-diff"):
@@ -107,7 +107,9 @@ def derivative_easting(grid, order=1, method="finite-diff"):
         for _ in range(order):
             grid = grid.differentiate(coord=coordinate)
     elif method == "fft":
-        grid = apply_filter(grid, derivative_easting_kernel, order=order)
+        grid = apply_filter(
+            grid, derivative_easting_kernel, filter_kwargs={"order": order}
+        )
     else:
         msg = (
             f"Invalid method '{method}'. Please select one from 'finite-diff' or 'fft'."
@@ -166,7 +168,9 @@ def derivative_northing(grid, order=1, method="finite-diff"):
         for _ in range(order):
             grid = grid.differentiate(coord=coordinate)
     elif method == "fft":
-        return apply_filter(grid, derivative_northing_kernel, order=order)
+        return apply_filter(
+            grid, derivative_northing_kernel, filter_kwargs={"order": order}
+        )
     else:
         msg = (
             f"Invalid method '{method}'. Please select one from 'finite-diff' or 'fft'."
@@ -209,7 +213,9 @@ def upward_continuation(grid, height_displacement):
     harmonica.filters.upward_continuation_kernel
     """
     return apply_filter(
-        grid, upward_continuation_kernel, height_displacement=height_displacement
+        grid,
+        upward_continuation_kernel,
+        filter_kwargs={"height_displacement": height_displacement},
     )
 
 
@@ -245,7 +251,12 @@ def gaussian_lowpass(grid, wavelength):
     --------
     harmonica.filters.gaussian_lowpass_kernel
     """
-    return apply_filter(grid, gaussian_lowpass_kernel, wavelength=wavelength)
+    return apply_filter(
+        grid,
+        gaussian_lowpass_kernel,
+        pad=False,
+        filter_kwargs={"wavelength": wavelength},
+    )
 
 
 def gaussian_highpass(grid, wavelength):
@@ -280,7 +291,12 @@ def gaussian_highpass(grid, wavelength):
     --------
     harmonica.filters.gaussian_highpass_kernel
     """
-    return apply_filter(grid, gaussian_highpass_kernel, wavelength=wavelength)
+    return apply_filter(
+        grid,
+        gaussian_highpass_kernel,
+        pad=False,
+        filter_kwargs={"wavelength": wavelength},
+    )
 
 
 def reduction_to_pole(
@@ -335,10 +351,12 @@ def reduction_to_pole(
     return apply_filter(
         grid,
         reduction_to_pole_kernel,
-        inclination=inclination,
-        declination=declination,
-        magnetization_inclination=magnetization_inclination,
-        magnetization_declination=magnetization_declination,
+        filter_kwargs={
+            "inclination": inclination,
+            "declination": declination,
+            "magnetization_inclination": magnetization_inclination,
+            "magnetization_declination": magnetization_declination,
+        },
     )
 
 
@@ -361,8 +379,8 @@ def total_gradient_amplitude(grid):
     Returns
     -------
     total_gradient_amplitude_grid : :class:`xarray.DataArray`
-        A :class:`xarray.DataArray` after calculating the
-        total gradient amplitude of the passed ``grid``.
+        A :class:`xarray.DataArray` after calculating the total gradient
+        amplitude of the passed ``grid``.
 
     Notes
     -----
@@ -398,10 +416,9 @@ def tilt_angle(grid):
     r"""
     Calculate the tilt angle of a potential field grid.
 
-    Compute the tilt of a regular gridded potential field
-    :math:`M`. The horizontal derivatives are calculated
-    through finite-differences while the upward derivative
-    is calculated using FFT.
+    Compute the tilt of a regular gridded potential field :math:`M`. The
+    horizontal derivatives are calculated through finite-differences while the
+    upward derivative is calculated using FFT.
 
     Parameters
     ----------
@@ -442,17 +459,12 @@ def tilt_angle(grid):
     [Blakely1995]_
     [MillerSingh1994]_
     """
-    # Run sanity checks on the grid
     grid_sanity_checks(grid)
-    # Calculate the gradients of the grid
-    gradient = (
-        derivative_easting(grid, order=1),
-        derivative_northing(grid, order=1),
-        derivative_upward(grid, order=1),
-    )
-    # Calculate and return the tilt
-    horiz_deriv = np.sqrt(gradient[0] ** 2 + gradient[1] ** 2)
-    tilt = np.arctan2(gradient[2], horiz_deriv)
+    deriv_east = derivative_easting(grid, order=1)
+    deriv_north = derivative_northing(grid, order=1)
+    deriv_up = derivative_upward(grid, order=1)
+    horiz_deriv = np.hypot(deriv_east, deriv_north)
+    tilt = np.arctan2(deriv_up, horiz_deriv)
     return tilt
 
 

harmonica/filters/_fft.py

@@ -8,11 +8,10 @@
 Wrap xrft functions to compute FFTs and inverse FFTs.
 """
 
-from xrft.xrft import fft as _fft
-from xrft.xrft import ifft as _ifft
+import xrft
 
 
-def fft(grid, true_phase=True, true_amplitude=True, drop_bad_coords=True, **kwargs):
+def fft(grid, true_phase=True, true_amplitude=True, **kwargs):
     """
     Compute Fast Fourier Transform of a 2D regular grid.
 
@@ -32,21 +31,15 @@ def fft(grid, true_phase=True, true_amplitude=True, drop_bad_coords=True, **kwar
         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)
+    return xrft.fft(
+        grid, true_phase=true_phase, true_amplitude=true_amplitude, **kwargs
+    )
 
 
 def ifft(fourier_transform, true_phase=True, true_amplitude=True, **kwargs):
@@ -75,9 +68,10 @@ def ifft(fourier_transform, true_phase=True, true_amplitude=True, **kwargs):
     grid : :class:`xarray.DataArray`
         Array with the inverse Fourier transform of the passed grid.
     """
-    return _ifft(
+    return xrft.ifft(
         fourier_transform,
         true_phase=true_phase,
         true_amplitude=true_amplitude,
+        lag=(None, None),  # Mutes an annoying FutureWarning from xrft
         **kwargs,
     )