diff --git a/doc/api/index.rst b/doc/api/index.rst index 75d65ce9a..9ecd77a0c 100644 --- a/doc/api/index.rst +++ b/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 diff --git a/harmonica/__init__.py b/harmonica/__init__.py index dc96420d8..e71435769 100644 --- a/harmonica/__init__.py +++ b/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, diff --git a/harmonica/synthetic/__init__.py b/harmonica/synthetic/__init__.py new file mode 100644 index 000000000..2c5b3423b --- /dev/null +++ b/harmonica/synthetic/__init__.py @@ -0,0 +1,2 @@ +# pylint: disable=missing-docstring +from .surveys import ground_survey, airborne_survey diff --git a/harmonica/synthetic/surveys.py b/harmonica/synthetic/surveys.py new file mode 100644 index 000000000..b1f6e27d6 --- /dev/null +++ b/harmonica/synthetic/surveys.py @@ -0,0 +1,161 @@ +""" +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, data_region=(-5.0, -4.0, 56.0, 56.5)): + """ + Create measurement locations for a synthetic airborne survey + + The observation points are sampled from the Great Britain total-field + magnetic anomaly dataset (see + :func:`harmonica.datasets.fetch_britain_magnetic`). A portion of the + original survey is cut (*data_region*) and the coordinates may be scaled to + the given *region*. + + Parameters + ---------- + region : tuple or list (optional) + Survey horizontal coordinates will be scaled to span this area. The + boundaries must be passed in the following order: (``east``, ``west``, + ``south``, ``north``, ...), defined on a geodetic coordinate system and + in degrees. Only the 4 horizontal boundaries are used. Subsequent + boundaries will be ignored. If ``None``, the survey points won't be + scaled. Default ``None``. + data_region : tuple or list (optional) + Subsection of the original Great Britain magnetic dataset that 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. The :class:`pandas.DataFrame` will have + the following columns: ``longitude``, ``latitude``, ``height``. + 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 data_region + if region is not None: + check_region(region[:4]) + check_region(data_region) + # 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 data_region and scale it to the passed region + survey = _cut_and_scale(survey, region, data_region) + return survey + + +def ground_survey(region=None, data_region=(13.60, 20.30, -24.20, -17.5)): + """ + Create measurement locations for 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``. + data_region : 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. The :class:`pandas.DataFrame` will have + the following columns: ``longitude``, ``latitude``, ``height``. + Longitudes and latitudes are in degrees, and heights in meters. + + See also + -------- + datasets.fetch_south_africa_gravity + """ + # Sanity checks for region and data_region + if region is not None: + check_region(region[:4]) + check_region(data_region) + # 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 data_region and scale it to the passed region + survey = _cut_and_scale(survey, region, data_region) + return survey + + +def _cut_and_scale(survey, region, data_region): + """ + Cut a subsection from the original survey and scale it to the given 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 to 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. + data_region : 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 data_region + inside_points = inside((survey.longitude, survey.latitude), data_region) + 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 diff --git a/harmonica/tests/test_synthetic_surveys.py b/harmonica/tests/test_synthetic_surveys.py new file mode 100644 index 000000000..551d36db3 --- /dev/null +++ b/harmonica/tests/test_synthetic_surveys.py @@ -0,0 +1,103 @@ +""" +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 + """ + # Expected region for the default data_region + 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) + + +def test_scale_ground_survey(): + """ + Test if the synthetic ground survey returns the expected survey after scaled + """ + 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 + """ + # Expected region for the default data_region + 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) + + +def test_scale_airborne_survey(): + """ + Test if the synthetic airborne survey returns the expected survey after scaled + """ + 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_data_region_ground_survey(): + """ + Test if ground survey is changed against a different data_region + """ + data_region = (10, 30, -30, -12) # a bigger data_region 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, data_region=data_region) + 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_data_region_airborne_survey(): + """ + Test if a different cut_region produces a different airborne survey + """ + data_region = (-7, -2, 53, 58) # a bigger data_region 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, data_region=data_region) + 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