diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index a369ce0ce..18d961e76 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -138,19 +138,8 @@ jobs: - name: List installed packages run: python -m pip freeze - - name: Copy test data to cache - run: | - echo "Copy data to " $HARMONICA_DATA_DIR/main - set -x -e - mkdir -p $HARMONICA_DATA_DIR/main - cp -r data/* $HARMONICA_DATA_DIR/main - env: - # Define directory where sample data will be copied - HARMONICA_DATA_DIR: ${{ runner.temp }}/cache/harmonica - - name: Run the tests run: | - ls $HARMONICA_DATA_DIR if [ ${{ matrix.os }} == 'ubuntu' ]; then # Set NUMBA_THREADING_LAYER to workqueue on Ubuntu to prevent # endless loop on some test functions that make use of Numba @@ -158,9 +147,6 @@ jobs: else make test fi - env: - # Define directory where sample data will be copied - HARMONICA_DATA_DIR: ${{ runner.temp }}/cache/harmonica - name: Convert coverage report to XML for codecov run: coverage xml diff --git a/Makefile b/Makefile index f797e369e..2affb97e4 100644 --- a/Makefile +++ b/Makefile @@ -3,7 +3,7 @@ PROJECT=harmonica TESTDIR=tmp-test-dir-with-unique-name PYTEST_ARGS=--cov-config=../.coveragerc --cov-report=term-missing --cov=$(PROJECT) --doctest-modules -v --pyargs NUMBATEST_ARGS=--doctest-modules -v --pyargs -m use_numba -STYLE_CHECK_FILES=$(PROJECT) examples data/examples doc/conf.py tools +STYLE_CHECK_FILES=$(PROJECT) examples doc/conf.py tools help: @echo "Commands:" diff --git a/data/README.md b/data/README.md deleted file mode 100644 index 56827e309..000000000 --- a/data/README.md +++ /dev/null @@ -1,33 +0,0 @@ -# Sample data sets - -These files are used as sample data in Harmonica: - -* `gravity-earth-0.5deg.nc.xz`: Global gravity of the Earth at 0.5 degree grid spacing. - The file contains the magnitude of the gradient of the gravity potential of the Earth - (including the centrifugal component) on a regular grid of points at a constant - geometric (ellipsoidal) height of 10 km. The grid was generated from the EIGEN-6C4 - gravity field model using the [ICGEM Calculation Service](http://icgem.gfz-potsdam.de) - and the WGS84 Reference System. The data are stored in a netCDF file and then `xz` - compressed. -* `etopo1-0.5deg.nc.xz`: Global relief (topography and bathymetry) of the Earth based on - the ETOPO1 model at 0.5 degree grid spacing. The grid was generated by the [ICGEM - Calculation Service](http://icgem.gfz-potsdam.de). The data are stored in a netCDF - file and then `xz` compressed. -* `britain-magnetic.csv.xz`: Total-field magnetic anomaly data from Great - Britain. Complete airborne survey of the entire Great Britain conducted - between 1955 and 1965. The data are made available by the British Geological - Survey (BGS) through their [geophysical data - portal](https://www.bgs.ac.uk/products/geophysics/aeromagneticRegional.html) - under an [Open Government - License](http://www.nationalarchives.gov.uk/doc/open-government-licence/version/3/). -* `south-africa-gravity.ast.xz`: - Land gravity survey performed in January 1986 within the boundaries of the Republic of - South Africa. Columns are longitude, latitude, elevation (above sea level) and gravity - (mGal). The observed gravity values are referenced to the International Gravity - Standardization Net 1971 (IGSN 71). The data was made available by the [National - Centers for Environmental Information (NCEI)](https://www.ngdc.noaa.gov/) (formerly - NGDC) and are in the [public - domain](https://www.ngdc.noaa.gov/ngdcinfo/privacy.html#copyright-notice). Original - data files can be found at: - https://www.ngdc.noaa.gov/mgg/gravity/1999/data/regional/africa/ The data are stored - in a netCDF file and then `xz` compressed. diff --git a/data/britain-magnetic.csv.xz b/data/britain-magnetic.csv.xz deleted file mode 100644 index fedfb1302..000000000 Binary files a/data/britain-magnetic.csv.xz and /dev/null differ diff --git a/data/etopo1-0.5deg.nc.xz b/data/etopo1-0.5deg.nc.xz deleted file mode 100644 index 4cd1e7bc3..000000000 Binary files a/data/etopo1-0.5deg.nc.xz and /dev/null differ diff --git a/data/examples/README.txt b/data/examples/README.txt deleted file mode 100644 index d0f2559bb..000000000 --- a/data/examples/README.txt +++ /dev/null @@ -1,32 +0,0 @@ -.. _sample_data: - -Sample Data -=========== - -Harmonica provides some sample data for testing through the -:mod:`harmonica.datasets` module. - -.. warning:: - - The :mod:`harmonica.datasets` module and every sample dataset a will be - deprecated in Harmonica v0.6.0. The examples and the user guide will - transition to using Ensaio (https://www.fatiando.org/ensaio/) instead. - - -Where is my data? ------------------ - -The sample data files are downloaded automatically by :mod:`pooch` the first -time you load them. The files are saved to the default cache location on your -operating system. The location varies depending on your system and -configuration. We provide the :func:`harmonica.datasets.locate` function if you -need to find the data storage location on your system. - -You can change the base data directory by setting the ``HARMONICA_DATA_DIR`` -environment variable to a different path. - - -Available datasets ------------------- - -These are the datasets currently available: diff --git a/data/examples/britain_magnetic.py b/data/examples/britain_magnetic.py deleted file mode 100644 index 447677bcb..000000000 --- a/data/examples/britain_magnetic.py +++ /dev/null @@ -1,71 +0,0 @@ -# Copyright (c) 2018 The Harmonica Developers. -# Distributed under the terms of the BSD 3-Clause License. -# SPDX-License-Identifier: BSD-3-Clause -# -# This code is part of the Fatiando a Terra project (https://www.fatiando.org) -# -""" -Total Field Magnetic Anomaly from Great Britain -================================================ - -These data are a complete airborne survey of the entire Great Britain -conducted between 1955 and 1965. The data are made available by the -British Geological Survey (BGS) through their `geophysical data portal -`__. - -License: `Open Government License -`__ - -The columns of the data table are longitude, latitude, total-field -magnetic anomaly (nanoTesla), observation height relative to the WGS84 datum -(in meters), survey area, and line number and line segment for each data point. - -Latitude, longitude, and elevation data converted from original OSGB36 -(epsg:27700) coordinate system to WGS84 (epsg:4326) using to_crs function in -GeoPandas. - -See the original data for more processing information. - -If the file isn't already in your data directory, it will be downloaded -automatically. -""" -import numpy as np -import pygmt -import verde as vd - -import harmonica as hm - -# Fetch the data in a pandas.DataFrame -data = hm.datasets.fetch_britain_magnetic() -print(data) - -# Get the region of the data -region = vd.get_region((data.longitude, data.latitude)) - -# Plot the observations in a Mercator map using PyGMT -fig = pygmt.Figure() - -title = "Magnetic data from Great Britain" - -# Make colormap scaled to 97 percentile of data. -maxabs = np.percentile(data.total_field_anomaly_nt, 97) -pygmt.makecpt(cmap="polar", series=(-maxabs, maxabs), background=True) - -fig.plot( - region=region, - projection="M10c", - frame=["ag", f"+t{title}"], - x=data.longitude, - y=data.latitude, - color=data.total_field_anomaly_nt, - style="c0.02c", - cmap=True, -) - -fig.colorbar( - cmap=True, position="JMR", frame=["a100f50", "x+ltotal field magnetic anomaly [nT]"] -) - -fig.coast(shorelines="1p,black") - -fig.show() diff --git a/data/examples/earth_geoid.py b/data/examples/earth_geoid.py deleted file mode 100644 index 3638ecc76..000000000 --- a/data/examples/earth_geoid.py +++ /dev/null @@ -1,44 +0,0 @@ -# Copyright (c) 2018 The Harmonica Developers. -# Distributed under the terms of the BSD 3-Clause License. -# SPDX-License-Identifier: BSD-3-Clause -# -# This code is part of the Fatiando a Terra project (https://www.fatiando.org) -# -""" -Earth Geoid -=========== - -The geoid is the equipotential surface of the Earth's gravity potential that -coincides with mean sea level. It's often represented by "geoid heights", which -indicate the height of the geoid relative to the reference ellipsoid (WGS84 in -this case). Negative values indicate that the geoid is below the ellipsoid -surface and positive values that it is above. The data are on a regular grid -with 0.5 degree spacing and was generated from the spherical harmonic model -EIGEN-6C4 [Forste_etal2014]_. -""" -import pygmt - -import harmonica as hm - -# Load the geoid grid -data = hm.datasets.fetch_geoid_earth() -print(data) - -# Make a plot using PyGMT -fig = pygmt.Figure() - -title = "Geoid heights (EIGEN-6C4)" - -fig.grdimage( - region="g", - projection="G100/0/15c", - frame=f"+t{title}", - grid=data.geoid, - cmap="vik", -) - -fig.coast(shorelines="0.5p,black", resolution="crude") - -fig.colorbar(cmap=True, frame=["a25f10", "x+lmeters"]) - -fig.show() diff --git a/data/examples/earth_gravity.py b/data/examples/earth_gravity.py deleted file mode 100644 index 25173eeda..000000000 --- a/data/examples/earth_gravity.py +++ /dev/null @@ -1,41 +0,0 @@ -# Copyright (c) 2018 The Harmonica Developers. -# Distributed under the terms of the BSD 3-Clause License. -# SPDX-License-Identifier: BSD-3-Clause -# -# This code is part of the Fatiando a Terra project (https://www.fatiando.org) -# -""" -Earth Gravity -============= - -This is the magnitude of the gravity vector of the Earth (gravitational -+ centrifugal) at 10 km height. The data is on a regular grid with 0.5 degree -spacing at 10km ellipsoidal height. It was generated from the spherical -harmonic model EIGEN-6C4 [Forste_etal2014]_. -""" -import pygmt - -import harmonica as hm - -# Load the gravity grid -data = hm.datasets.fetch_gravity_earth() -print(data) - -# Make a plot using Pygmt -fig = pygmt.Figure() - -title = "Gravity of the Earth (EIGEN-6C4)" - -fig.grdimage( - region="g", - projection="G150/0/15c", - frame=f"+t{title}", - grid=data.gravity, - cmap="viridis", -) - -fig.coast(shorelines="0.5p,black", resolution="crude") - -fig.colorbar(cmap=True, frame=["a1000f250", "x+lmGal"]) - -fig.show() diff --git a/data/examples/earth_topography.py b/data/examples/earth_topography.py deleted file mode 100644 index cbc3575f0..000000000 --- a/data/examples/earth_topography.py +++ /dev/null @@ -1,41 +0,0 @@ -# Copyright (c) 2018 The Harmonica Developers. -# Distributed under the terms of the BSD 3-Clause License. -# SPDX-License-Identifier: BSD-3-Clause -# -# This code is part of the Fatiando a Terra project (https://www.fatiando.org) -# -""" -Earth Topography -================ - -The topography and bathymetry of the Earth according to the ETOPO1 model -[AmanteEakins2009]_. The original model has 1 arc-minute grid spacing but here -we downsampled to 0.5 degree grid spacing to save space and download times. -Heights are referenced to sea level. -""" -import pygmt - -import harmonica as hm - -# Load the topography grid -data = hm.datasets.fetch_topography_earth() -print(data) - -# Make a plot using PyGMT -fig = pygmt.Figure() - -title = "Topography of the Earth (ETOPO1)" - -fig.grdimage( - region="g", - projection="G-30/0/15c", - frame=f"+t{title}", - grid=data.topography, - cmap="globe", -) - -fig.coast(shorelines="0.5p,black", resolution="crude") - -fig.colorbar(cmap=True, frame=["a2000f500", "x+lmeters"]) - -fig.show() diff --git a/data/examples/south_africa_gravity.py b/data/examples/south_africa_gravity.py deleted file mode 100644 index 7a459e19a..000000000 --- a/data/examples/south_africa_gravity.py +++ /dev/null @@ -1,55 +0,0 @@ -# Copyright (c) 2018 The Harmonica Developers. -# Distributed under the terms of the BSD 3-Clause License. -# SPDX-License-Identifier: BSD-3-Clause -# -# This code is part of the Fatiando a Terra project (https://www.fatiando.org) -# -""" -Land Gravity Data from South Africa -=================================== - -Land gravity survey performed in January 1986 within the boundaries of the -Republic of South Africa. The data was made available by the `National Centers -for Environmental Information (NCEI) `__ (formerly -NGDC) and are in the `public domain -`__. The -entire dataset is stored in a :class:`pandas.DataFrame` with columns: -longitude, latitude, elevation (above sea level) and gravity(mGal). See the -documentation for :func:`harmonica.datasets.fetch_south_africa_gravity` for -more information. -""" -import pygmt -import verde as vd - -import harmonica as hm - -# Fetch the data in a pandas.DataFrame -data = hm.datasets.fetch_south_africa_gravity() -print(data) - -# Get the region of the grid -region = vd.get_region((data.longitude.values, data.latitude.values)) - -# Make a plot using PyGMT -fig = pygmt.Figure() - -title = "Observed gravity data from South Africa" - -pygmt.makecpt(cmap="viridis", series=(data.gravity.min(), data.gravity.max())) - -fig.plot( - region=region, - projection="M15c", - frame=["ag", f"+t{title}"], - x=data.longitude, - y=data.latitude, - color=data.gravity, - style="c0.1c", - cmap=True, -) - -fig.coast(shorelines="1p,black") - -fig.colorbar(cmap=True, frame=["a200f50", "x+lobserved gravity [mGal]"], position="JMR") - -fig.show() diff --git a/data/examples/south_africa_topography.py b/data/examples/south_africa_topography.py deleted file mode 100644 index 4067d62e9..000000000 --- a/data/examples/south_africa_topography.py +++ /dev/null @@ -1,45 +0,0 @@ -# Copyright (c) 2018 The Harmonica Developers. -# Distributed under the terms of the BSD 3-Clause License. -# SPDX-License-Identifier: BSD-3-Clause -# -# This code is part of the Fatiando a Terra project (https://www.fatiando.org) -# -""" -South Africa Topography -======================= - -The topography and bathymetry of South Africa according to the ETOPO1 model -[AmanteEakins2009]_. The original model has 1 arc-minute grid spacing but here -we downsampled to 0.1 degree grid spacing to save space and download times. -Heights are referenced to sea level. -""" -import pygmt -import verde as vd - -import harmonica as hm - -# Load the topography grid -data = hm.datasets.fetch_south_africa_topography() -print(data) - -# Get the region of the grid -region = vd.get_region((data.longitude.values, data.latitude.values)) - -# Make a plot using PyGMT -fig = pygmt.Figure() - -title = "Topography of South africa (ETOPO1)" - -fig.grdimage( - region=region, - projection="M15c", - grid=data.topography, - frame=["ag", f"+t{title}"], - cmap="earth", -) - -fig.colorbar(cmap=True, frame=["a2000f500", "x+lmeters"]) - -fig.coast(shorelines="1p,black") - -fig.show() diff --git a/data/geoid-earth-0.5deg.nc.xz b/data/geoid-earth-0.5deg.nc.xz deleted file mode 100644 index 9fc2ea1b7..000000000 Binary files a/data/geoid-earth-0.5deg.nc.xz and /dev/null differ diff --git a/data/gravity-earth-0.5deg.nc.xz b/data/gravity-earth-0.5deg.nc.xz deleted file mode 100644 index 962a7c919..000000000 Binary files a/data/gravity-earth-0.5deg.nc.xz and /dev/null differ diff --git a/data/south-africa-gravity.ast.xz b/data/south-africa-gravity.ast.xz deleted file mode 100644 index a218aa090..000000000 Binary files a/data/south-africa-gravity.ast.xz and /dev/null differ diff --git a/data/south-africa-topography.nc.xz b/data/south-africa-topography.nc.xz deleted file mode 100644 index d5279a3a1..000000000 Binary files a/data/south-africa-topography.nc.xz and /dev/null differ diff --git a/doc/api/index.rst b/doc/api/index.rst index 9347a2775..6a00a12bf 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -101,50 +101,6 @@ Visualization visualization.prism_to_pyvista -Synthetic models and surveys ----------------------------- - -.. warning:: - - The :mod:`harmonica.synthetic` module will be deprecated in Harmonica - v0.6.0 - -.. autosummary:: - :toctree: generated/ - - synthetic.airborne_survey - synthetic.ground_survey - - -.. automodule:: harmonica.datasets - -.. currentmodule:: harmonica - -Datasets --------- - -.. warning:: - - The :mod:`harmonica.datasets` module and every sample dataset a will be - deprecated in Harmonica v0.6.0. The examples and the user guide will - transition to using Ensaio (https://www.fatiando.org/ensaio/) instead. - -.. warning:: - - The :mod:`harmonica.datasets` module will be deprecated in Harmonica - v0.6.0 - - -.. autosummary:: - :toctree: generated/ - - datasets.locate - datasets.fetch_gravity_earth - datasets.fetch_geoid_earth - datasets.fetch_topography_earth - datasets.fetch_britain_magnetic - datasets.fetch_south_africa_gravity - Utilities --------- diff --git a/doc/conf.py b/doc/conf.py index 5a2f6dd0b..d5ede5f39 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -91,9 +91,9 @@ # ----------------------------------------------------------------------------- sphinx_gallery_conf = { # path to your examples scripts - "examples_dirs": ["../examples", "../data/examples"], + "examples_dirs": ["../examples"], # path where to save gallery generated examples - "gallery_dirs": ["gallery", "sample_data"], + "gallery_dirs": ["gallery"], "filename_pattern": r"\.py", # Remove the "Download all examples" button from the top level gallery "download_all_examples": False, diff --git a/doc/install.rst b/doc/install.rst index afdf0fd5a..9133bbdb3 100644 --- a/doc/install.rst +++ b/doc/install.rst @@ -53,6 +53,7 @@ Optional: The examples in the :ref:`gallery` also use: * `boule `__ +* `ensaio `__ for downloading sample datasets * `pygmt `__ for plotting maps * `pyproj `__ for cartographic projections * `ensaio `__ for downloading sample datasets diff --git a/doc/references.rst b/doc/references.rst index 43de52c7d..25dc76f17 100644 --- a/doc/references.rst +++ b/doc/references.rst @@ -1,13 +1,11 @@ References ========== -.. [AmanteEakins2009] Amante, C. and B.W. Eakins, 2009. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS NGDC-24. National Geophysical Data Center, NOAA. doi:`10.7289/V5C8276M `__ .. [Balmino1973] Balmino, G., Lambeck, K., & Kaula, W. M. (1973). A spherical harmonic analysis of the Earth's topography. Journal of Geophysical Research, 78(2), 478-481. .. [BarthelmesKohler2016] Barthelmes, F. and Kohler, W. (2016), International Centre for Global Earth Models (ICGEM), in: Drewes, H., Kuglitsch, F., Adam, J. et al., The Geodesists Handbook 2016, Journal of Geodesy (2016), 90(10), pp 907-1205, doi:`10.1007/s00190-016-0948-z `__ .. [Blakely1995] Blakely, R. (1995). Potential Theory in Gravity and Magnetic Applications. Cambridge: Cambridge University Press. doi:`10.1017/CBO9780511549816 `__ .. [Cooper2000] Cooper, G.R.J. (2000), Gridding gravity data using an equivalent layer, Computers & Geosciences, Computers & Geosciences, doi:`10.1016/S0098-3004(99)00089-8 `__ .. [Dampney1969] Dampney, C. N. G. (1969). The equivalent source technique. Geophysics, 34(1), 39–53. doi:`10.1190/1.1439996 `__ -.. [Forste_etal2014] Förste, Christoph; Bruinsma, Sean.L.; Abrikosov, Oleg; Lemoine, Jean-Michel; Marty, Jean Charles; Flechtner, Frank; Balmino, G.; Barthelmes, F.; Biancale, R. (2014): EIGEN-6C4 The latest combined global gravity field model including GOCE data up to degree and order 2190 of GFZ Potsdam and GRGS Toulouse. GFZ Data Services. doi:`10.5880/icgem.2015.1 `__ .. [Fukushima2020] Fukushima, T. (2020). Speed and accuracy improvements in standard algorithm for prismatic gravitational field. Geophysical Journal International. doi:`10.1093/gji/ggaa240 `__ .. [Geosoft1999] Geosoft Incorporated. (1999). Montaj MAGMAP filtering; 2–D frequency domain of potential field data extension for Oasis Montaj v. 6.1. .. [Grombein2013] Grombein, T., Seitz, K., Heck, B. (2013), Optimized formulas for the gravitational field of a tesseroid, Journal of Geodesy. doi:`10.1007/s00190-013-0636-1 `__ diff --git a/examples/equivalent_sources/block_averaged_sources.py b/examples/equivalent_sources/block_averaged_sources.py index 9758a471e..ce5ec540d 100644 --- a/examples/equivalent_sources/block_averaged_sources.py +++ b/examples/equivalent_sources/block_averaged_sources.py @@ -40,6 +40,8 @@ The depth of the sources and which strategy to use can be set up through the ``depth`` and the ``depth_type`` parameters, respectively. """ +import ensaio +import pandas as pd import pygmt import pyproj import verde as vd @@ -47,7 +49,8 @@ import harmonica as hm # Fetch the sample total-field magnetic anomaly data from Great Britain -data = hm.datasets.fetch_britain_magnetic() +fname = ensaio.fetch_britain_magnetic(version=1) +data = pd.read_csv(fname) # Slice a smaller portion of the survey data to speed-up calculations for this # example @@ -55,13 +58,13 @@ inside = vd.inside((data.longitude, data.latitude), region) data = data[inside] print("Number of data points:", data.shape[0]) -print("Mean height of observations:", data.altitude_m.mean()) +print("Mean height of observations:", data.height_m.mean()) # Since this is a small area, we'll project our data and use Cartesian # coordinates projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean()) easting, northing = projection(data.longitude.values, data.latitude.values) -coordinates = (easting, northing, data.altitude_m) +coordinates = (easting, northing, data.height_m) xy_region = vd.get_region((easting, northing)) # Create the equivalent sources. diff --git a/examples/equivalent_sources/cartesian.py b/examples/equivalent_sources/cartesian.py index 9e9a0f846..8cbe9a6f0 100644 --- a/examples/equivalent_sources/cartesian.py +++ b/examples/equivalent_sources/cartesian.py @@ -30,6 +30,8 @@ least-squares fitting process. The main disadvantage is the increased computational load (both in terms of time and memory). """ +import ensaio +import pandas as pd import pygmt import pyproj import verde as vd @@ -37,7 +39,8 @@ import harmonica as hm # Fetch the sample total-field magnetic anomaly data from Great Britain -data = hm.datasets.fetch_britain_magnetic() +fname = ensaio.fetch_britain_magnetic(version=1) +data = pd.read_csv(fname) # Slice a smaller portion of the survey data to speed-up calculations for this # example @@ -45,13 +48,13 @@ inside = vd.inside((data.longitude, data.latitude), region) data = data[inside] print("Number of data points:", data.shape[0]) -print("Mean height of observations:", data.altitude_m.mean()) +print("Mean height of observations:", data.height_m.mean()) # Since this is a small area, we'll project our data and use Cartesian # coordinates projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean()) easting, northing = projection(data.longitude.values, data.latitude.values) -coordinates = (easting, northing, data.altitude_m) +coordinates = (easting, northing, data.height_m) xy_region = vd.get_region((easting, northing)) # Create the equivalent sources. diff --git a/examples/equivalent_sources/gradient_boosted.py b/examples/equivalent_sources/gradient_boosted.py index 62b86badc..6884715f6 100644 --- a/examples/equivalent_sources/gradient_boosted.py +++ b/examples/equivalent_sources/gradient_boosted.py @@ -29,6 +29,8 @@ """ import boule as bl +import ensaio +import pandas as pd import pygmt import pyproj import verde as vd @@ -36,7 +38,8 @@ import harmonica as hm # Fetch the sample gravity data from South Africa -data = hm.datasets.fetch_south_africa_gravity() +fname = ensaio.fetch_southern_africa_gravity(version=1) +data = pd.read_csv(fname) # Slice a smaller portion of the survey data to speed-up calculations for this # example @@ -44,19 +47,19 @@ inside = vd.inside((data.longitude, data.latitude), region) data = data[inside] print("Number of data points:", data.shape[0]) -print("Mean height of observations:", data.elevation.mean()) +print("Mean height of observations:", data.height_sea_level_m.mean()) # Since this is a small area, we'll project our data and use Cartesian # coordinates projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean()) easting, northing = projection(data.longitude.values, data.latitude.values) -coordinates = (easting, northing, data.elevation) +coordinates = (easting, northing, data.height_sea_level_m) xy_region = vd.get_region((easting, northing)) # Compute the gravity disturbance ellipsoid = bl.WGS84 -data["gravity_disturbance"] = data.gravity - ellipsoid.normal_gravity( - data.latitude, data.elevation +data["gravity_disturbance"] = data.gravity_mgal - ellipsoid.normal_gravity( + data.latitude, data.height_sea_level_m ) # Create the equivalent sources diff --git a/examples/equivalent_sources/spherical.py b/examples/equivalent_sources/spherical.py index d265d0347..1a4727a31 100644 --- a/examples/equivalent_sources/spherical.py +++ b/examples/equivalent_sources/spherical.py @@ -19,22 +19,25 @@ account the curvature of the Earth. """ import boule as bl +import ensaio import numpy as np +import pandas as pd import pygmt import verde as vd import harmonica as hm # Fetch the sample gravity data from South Africa -data = hm.datasets.fetch_south_africa_gravity() +fname = ensaio.fetch_southern_africa_gravity(version=1) +data = pd.read_csv(fname) # Downsample the data using a blocked mean to speed-up the computations # for this example. This is preferred over simply discarding points to avoid # aliasing effects. blocked_mean = vd.BlockReduce(np.mean, spacing=0.2, drop_coords=False) (longitude, latitude, elevation), gravity_data = blocked_mean.filter( - (data.longitude, data.latitude, data.elevation), - data.gravity, + (data.longitude, data.latitude, data.height_sea_level_m), + data.gravity_mgal, ) # Compute gravity disturbance by removing the gravity of normal Earth diff --git a/examples/forward/prisms_topo_gravity.py b/examples/forward/prisms_topo_gravity.py index d43e740df..44b751300 100644 --- a/examples/forward/prisms_topo_gravity.py +++ b/examples/forward/prisms_topo_gravity.py @@ -18,18 +18,29 @@ compute the gravity effect of the model on a regular grid of observation points. """ +import ensaio import pygmt import pyproj import verde as vd +import xarray as xr import harmonica as hm -# Read South Africa topography -south_africa_topo = hm.datasets.fetch_south_africa_topography() +# Read Earth's topography grid +fname = ensaio.fetch_earth_topography(version=1) +topography = xr.load_dataset(fname) + +# Crop the topography limited to South Africa +region = (12, 33, -35, -18) +region_padded = vd.pad_region(region, pad=5) # pad the original region +topography = topography.sel( + longitude=slice(*region_padded[:2]), + latitude=slice(*region_padded[2:]), +) # Project the grid -projection = pyproj.Proj(proj="merc", lat_ts=south_africa_topo.latitude.values.mean()) -south_africa_topo = vd.project_grid(south_africa_topo.topography, projection=projection) +projection = pyproj.Proj(proj="merc", lat_ts=topography.latitude.values.mean()) +south_africa_topo = vd.project_grid(topography.topography, projection=projection) # Create a 2d array with the density of the prisms Points above the geoid will # have a density of 2670 kg/m^3 Points below the geoid will have a density @@ -49,9 +60,7 @@ ) # Compute gravity field on a regular grid located at 4000m above the ellipsoid -coordinates = vd.grid_coordinates( - region=(12, 33, -35, -18), spacing=0.2, extra_coords=4000 -) +coordinates = vd.grid_coordinates(region=region, spacing=0.2, extra_coords=4000) easting, northing = projection(*coordinates[:2]) coordinates_projected = (easting, northing, coordinates[-1]) prisms_gravity = prisms.prism_layer.gravity(coordinates_projected, field="g_z") diff --git a/examples/gravity_disturbance.py b/examples/gravity_disturbance.py index b279007a9..b976f5dee 100644 --- a/examples/gravity_disturbance.py +++ b/examples/gravity_disturbance.py @@ -16,17 +16,18 @@ data. """ import boule as bl +import ensaio import pygmt - -import harmonica as hm +import xarray as xr # Load the global gravity grid -data = hm.datasets.fetch_gravity_earth() +fname = ensaio.fetch_earth_gravity(version=1) +data = xr.load_dataset(fname) print(data) # Calculate normal gravity using the WGS84 ellipsoid ellipsoid = bl.WGS84 -gamma = ellipsoid.normal_gravity(data.latitude, data.height_over_ell) +gamma = ellipsoid.normal_gravity(data.latitude, data.height) # The disturbance is the observed minus normal gravity (calculated at the # observation point) disturbance = data.gravity - gamma diff --git a/examples/gravity_disturbance_topofree.py b/examples/gravity_disturbance_topofree.py index 4d9fbfe87..737fbf6f7 100644 --- a/examples/gravity_disturbance_topofree.py +++ b/examples/gravity_disturbance_topofree.py @@ -26,24 +26,28 @@ correction. """ import boule as bl +import ensaio import pygmt import xarray as xr import harmonica as hm # Load the global gravity, topography, and geoid grids +fname_gravity = ensaio.fetch_earth_gravity(version=1) +fname_geoid = ensaio.fetch_earth_geoid(version=1) +fname_topo = ensaio.fetch_earth_topography(version=1) data = xr.merge( [ - hm.datasets.fetch_gravity_earth(), - hm.datasets.fetch_geoid_earth(), - hm.datasets.fetch_topography_earth(), + xr.load_dataarray(fname_gravity), + xr.load_dataarray(fname_geoid), + xr.load_dataarray(fname_topo), ] ) print(data) # Calculate normal gravity and the disturbance ellipsoid = bl.WGS84 -gamma = ellipsoid.normal_gravity(data.latitude, data.height_over_ell) +gamma = ellipsoid.normal_gravity(data.latitude, data.height) disturbance = data.gravity - gamma # Reference the topography to the ellipsoid diff --git a/examples/isostatic_moho_airy.py b/examples/isostatic_moho_airy.py index 6b37fd176..2ee8ba648 100644 --- a/examples/isostatic_moho_airy.py +++ b/examples/isostatic_moho_airy.py @@ -24,14 +24,17 @@ (:func:`harmonica.datasets.fetch_topography_earth`) to calculate the Airy isostatic Moho depth of Africa. """ +import ensaio import numpy as np import pygmt +import xarray as xr import harmonica as hm # Load the elevation model and cut out the portion of the data corresponding to # Africa -data = hm.datasets.fetch_topography_earth() +fname = ensaio.fetch_earth_topography(version=1) +data = xr.load_dataset(fname) region = (-20, 60, -40, 45) data_africa = data.sel(latitude=slice(*region[2:]), longitude=slice(*region[:2])) print("Topography/bathymetry grid:") diff --git a/examples/visualization/prism_layer_pyvista.py b/examples/visualization/prism_layer_pyvista.py index 740609515..5db906ab7 100644 --- a/examples/visualization/prism_layer_pyvista.py +++ b/examples/visualization/prism_layer_pyvista.py @@ -28,14 +28,17 @@ plot using ``pyvista``. """ +import ensaio import pyproj import pyvista as pv import verde as vd +import xarray as xr import harmonica as hm # Read South Africa topography -south_africa_topo = hm.datasets.fetch_south_africa_topography() +fname = ensaio.fetch_southern_africa_topography(version=1) +south_africa_topo = xr.load_dataset(fname) # Project the grid projection = pyproj.Proj(proj="merc", lat_ts=south_africa_topo.latitude.values.mean()) diff --git a/harmonica/__init__.py b/harmonica/__init__.py index d50c7a998..8bd122a3c 100644 --- a/harmonica/__init__.py +++ b/harmonica/__init__.py @@ -6,7 +6,6 @@ # # # Import functions/classes to make the public API -from . import datasets, synthetic from ._equivalent_sources.cartesian import EquivalentSources from ._equivalent_sources.gradient_boosted import EquivalentSourcesGB from ._equivalent_sources.spherical import EquivalentSourcesSph diff --git a/harmonica/datasets/__init__.py b/harmonica/datasets/__init__.py deleted file mode 100644 index 72bdc4f99..000000000 --- a/harmonica/datasets/__init__.py +++ /dev/null @@ -1,15 +0,0 @@ -# Copyright (c) 2018 The Harmonica Developers. -# Distributed under the terms of the BSD 3-Clause License. -# SPDX-License-Identifier: BSD-3-Clause -# -# This code is part of the Fatiando a Terra project (https://www.fatiando.org) -# -from .sample_data import ( - fetch_britain_magnetic, - fetch_geoid_earth, - fetch_gravity_earth, - fetch_south_africa_gravity, - fetch_south_africa_topography, - fetch_topography_earth, - locate, -) diff --git a/harmonica/datasets/registry.txt b/harmonica/datasets/registry.txt deleted file mode 100644 index 8917d45de..000000000 --- a/harmonica/datasets/registry.txt +++ /dev/null @@ -1,6 +0,0 @@ -etopo1-0.5deg.nc.xz d975b2f90111043744f70eb1382ed4c654c79065cabe71d067ea8fde7010be2b -gravity-earth-0.5deg.nc.xz 02c62a251225e2e76722a80d87c488160fbdcae2c1a8bbc2386c32e68ea8f43a -geoid-earth-0.5deg.nc.xz 1fe25866609bde37813599e34d096db0df27f5e335a51c6db0c6613eb64dc62f -britain-magnetic.csv.xz 2e6e42095dc36022c66aa3a0322ded9f0e5a5f6eabe826d5c646c70394e8872d -south-africa-gravity.ast.xz 86401cb8d2142bca340b984469f3f36dd0b40362ba028b6989ed39964ee50f6e -south-africa-topography.nc.xz b94a08b5211d318d9a84c237849d60610f52a2eb5acc74f0c755b1501d8daf53 diff --git a/harmonica/datasets/sample_data.py b/harmonica/datasets/sample_data.py deleted file mode 100644 index 8bddb45b9..000000000 --- a/harmonica/datasets/sample_data.py +++ /dev/null @@ -1,265 +0,0 @@ -# Copyright (c) 2018 The Harmonica Developers. -# Distributed under the terms of the BSD 3-Clause License. -# SPDX-License-Identifier: BSD-3-Clause -# -# This code is part of the Fatiando a Terra project (https://www.fatiando.org) -# -""" -Functions to load sample datasets used in the Harmonica docs. -""" -import warnings - -import pandas as pd -import pkg_resources -import pooch -import xarray as xr - -from .._version import __version__ as version - -REGISTRY = pooch.create( - path=pooch.os_cache("harmonica"), - base_url="https://github.com/fatiando/harmonica/raw/{version}/data/", - version=version, - version_dev="main", - env="HARMONICA_DATA_DIR", -) -with pkg_resources.resource_stream( - "harmonica.datasets", "registry.txt" -) as registry_file: - REGISTRY.load_registry(registry_file) - - -def _deprecation_warning(): - """ - Raise a FutureWarning about deprecation of synthetic module - """ - msg = ( - "The 'datasets' module will be deprecated in Harmonica v0.6.0. " - + "Harmonica will transition to using " - + "Ensaio (https://www.fatiando.org/ensaio/) instead." - ) - warnings.warn(msg, FutureWarning) - - -def locate(): - r""" - The absolute path to the sample data storage location on disk. - - This is where the data are saved on your computer. The location is - dependent on the operating system. The folder locations are defined by the - ``appdirs`` package (see the `appdirs documentation - `__). - - The location can be overwritten by the ``HARMONICA_DATA_DIR`` environment - variable to the desired destination. - - Returns - ------- - path : str - The local data storage location. - - """ - return str(REGISTRY.abspath) - - -def fetch_geoid_earth(): - """ - Fetch a global grid of the geoid height. - - The geoid height is the height of the geoid above (positive) or below - (negative) the ellipsoid (WGS84). The data are on a regular grid with 0.5 - degree spacing, which was generated from the spherical harmonic model - EIGEN-6C4 [Forste_etal2014]_ using the `ICGEM Calculation Service - `__. See the ``attrs`` attribute of the - :class:`xarray.Dataset` for information regarding the grid generation. - - If the file isn't already in your data directory, it will be downloaded - automatically. - - Returns - ------- - grid : :class:`xarray.Dataset` - The geoid grid (in meters). Coordinates are geodetic latitude and - longitude. - - """ - _deprecation_warning() - fname = REGISTRY.fetch("geoid-earth-0.5deg.nc.xz", processor=pooch.Decompress()) - data = xr.open_dataset(fname, engine="scipy") - # Capture attributes dict because it's removed after converting the data to - # float64 - attrs = data.attrs.copy() - # The data are stored as ints and data as float32 to save space on the - # data file. Cast them to float64 to avoid integer division errors. - data = data.astype("float64") - data.attrs = attrs - return data - - -def fetch_gravity_earth(): - """ - Fetch a global grid of Earth gravity. - - Gravity is the magnitude of the gravity vector of the Earth (gravitational - + centrifugal). The gravity observations are at 10 km (geometric) height - and on a regular grid with 0.5 degree spacing. The grid was generated from - the spherical harmonic model EIGEN-6C4 [Forste_etal2014]_ using the `ICGEM - Calculation Service `__. See the ``attrs`` - attribute of the :class:`xarray.Dataset` for information regarding the grid - generation. - - If the file isn't already in your data directory, it will be downloaded - automatically. - - Returns - ------- - grid : :class:`xarray.Dataset` - The gravity grid (in mGal). Includes a computation (geometric) height - grid (``height_over_ell``). Coordinates are geodetic latitude and - longitude. - - """ - _deprecation_warning() - fname = REGISTRY.fetch("gravity-earth-0.5deg.nc.xz", processor=pooch.Decompress()) - data = xr.open_dataset(fname, engine="scipy") - # Capture attributes dict because it's removed after converting the data to - # float64 - attrs = data.attrs.copy() - # The data are stored as ints and data as float32 to save space on the - # data file. Cast them to float64 to avoid integer division errors. - data = data.astype("float64") - data.attrs = attrs - return data - - -def fetch_topography_earth(): - """ - Fetch a global grid of Earth relief (topography and bathymetry). - - The grid is based on the ETOPO1 model [AmanteEakins2009]_. The original - model has 1 arc-minute grid spacing but here we downsampled to 0.5 degree - grid spacing to save space and download times. The downsampled grid was - generated from a spherical harmonic model using the `ICGEM Calculation - Service `__. See the ``attrs`` attribute of - the returned :class:`xarray.Dataset` for information regarding the grid - generation. - - ETOPO1 heights are referenced to "sea level". - - If the file isn't already in your data directory, it will be downloaded - automatically. - - Returns - ------- - grid : :class:`xarray.Dataset` - The topography grid (in meters) relative to sea level. Coordinates are - geodetic latitude and longitude. - - """ - _deprecation_warning() - fname = REGISTRY.fetch("etopo1-0.5deg.nc.xz", processor=pooch.Decompress()) - data = xr.open_dataset(fname, engine="scipy") - # Capture attributes dict because it's removed after converting the data to - # float64 - attrs = data.attrs.copy() - # The data are stored as int16 to save disk space. Cast them to floats to - # avoid integer division problems when processing. - data = data.astype("float64") - data.attrs = attrs - return data - - -def fetch_britain_magnetic(): - """ - Fetch total-field magnetic anomaly data of Great Britain. - - These data are a complete airborne survey of the entire Great Britain - conducted between 1955 and 1965. The data are made available by the - British Geological Survey (BGS) through their `geophysical data portal - `__. - - License: `Open Government License - `__ - - The columns of the data table are longitude, latitude, total-field magnetic - anomaly (nanoTesla), observation height relative to the WGS84 datum (in - meters), survey area, and line number and line segment for each data point. - - Latitude, longitude, and elevation data converted from original OSGB36 - (epsg:27700) coordinate system to WGS84 (epsg:4326) using to_crs function - in GeoPandas. - - See the original data for more processing information. - - If the file isn't already in your data directory, it will be downloaded - automatically. - - Returns - ------- - data : :class:`pandas.DataFrame` - The magnetic anomaly data. - """ - _deprecation_warning() - return pd.read_csv(REGISTRY.fetch("britain-magnetic.csv.xz"), compression="xz") - - -def fetch_south_africa_gravity(): - """ - Fetch gravity station data from South Africa - - This data base (14559 records), received in January 1986, consists in land - gravity surveys within the boundaries of the Republic of South Africa. The - survey was conducted by Dr. R.J. Kleywegt with the contribution of the - Republic of South Africa, the Department of Mineral and Energy Affairs and - the Geological Survey. The data was made available by the `National Centers - for Environmental Information (NCEI) `__ - (formerly NGDC) and are in the `public domain - `__. - Original data files can be found at: - https://www.ngdc.noaa.gov/mgg/gravity/1999/data/regional/africa/ - - Principal gravity parameters include elevation and observed gravity. The - observed gravity values are referenced to the International Gravity - Standardization Net 1971 (IGSN 71). Elevations are referenced above the sea - level. Both ``longitude`` and ``latitude`` are given in degrees, - ``elevation`` in meters and ``gravity`` in mGal. - - Returns - ------- - data : :class:`pandas.DataFrame` - The gravity data. - - """ - _deprecation_warning() - fname = REGISTRY.fetch("south-africa-gravity.ast.xz") - columns = ["latitude", "longitude", "elevation", "gravity"] - return pd.read_csv(fname, sep=r"\s+", names=columns, compression="xz") - - -def fetch_south_africa_topography(): - """ - Fetch downsampled ETOPO1 topography grid for South Africa - - The grid is based on the ETOPO1 model [AmanteEakins2009]_. The original - model has 1 arc-minute grid spacing but here we downsampled to 0.1 degree - grid spacing to save space and download times and cut it to the South - Africa region. - - ETOPO1 heights are referenced to "sea level". - - If the file isn't already in your data directory, it will be downloaded - automatically. - - Returns - ------- - grid : :class:`xarray.Dataset` - The topography grid (in meters) relative to sea level. Coordinates are - geodetic latitude and longitude. - - """ - _deprecation_warning() - fname = REGISTRY.fetch( - "south-africa-topography.nc.xz", processor=pooch.Decompress() - ) - data = xr.open_dataset(fname, engine="scipy") - return data diff --git a/harmonica/synthetic/__init__.py b/harmonica/synthetic/__init__.py deleted file mode 100644 index d253dbb99..000000000 --- a/harmonica/synthetic/__init__.py +++ /dev/null @@ -1,7 +0,0 @@ -# Copyright (c) 2018 The Harmonica Developers. -# Distributed under the terms of the BSD 3-Clause License. -# SPDX-License-Identifier: BSD-3-Clause -# -# This code is part of the Fatiando a Terra project (https://www.fatiando.org) -# -from .surveys import airborne_survey, ground_survey diff --git a/harmonica/synthetic/surveys.py b/harmonica/synthetic/surveys.py deleted file mode 100644 index 6b1a2686b..000000000 --- a/harmonica/synthetic/surveys.py +++ /dev/null @@ -1,180 +0,0 @@ -# Copyright (c) 2018 The Harmonica Developers. -# Distributed under the terms of the BSD 3-Clause License. -# SPDX-License-Identifier: BSD-3-Clause -# -# This code is part of the Fatiando a Terra project (https://www.fatiando.org) -# -""" -Create synthetic surveys for gravity and magnetic observations -""" -import warnings - -from verde import get_region, inside -from verde.coordinates import check_region - -from ..datasets import fetch_britain_magnetic, fetch_south_africa_gravity - - -def _deprecation_warning(): - """ - Raise a FutureWarning about deprecation of synthetic module - """ - warnings.warn( - "The 'synthetic' module will be deprecated in Harmonica v0.6.0.", - FutureWarning, - ) - - -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 - -------- - harmonica.datasets.fetch_britain_magnetic - """ - _deprecation_warning() - # 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 - -------- - harmonica.datasets.fetch_south_africa_gravity - """ - _deprecation_warning() - # 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_sample_data.py b/harmonica/tests/test_sample_data.py deleted file mode 100644 index fa361cec9..000000000 --- a/harmonica/tests/test_sample_data.py +++ /dev/null @@ -1,172 +0,0 @@ -# Copyright (c) 2018 The Harmonica Developers. -# Distributed under the terms of the BSD 3-Clause License. -# SPDX-License-Identifier: BSD-3-Clause -# -# This code is part of the Fatiando a Terra project (https://www.fatiando.org) -# -""" -Test the sample data loading functions. -""" -import os - -import numpy.testing as npt -import pytest - -from ..datasets.sample_data import ( - fetch_britain_magnetic, - fetch_geoid_earth, - fetch_gravity_earth, - fetch_south_africa_gravity, - fetch_south_africa_topography, - fetch_topography_earth, - locate, -) - - -@pytest.mark.parametrize( - "fetch_function", - ( - fetch_britain_magnetic, - fetch_geoid_earth, - fetch_gravity_earth, - fetch_south_africa_gravity, - fetch_south_africa_topography, - fetch_topography_earth, - ), -) -def test_deprecation_warning(fetch_function): - """ - Checks if deprecation warning is raised - """ - message = "The 'datasets' module will be deprecated in Harmonica v0.6.0.*" - with pytest.warns(FutureWarning, match=message): - fetch_function() - - -@pytest.mark.filterwarnings( - "ignore:The 'datasets' module will be deprecated:FutureWarning" -) -def test_datasets_locate(): - "Make sure the data cache location has the right package name" - # Fetch a dataset first to make sure that the cache folder exists. Since - # Pooch 1.1.1 the cache isn't created until a download is requested. - fetch_gravity_earth() - path = locate() - assert os.path.exists(path) - # This is the most we can check in a platform independent way without - # testing appdirs itself. - assert "harmonica" in path - - -@pytest.mark.filterwarnings( - "ignore:The 'datasets' module will be deprecated:FutureWarning" -) -def test_geoid_earth(): - "Sanity checks for the loaded grid" - grid = fetch_geoid_earth() - assert grid.geoid.shape == (361, 721) - npt.assert_allclose(grid.geoid.min(), -106.257344) - npt.assert_allclose(grid.geoid.max(), 84.722744) - assert grid.attrs - assert grid.attrs.get("refsysname") == "WGS84" - assert grid.attrs.get("max_used_degree") == "1277" - assert grid.attrs.get("tide_system") == "tide_free" - assert grid.attrs.get("modelname") == "EIGEN-6C4" - - -@pytest.mark.filterwarnings( - "ignore:The 'datasets' module will be deprecated:FutureWarning" -) -def test_gravity_earth(): - "Sanity checks for the loaded grid" - grid = fetch_gravity_earth() - assert grid.gravity.shape == (361, 721) - npt.assert_allclose(grid.gravity.max(), 9.8018358e05) - npt.assert_allclose(grid.gravity.min(), 9.7476403e05) - assert grid.height_over_ell.shape == (361, 721) - npt.assert_allclose(grid.height_over_ell, 10000) - assert grid.attrs - assert grid.attrs.get("refsysname") == "WGS84" - assert grid.attrs.get("max_used_degree") == "1277" - assert grid.attrs.get("tide_system") == "tide_free" - assert grid.attrs.get("modelname") == "EIGEN-6C4" - - -@pytest.mark.filterwarnings( - "ignore:The 'datasets' module will be deprecated:FutureWarning" -) -def test_topography_earth(): - "Sanity checks for the loaded grid" - grid = fetch_topography_earth() - assert grid.topography.shape == (361, 721) - npt.assert_allclose(grid.topography.max(), 5651, atol=1) - npt.assert_allclose(grid.topography.min(), -8409, atol=1) - assert grid.attrs - assert grid.attrs.get("refsysname") == "WGS84" - assert grid.attrs.get("max_used_degree") == "1277" - assert grid.attrs.get("modelname") == "etopo1-2250" - - -@pytest.mark.filterwarnings( - "ignore:The 'datasets' module will be deprecated:FutureWarning" -) -def test_britain_magnetic(): - "Sanity checks for the loaded dataset" - data = fetch_britain_magnetic() - assert data.shape == (541508, 6) - npt.assert_allclose(data.longitude.min(), -8.65338) - npt.assert_allclose(data.longitude.max(), 1.92441) - npt.assert_allclose(data.latitude.min(), 49.81407) - npt.assert_allclose(data.latitude.max(), 60.97483) - npt.assert_allclose(data.total_field_anomaly_nt.min(), -3735) - npt.assert_allclose(data.total_field_anomaly_nt.max(), 2792) - npt.assert_allclose(data.altitude_m.min(), 201.0) - npt.assert_allclose(data.altitude_m.max(), 1545.0) - assert set(data.survey_area.unique()) == { - "CA55_NORTH", - "CA55_SOUTH", - "CA57", - "CA58", - "CA59", - "CA60", - "CA63", - "HG56", - "HG57", - "HG58", - "HG61", - "HG62", - "HG64", - "HG65", - } - - -@pytest.mark.filterwarnings( - "ignore:The 'datasets' module will be deprecated:FutureWarning" -) -def test_south_africa_gravity(): - "Sanity checks for the loaded dataset" - data = fetch_south_africa_gravity() - assert data.shape == (14559, 4) - npt.assert_allclose(data.longitude.min(), 11.90833) - npt.assert_allclose(data.longitude.max(), 32.74667) - npt.assert_allclose(data.latitude.min(), -34.996) - npt.assert_allclose(data.latitude.max(), -17.33333) - npt.assert_allclose(data.elevation.min(), -1038.00) - npt.assert_allclose(data.elevation.max(), 2622.17) - npt.assert_allclose(data.gravity.min(), 978131.3) - npt.assert_allclose(data.gravity.max(), 979766.65) - - -@pytest.mark.filterwarnings( - "ignore:The 'datasets' module will be deprecated:FutureWarning" -) -def test_south_africa_topography(): - "Sanity checks for the loaded dataset" - data = fetch_south_africa_topography() - assert data.topography.shape == (171, 211) - npt.assert_allclose(data.topography.min(), -5007.0, atol=0.1) - npt.assert_allclose(data.topography.max(), 3286.0, atol=0.1) - npt.assert_allclose(data.longitude.min(), 12) - npt.assert_allclose(data.longitude.max(), 33) - npt.assert_allclose(data.latitude.min(), -35) - npt.assert_allclose(data.latitude.max(), -18) diff --git a/harmonica/tests/test_synthetic_surveys.py b/harmonica/tests/test_synthetic_surveys.py deleted file mode 100644 index e68414bf9..000000000 --- a/harmonica/tests/test_synthetic_surveys.py +++ /dev/null @@ -1,139 +0,0 @@ -# Copyright (c) 2018 The Harmonica Developers. -# Distributed under the terms of the BSD 3-Clause License. -# SPDX-License-Identifier: BSD-3-Clause -# -# This code is part of the Fatiando a Terra project (https://www.fatiando.org) -# -""" -Test functions to create synthetic surveys -""" -import numpy.testing as npt -import pytest - -from ..synthetic import airborne_survey, ground_survey - - -@pytest.mark.parametrize("synthetic_survey", (airborne_survey, ground_survey)) -def test_deprecation_warning(synthetic_survey): - """ - Checks if deprecation warning is raised - """ - message = "The 'synthetic' module will be deprecated in Harmonica v0.6.0." - with pytest.warns(FutureWarning, match=message): - synthetic_survey() - - -@pytest.mark.filterwarnings( - "ignore:The '(synthetic|datasets)' module will be deprecated:FutureWarning" -) -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) - - -@pytest.mark.filterwarnings( - "ignore:The '(synthetic|datasets)' module will be deprecated:FutureWarning" -) -def test_scale_ground_survey(): - """ - Test if 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) - - -@pytest.mark.filterwarnings( - "ignore:The '(synthetic|datasets)' module will be deprecated:FutureWarning" -) -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) - - -@pytest.mark.filterwarnings( - "ignore:The '(synthetic|datasets)' module will be deprecated:FutureWarning" -) -def test_scale_airborne_survey(): - """ - Test if 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) - - -@pytest.mark.filterwarnings( - "ignore:The '(synthetic|datasets)' module will be deprecated:FutureWarning" -) -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 - - -@pytest.mark.filterwarnings( - "ignore:The '(synthetic|datasets)' module will be deprecated:FutureWarning" -) -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 diff --git a/setup.cfg b/setup.cfg index 11ba9c675..a9d43db09 100644 --- a/setup.cfg +++ b/setup.cfg @@ -45,7 +45,6 @@ install_requires = scipy>=1.5 scikit-learn>=0.24 numba>=0.52 - pooch>=1.2 xarray>=0.16 verde>=1.7 xrft>=1.0