Merge branch 'main' into extend-py-3.12

.flake8

@@ -22,6 +22,8 @@ exclude =
 per-file-ignores =
     # disable unused-imports errors on __init__.py
     __init__.py: F401
+    # ignore shadowing of builtin for harmonica/_io
+    harmonica/_io/__init__.py: A005
 
 # Configure flake8-rst-docstrings
 # -------------------------------

.github/workflows/test.yml

@@ -44,15 +44,15 @@ jobs:
           - optional
         include:
           - dependencies: oldest
-            python: "3.8"
+            python: "3.9"
           - dependencies: latest
             python: "3.12"
           - dependencies: optional
             python: "3.12"
           # test on macos-13 (x86) using oldest dependencies and python 3.8
           - os: macos-13
             dependencies: oldest
-            python: "3.8"
+            python: "3.9"
         exclude:
           # don't test on macos-latest (arm64) with oldest dependencies
           - os: macos-latest

doc/compatibility.rst

@@ -64,3 +64,5 @@ following releases to ensure compatibility:
       - 0.4.0
     * - 3.7
       - 0.6.0
+    * - 3.8
+      - 0.6.0

doc/install.rst

@@ -42,7 +42,7 @@ There are different ways to install Harmonica:
 Which Python?
 -------------
 
-You'll need **Python 3.8 or greater**.
+You'll need **Python 3.9 or greater**.
 See :ref:`python-versions` if you require support for older versions.
 
 Dependencies

doc/user_guide/equivalent_sources/index.rst

@@ -72,7 +72,7 @@ Now we can initialize the :class:`harmonica.EquivalentSources` class.
 
    import harmonica as hm
 
-   equivalent_sources = hm.EquivalentSources(depth=10e3, damping=10)
+   equivalent_sources = hm.EquivalentSources(damping=10)
    equivalent_sources
 
 By default, it places the sources one beneath each data point at a relative
@@ -81,6 +81,16 @@ This *relative depth* can be set through the ``depth`` argument.
 Deepest sources generate smoother predictions (*underfitting*), while shallow
 ones tend to overfit the data.
 
+.. hint::
+
+   By default, since Harmonica v0.7.0, the sources will be located at a depth
+   below the data points estimated as 4.5 times the mean distance between
+   first neighboring sources. Alternatively, we can set a value for this depth
+   below the data points through the ``depth`` argument.
+
+   The estimated value for the depth of the sources can be explored through the
+   :attr:`harmonica.EquivalentSources.depth_` attribute.
+
 The ``damping`` parameter is used to smooth the coefficients of the sources and
 stabilize the least square problem. A higher ``damping`` will create smoother
 predictions, while a lower one could overfit the data and create artifacts.

env/requirements-style.txt

@@ -1,14 +1,14 @@
 # Requirements for checking code style
 pathspec
-black==22.10.*
-isort==5.10.*
-flake8==6.0.*
-flake8-bugbear==22.10.*
-flake8-builtins==2.0.*
-flake8-functions==0.0.7
+black==24.4.*
+isort==5.13.*
+flake8==7.0.*
+flake8-bugbear==24.4.*
+flake8-builtins==2.5.*
+flake8-functions==0.0.8
 flake8-mutable==1.2.*
 flake8-rst-docstrings==0.3.*
-flake8-simplify==0.19.*
-flake8-unused-arguments==0.0.12
-pep8-naming==0.13.*
+flake8-simplify==0.21.*
+flake8-unused-arguments==0.0.13
+pep8-naming==0.14.*
 burocrata==0.2.*

environment.yml

@@ -44,16 +44,16 @@ dependencies:
   - gdal
   - trame # needed by pyvista
   # Code style checks and autoformat
-  - black==22.10.*
-  - isort==5.10.*
-  - flake8==6.0.*
-  - flake8-bugbear==22.10.*
-  - flake8-builtins==2.0.*
-  - flake8-functions==0.0.7
+  - black==24.4.*
+  - isort==5.13.*
+  - flake8==7.0.*
+  - flake8-bugbear==24.4.*
+  - flake8-builtins==2.5.*
+  - flake8-functions==0.0.8
   - flake8-mutable==1.2.*
   - flake8-rst-docstrings==0.3.*
-  - flake8-simplify==0.19.*
-  - flake8-unused-arguments==0.0.12
-  - pep8-naming==0.13.*
+  - flake8-simplify==0.21.*
+  - flake8-unused-arguments==0.0.13
+  - pep8-naming==0.14.*
   - pip:
     - burocrata==0.2.*

harmonica/_equivalent_sources/cartesian.py

@@ -7,6 +7,8 @@
 """
 Equivalent sources for generic harmonic functions in Cartesian coordinates
 """
+from __future__ import annotations
+
 import warnings
 
 import numpy as np
@@ -64,11 +66,16 @@ class EquivalentSources(vdb.BaseGridder):
 
     The depth of the sources can be controlled by the ``depth`` argument.
     Each source is located beneath each data point or block-averaged location
-    at a depth equal to its elevation minus the value of the ``depth``
-    argument.
-    In both cases a positive value of ``depth`` locates sources _beneath_ the
-    data points or the block-averaged locations, thus a negative ``depth`` will
-    put the sources _above_ them.
+    at a depth equal to its elevation minus the value of the ``depth_``
+    attribute.
+    If ``"default"`` is passed to the ``depth`` argument, then the ``depth_``
+    attribute is set to 4.5 times the mean distance between first neighboring
+    sources.
+    If a numerical value is passed to the ``depth`` argument, then this is the
+    one used for the ``depth_`` attribute.
+    A positive value of ``depth_`` locates sources _beneath_ the data points or
+    the block-averaged locations, thus a negative ``depth_`` will put the
+    sources _above_ them.
 
     Custom source locations can be chosen by specifying the ``points``
     argument, in which case the ``block_size`` and ``depth`` arguments will be
@@ -100,13 +107,16 @@ class EquivalentSources(vdb.BaseGridder):
         If None, will place one point source below each observation point at
         a fixed relative depth below the observation point [Cooper2000]_.
         Defaults to None.
-    depth : float
+    depth : float or "default"
         Parameter used to control the depth at which the point sources will be
         located.
-        Each source is located beneath each data point (or block-averaged
-        location) at a depth equal to its elevation minus the ``depth`` value.
+        If a value is provided, each source is located beneath each data point
+        (or block-averaged location) at a depth equal to its elevation minus
+        the ``depth`` value.
+        If set to ``"default"``, the depth of the sources will be estimated as
+        4.5 times the mean distance between first neighboring sources.
         This parameter is ignored if *points* is specified.
-        Defaults to 500.
+        Defaults to ``"default"``.
     block_size: float, tuple = (s_north, s_east) or None
         Size of the blocks used on block-averaged equivalent sources.
         If a single value is passed, the blocks will have a square shape.
@@ -129,6 +139,10 @@ class EquivalentSources(vdb.BaseGridder):
         Coordinates of the equivalent point sources.
     coefs_ : array
         Estimated coefficients of every point source.
+    depth_ : float or None
+        Estimated depth of the sources calculated as 4.5 times the mean
+        distance between first neighboring sources. This attribute is set to
+        None if ``points`` is passed.
     region_ : tuple
         The boundaries (``[W, E, S, N]``) of the data used to fit the
         interpolator. Used as the default region for the
@@ -154,11 +168,16 @@ def __init__(
         self,
         damping=None,
         points=None,
-        depth=500,
+        depth: float | str = "default",
         block_size=None,
         parallel=True,
         dtype="float64",
     ):
+        if isinstance(depth, str) and depth != "default":
+            raise ValueError(
+                f"Found invalid 'depth' value equal to '{depth}'. "
+                "It should be 'default' or a numeric value."
+            )
         self.damping = damping
         self.points = points
         self.depth = depth
@@ -205,6 +224,7 @@ def fit(self, coordinates, data, weights=None):
         if self.points is None:
             self.points_ = self._build_points(coordinates)
         else:
+            self.depth_ = None  # set depth_ to None so we don't leave it unset
             self.points_ = tuple(
                 p.astype(self.dtype) for p in vdb.n_1d_arrays(self.points, 3)
             )
@@ -220,7 +240,12 @@ def _build_points(self, coordinates):
         and apply block-averaging if ``block_size`` is not None.
         The point sources will be placed beneath the (averaged) observation
         points at a depth calculated as the elevation of the data point minus
-        the ``depth``.
+        the ``depth_`` attribute.
+
+        If ``depth`` is set to ``"default"``, the ``depth_`` attribute is set
+        as 4.5 times the mean distance between first neighboring sources.
+        If ``depth`` is set to a numerical value, this is used for the
+        ``depth_`` attribute.
 
         Parameters
         ----------
@@ -238,10 +263,14 @@ def _build_points(self, coordinates):
         """
         if self.block_size is not None:
             coordinates = self._block_average_coordinates(coordinates)
+        if self.depth == "default":
+            self.depth_ = 4.5 * np.mean(vd.median_distance(coordinates, k_nearest=1))
+        else:
+            self.depth_ = self.depth
         return (
             coordinates[0],
             coordinates[1],
-            coordinates[2] - self.depth,
+            coordinates[2] - self.depth_,
         )
 
     def _block_average_coordinates(self, coordinates):
@@ -425,6 +454,7 @@ def grid(
                 f"The {args} arguments are being ignored. The 'grid' method "
                 + "will not take any keyword arguments in the next Harmonica release",
                 FutureWarning,
+                stacklevel=1,
             )
 
         # Grid data

harmonica/_equivalent_sources/gradient_boosted.py

@@ -316,7 +316,9 @@ def _create_windows(self, coordinates, shuffle=True):
             ndata = coordinates[0].size
             if ndata <= 5e3:
                 warnings.warn(
-                    f"Found {ndata} number of coordinates (<= 5e3). Only one window will be used."
+                    f"Found {ndata} number of coordinates (<= 5e3). "
+                    "Only one window will be used.",
+                    stacklevel=1,
                 )
                 source_windows_nonempty = [np.arange(self.points_[0].size)]
                 data_windows_nonempty = [np.arange(ndata)]

harmonica/_equivalent_sources/spherical.py

@@ -311,6 +311,7 @@ def grid(
                 f"The {args} arguments are being ignored. The 'grid' method "
                 + "will not take any keyword arguments in the next Harmonica release",
                 FutureWarning,
+                stacklevel=1,
             )
 
         # Grid data

harmonica/_forward/prism_gravity.py

@@ -261,7 +261,9 @@ def _check_singular_points(coordinates, prisms, field):
         return None
     if functions[field](coordinates, prisms):
         warnings.warn(
-            "Found observation point on singular point of a prism.", UserWarning
+            "Found observation point on singular point of a prism.",
+            UserWarning,
+            stacklevel=1,
         )
 
 

harmonica/_forward/prism_layer.py

@@ -420,10 +420,10 @@ def _get_nonans_mask(self, property_name=None):
             # Warn if a nan is found within the masked property
             if not mask_property[mask].all():
                 warnings.warn(
-                    'Found missing values in "{}" property '.format(property_name)
-                    + "of the prisms layer. "
-                    + "The prisms with a nan as "
-                    + '"{}" will be ignored.'.format(property_name)
+                    f"Found missing values in '{property_name}' property "
+                    + "of the prisms layer. The prisms with a nan as "
+                    + f"'{property_name}' will be ignored.",
+                    stacklevel=1,
                 )
             mask = np.logical_and(mask, mask_property)
         return mask

harmonica/_forward/tesseroid_layer.py

@@ -329,10 +329,10 @@ def _get_nonans_mask(self, property_name=None):
             # Warn if a nan is found within the masked property
             if not mask_property[mask].all():
                 warnings.warn(
-                    'Found missing values in "{}" property '.format(property_name)
-                    + "of the tesseroid layer. "
-                    + "The tesseroids with nan as "
-                    + '"{}" will be ignored.'.format(property_name)
+                    f"Found missing values in '{property_name}' property "
+                    + "of the tesseroid layer. The tesseroids with nan as "
+                    + f"'{property_name}' will be ignored.",
+                    stacklevel=1,
                 )
             mask = np.logical_and(mask, mask_property)
         return mask

harmonica/_io/icgem_gdf.py

@@ -123,7 +123,7 @@ def _read_gdf_file(fname, **kwargs):
         rawdata = np.loadtxt(gdf_file, ndmin=2, unpack=True, **kwargs)
     _check_gdf_integrity(metadata)
     # Remove column names from the metadata if they weren't read
-    if kwargs.get("usecols", None) is not None:
+    if kwargs.get("usecols") is not None:
         metadata["attributes"] = [metadata["attributes"][i] for i in kwargs["usecols"]]
     if len(metadata["attributes"]) != rawdata.shape[0]:
         raise IOError(