Add ground and airborne synthetic surveys

doc/api/index.rst

@@ -70,6 +70,14 @@ Input and Output
 
  load_icgem_gdf
 
+Synthetic models and surveys
+----------------------------
+.. autosummary::
+ :toctree: generated/
+
+ synthetic.airborne_survey
+ synthetic.ground_survey
+
 
 .. automodule:: harmonica.datasets
 

harmonica/__init__.py

@@ -2,6 +2,7 @@
 # Import functions/classes to make the public API
 from . import version
 from . import datasets
+from . import synthetic
 from .ellipsoid import (
  set_ellipsoid,
  get_ellipsoid,

harmonica/synthetic/__init__.py

@@ -0,0 +1,2 @@
+# pylint: disable=missing-docstring
+from .surveys import ground_survey, airborne_survey

harmonica/synthetic/surveys.py

@@ -0,0 +1,157 @@
+"""
+Create synthetic surveys for gravity and magnetic observations
+"""
+from verde import get_region, inside
+from verde.coordinates import check_region
+
+from ..datasets import fetch_britain_magnetic, fetch_south_africa_gravity
+
+
+def airborne_survey(region=None, subsection=(-5.0, -4.0, 56.0, 56.5)):
+ """
+ Create a synthetic ground survey
+
+ The observation points are sampled from the Great Britain total-field magnetic
+ anomaly dataset (see :func:`harmonica.datasets.fetch_britain_magnetic`).
+ Only a portion of the original survey is sampled and its region may be
+ scaled to the passed ``region``.
+
+ Parameters
+ ----------
+ region : tuple or list (optional)
+ Region at which the survey points coordinates will be scaled.
+ The boundaries must be passed in the following order:
+ (``east``, ``west``, ``south``, ``north``, ...), defined on a geodetic
+ coordinate system and in degrees.
+ All subsequent boundaries will be ignored.
+ If ``None``, the survey points won't be scaled. Default ``None``.
+ subsection : tuple or list (optional)
+ Region where the original Great Britain magnetic dataset will be sampled.
+ The boundaries must be passed in the following order:
+ (``east``, ``west``, ``south``, ``north``, ...), defined on a geodetic
+ coordinate system and in degrees. All subsequent boundaries will be ignored.
+
+ Returns
+ -------
+ survey : :class:`pandas.DataFrame`
+ Dataframe containing the coordinates of the observation points on a geodetic
+ coordinate system. Longitudes and latitudes are in degrees, and heights in
+ meters.
+
+ See also
+ --------
+ datasets.fetch_britain_magnetic:
+ Fetch total-field magnetic anomaly data of Great Britain.
+ """
+ # Sanity checks for region and subsection
+ if region is not None:
+ check_region(region[:4])
+ check_region(subsection)
+ # Fetch airborne magnetic survey from Great Britain
+ survey = fetch_britain_magnetic()
+ # Rename the "elevation" column to "height" and
+ # keep only the longitude, latitude and height
+ survey = survey.rename(columns={"altitude_m": "height"}).filter(
+ ["longitude", "latitude", "height"]
+ )
+ # Cut the survey into the subsection and scale it to the passed region
+ survey = _cut_and_scale(survey, region, subsection)
+ return survey
+
+
+def ground_survey(region=None, subsection=(13.60, 20.30, -24.20, -17.5)):
+ """
+ Create a synthetic ground survey
+
+ The observation points are sampled from the South Africa gravity dataset
+ (see :func:`harmonica.datasets.fetch_south_africa_gravity`).
+ Only a portion of the original survey is sampled and its region may be scaled to the
+ passed ``region``.
+
+ Parameters
+ ----------
+ region : tuple or list (optional)
+ Region at which the survey points coordinates will be scaled.
+ The boundaries must be passed in the following order:
+ (``east``, ``west``, ``south``, ``north``, ...), defined on a geodetic
+ coordinate system and in degrees.
+ All subsequent boundaries will be ignored.
+ If ``None``, the survey points won't be scaled. Default ``None``.
+ subsection : tuple or list (optional)
+ Region where the original Great Britain magnetic dataset will be sampled.
+ The boundaries must be passed in the following order:
+ (``east``, ``west``, ``south``, ``north``, ...), defined on a geodetic
+ coordinate system and in degrees. All subsequent boundaries will be ignored.
+
+ Returns
+ -------
+ survey : :class:`pandas.DataFrame`
+ Dataframe containing the coordinates of the observation points on a geodetic
+ coordinate system. Longitudes and latitudes are in degrees, and heights in
+ meters.
+
+ See also
+ --------
+ datasets.fetch_south_africa_gravity: Fetch gravity station data from South Africa.
+ """
+ # Sanity checks for region and subsection
+ if region is not None:
+ check_region(region[:4])
+ check_region(subsection)
+ # Fetch ground gravity survey from South Africa
+ survey = fetch_south_africa_gravity()
+ # Rename the "elevation" column to "height" and
+ # keep only the longitude, latitude and height
+ survey = survey.rename(columns={"elevation": "height"}).filter(
+ ["longitude", "latitude", "height"]
+ )
+ # Cut the survey into the subsection and scale it to the passed region
+ survey = _cut_and_scale(survey, region, subsection)
+ return survey
+
+
+def _cut_and_scale(survey, region, subsection):
+ """
+ Cut the original survey to the subsection and scale it to the passed region
+
+ Parameters
+ ----------
+ survey : :class:`pandas.DataFrame`
+ Original survey as a :class:`pandas.DataFrame` containing the following columns:
+ ``longitude``, ``latitude`` and ``height``.
+ region : tuple or list (optional)
+ Region at which the survey points coordinates will be scaled.
+ The boundaries must be passed in the following order:
+ (``east``, ``west``, ``south``, ``north``, ...), defined on a geodetic
+ coordinate system and in degrees.
+ All subsequent boundaries will be ignored.
+ If ``None``, the survey points won't be scaled.
+ subsection : tuple or list (optional)
+ Region where the original Great Britain magnetic dataset will be sampled.
+ The boundaries must be passed in the following order:
+ (``east``, ``west``, ``south``, ``north``, ...), defined on a geodetic
+ coordinate system and in degrees. All subsequent boundaries will be ignored.
+
+ Returns
+ -------
+ survey : :class:`pandas.DataFrame`
+ Dataframe containing the coordinates of the observation points on a geodetic
+ coordinate system. Longitudes and latitudes are in degrees, and heights in
+ meters.
+ """
+ # Cut the data into the subsection
+ inside_points = inside((survey.longitude, survey.latitude), subsection)
+ survey = survey[inside_points].copy()
+ # Scale survey coordinates to the passed region
+ if region is not None:
+ w, e, s, n = region[:4]
+ longitude_min, longitude_max, latitude_min, latitude_max = get_region(
+ (survey.longitude, survey.latitude)
+ )
+ survey["longitude"] = (e - w) / (longitude_max - longitude_min) * (
+ survey.longitude - longitude_min
+ ) + w
+ survey["latitude"] = (n - s) / (latitude_max - latitude_min) * (
+ survey.latitude - latitude_min
+ ) + s
+ return survey

harmonica/tests/test_synthetic_surveys.py

@@ -0,0 +1,94 @@
+"""
+Test functions to create synthetic surveys
+"""
+import numpy.testing as npt
+from ..synthetic import airborne_survey, ground_survey
+
+
+def test_ground_survey():
+ """
+ Test if the sythetic ground survey returns the expected survey
+ """
+ # Ground survey without scaling
+ subsection = (13.60, 20.30, -24.20, -17.5) # this is the default subsection
+ expected_region = (13.60833, 20.28333, -24.2, -17.50333)
+ survey = ground_survey()
+ assert set(survey.columns) == set(["longitude", "latitude", "height"])
+ assert survey.longitude.size == 963
+ npt.assert_allclose(survey.longitude.min(), expected_region[0])
+ npt.assert_allclose(survey.longitude.max(), expected_region[1])
+ npt.assert_allclose(survey.latitude.min(), expected_region[2])
+ npt.assert_allclose(survey.latitude.max(), expected_region[3])
+ npt.assert_allclose(survey.height.min(), 0.0)
+ npt.assert_allclose(survey.height.max(), 2052.2)
+ # Ground survey with scaling
+ region = (-10.1, 9.7, -20.3, -10.5) # a random region to scale the survey
+ survey = ground_survey(region=region)
+ assert set(survey.columns) == set(["longitude", "latitude", "height"])
+ assert survey.longitude.size == 963
+ npt.assert_allclose(survey.longitude.min(), region[0])
+ npt.assert_allclose(survey.longitude.max(), region[1])
+ npt.assert_allclose(survey.latitude.min(), region[2])
+ npt.assert_allclose(survey.latitude.max(), region[3])
+ npt.assert_allclose(survey.height.min(), 0.0)
+ npt.assert_allclose(survey.height.max(), 2052.2)
+
+
+def test_airborne_survey():
+ """
+ Test if the synthetic airborne survey returns the expected survey
+ """
+ subsection = (-5.0, -4.0, 56.0, 56.5) # this is the default subsection
+ expected_region = (-4.99975, -4.00003, 56.00011, 56.49997)
+ survey = airborne_survey()
+ assert set(survey.columns) == set(["longitude", "latitude", "height"])
+ assert survey.longitude.size == 5673
+ npt.assert_allclose(survey.longitude.min(), expected_region[0])
+ npt.assert_allclose(survey.longitude.max(), expected_region[1])
+ npt.assert_allclose(survey.latitude.min(), expected_region[2])
+ npt.assert_allclose(survey.latitude.max(), expected_region[3])
+ npt.assert_allclose(survey.height.min(), 359.0)
+ npt.assert_allclose(survey.height.max(), 1255.0)
+ # Ground survey with scaling
+ region = (-10.1, 9.7, -20.3, -10.5) # a random region to scale the survey
+ survey = airborne_survey(region=region)
+ assert set(survey.columns) == set(["longitude", "latitude", "height"])
+ assert survey.longitude.size == 5673
+ npt.assert_allclose(survey.longitude.min(), region[0])
+ npt.assert_allclose(survey.longitude.max(), region[1])
+ npt.assert_allclose(survey.latitude.min(), region[2])
+ npt.assert_allclose(survey.latitude.max(), region[3])
+ npt.assert_allclose(survey.height.min(), 359.0)
+ npt.assert_allclose(survey.height.max(), 1255.0)
+
+
+def test_subsection_ground_survey():
+ """
+ Test if ground survey is changed against a different subsection
+ """
+ subsection = (10, 30, -30, -12) # a bigger subsection than the default one
+ region = (-10.1, 9.7, -20.3, -10.5) # a random region to scale the survey
+ survey = ground_survey(region=region, subsection=subsection)
+ assert survey.longitude.size > 963
+ npt.assert_allclose(survey.longitude.min(), region[0])
+ npt.assert_allclose(survey.longitude.max(), region[1])
+ npt.assert_allclose(survey.latitude.min(), region[2])
+ npt.assert_allclose(survey.latitude.max(), region[3])
+ assert survey.height.min() <= 0.0
+ assert survey.height.max() >= 2052.2
+
+
+def test_subsection_airborne_survey():
+ """
+ Test if a different cut_region produces a different airborne survey
+ """
+ subsection = (-7, -2, 53, 58) # a bigger subsection than the default one
+ region = (-10.1, 9.7, -20.3, -10.5) # a random region to scale the survey
+ survey = airborne_survey(region=region, subsection=subsection)
+ assert survey.longitude.size > 5673
+ npt.assert_allclose(survey.longitude.min(), region[0])
+ npt.assert_allclose(survey.longitude.max(), region[1])
+ npt.assert_allclose(survey.latitude.min(), region[2])
+ npt.assert_allclose(survey.latitude.max(), region[3])
+ assert survey.height.min() <= 359.0
+ assert survey.height.max() >= 1255.0