Add theta_map function

Implements the theta map filter.
The function computes the arctangent of the horizontal gradient magnitude over
the total gradient amplitude. Includes a detailed docstring with formula and references.

Author: ErosKerouak

harmonica/_transformations.py

@@ -507,6 +507,72 @@ def tilt_angle(grid):
     return tilt
 
 
+def theta_map(grid):
+    r"""
+    Calculate the theta map of a potential field grid
+
+    Compute the theta map of a regularly gridded potential field
+    :math:`M`, used to enhance edges in gravity and magnetic data.
+    The horizontal and vertical derivatives are calculated from the
+    input grid, and the theta angle is defined as the arctangent
+    of the ratio between the total horizontal derivative and the
+    total gradient amplitude.
+
+    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.
+
+    Returns
+    -------
+    theta_grid : :class:`xarray.DataArray`
+       A :class:`xarray.DataArray` with the calculated theta angle
+        in radians.
+
+    Notes
+    -----
+    The theta angle is calculated as:
+
+    .. math::
+
+        \theta(M) = \tan^{-1} \left(
+            \frac{
+                \sqrt{
+                    \left( \frac{\partial M}{\partial x} \right)^2 +
+                    \left( \frac{\partial M}{\partial y} \right)^2
+                }
+            }{
+                \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
+                }
+            }
+        \right)
+
+    where :math:`M` is the regularly gridded potential field.
+
+    References
+    ----------
+    [Wijns et al., 2005]_
+    """
+    # 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 theta map
+    total_gradient = np.sqrt(gradient[0] ** 2 + gradient[1] ** 2 + gradient[2] ** 2)
+    horiz_deriv = np.sqrt(gradient[0] ** 2 + gradient[1] ** 2)
+    return np.arctan2(horiz_deriv, total_gradient)
+
+
 def _get_dataarray_coordinate(grid, dimension_index):
     """
     Return the name of the easting or northing coordinate in the grid