Add global topography data and update gravity data

Use a 0.5 degree ETOPO1 dataset as Earth topography (will be used in
isostasy examples and Bouguer gravity). Update the global gravity data
to use EIGEN-6C4 at constant height to avoid dealing with conversion
between orthometric and geometric heights. The topography grid is in
geometric heights converted using an EIGEN-6C4 geoid.

Add functions to load these datasets using Pooch and examples for each.
Add an example calculating the gravity disturbance using the global
gravity data.

MANIFEST.in

@@ -5,5 +5,6 @@ include CONTRIBUTING.md
 include MAINTENANCE.md
 include AUTHORS.md
 include versioneer.py
+include harmonica/datasets/registry.txt
 recursive-include harmonica/tests/data *
 recursive-include harmonica/tests/baseline *

Makefile

@@ -3,8 +3,8 @@ PROJECT=harmonica
 TESTDIR=tmp-test-dir-with-unique-name
 PYTEST_ARGS=--cov-config=../.coveragerc --cov-report=term-missing --cov=$(PROJECT) --doctest-modules -v --pyargs
 LINT_FILES=setup.py $(PROJECT)
-BLACK_FILES=setup.py $(PROJECT) examples
-FLAKE8_FILES=setup.py $(PROJECT) examples
+BLACK_FILES=setup.py $(PROJECT) examples data/examples
+FLAKE8_FILES=setup.py $(PROJECT) examples data/examples
 
 help:
  @echo "Commands:"

data/README.md

@@ -2,12 +2,16 @@
 
 These files are used as sample data in Harmonica:
 
-* `global_gravity_earth.nc.xz`: Global gravity data set. The file contains the
- magnitude of the gradient of the gravitational potential of the Earth
- (including the centrifugal potential) on a regular grid of points located on
- the Earth surface. The file was downloaded from the [ICGEM Calculation
- Service](http://icgem.gfz-potsdam.de/calcgrid) using the EIGEN-6C4 gravity
- field model and the WGS84 Reference System. The data is stored in a netCDF
- file and `xz` compressed. The netCDF file can be loaded through the
- `xarray.open_dataset()` function, obtaining a `xarray.Dataset` object.
-
+* `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 using Python's `lzma` library.
+* `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 original heights referenced to the
+ sea level were converted to geometric (ellipsoidal) heights using a geoid grid derived
+ from the EIGEN-6C4 gravity field model. Both grids were generated by the
+ [ICGEM Calculation Service](http://icgem.gfz-potsdam.de). The data are stored in a
+ netCDF file and then `xz` compressed using Python's `lzma` library.

data/examples/README.txt

@@ -0,0 +1,29 @@
+.. _sample_data:
+
+Sample Data
+===========
+
+Harmonica provides some sample data for testing through the :mod:`harmonica.datasets`
+module. The sample data are automatically downloaded from the `Github repository
+<https:/fatiando/harmonica>`__ to a folder on your computer the first time
+you use them. After that, the data are loaded from this folder. The download is managed
+by the :mod:`pooch` package.
+
+
+Where is my data?
+-----------------
+
+The data files are downloaded to a folder ``~/.harmonica/data/`` by default. This is the
+*base data directory*. :mod:`pooch` will create a separate folder in the base directory
+for each version of Harmonica. For example, the base data dir for v0.1.0 is
+``~/.harmonica/data/v0.1.0``. If you're using the latest development version from
+Github, the version is ``master``.
+
+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:

data/examples/gravity_earth.py

@@ -0,0 +1,29 @@
+"""
+Gravity of the Earth
+====================
+
+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 which was
+generated from the spherical harmonic model EIGEN-6C4 [Forste_etal2014]_.
+"""
+import matplotlib.pyplot as plt
+import cartopy.crs as ccrs
+import harmonica as hm
+
+# Load the gravity grid
+data = hm.datasets.fetch_gravity_earth()
+print(data)
+
+# Make a plot of data using Cartopy
+plt.figure(figsize=(6, 7))
+ax = plt.axes(projection=ccrs.Orthographic(central_longitude=150))
+pc = data.gravity.plot.pcolormesh(
+ ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False
+)
+plt.colorbar(
+ pc, label="mGal", orientation="horizontal", aspect=50, pad=0.005, shrink=0.7
+)
+ax.set_title("Gravity of the Earth (EIGEN-6C4)")
+ax.coastlines()
+plt.tight_layout()
+plt.show()

data/examples/topography_earth.py

@@ -0,0 +1,29 @@
+"""
+Topography of the Earth
+=======================
+
+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.
+"""
+import matplotlib.pyplot as plt
+import cartopy.crs as ccrs
+import harmonica as hm
+
+# Load the topography grid
+data = hm.datasets.fetch_topography_earth()
+print(data)
+
+# Make a plot of data using Cartopy
+plt.figure(figsize=(6, 7))
+ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-30))
+pc = data.topography.plot.pcolormesh(
+ ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False, cmap="terrain"
+)
+plt.colorbar(
+ pc, label="meters", orientation="horizontal", aspect=50, pad=0.005, shrink=0.7
+)
+ax.set_title("Topography of the Earth (ETOPO1)")
+ax.coastlines()
+plt.tight_layout()
+plt.show()

doc/api/index.rst

@@ -34,6 +34,20 @@ Normal Gravity
 
  normal_gravity
 
+
+.. automodule:: harmonica.datasets
+
+.. currentmodule:: harmonica
+
+Datasets
+--------
+
+.. autosummary::
+ :toctree: generated/
+
+ datasets.fetch_gravity_earth
+ datasets.fetch_topography_earth
+
 Utilities
 ---------
 

doc/conf.py

@@ -46,9 +46,9 @@
 
 sphinx_gallery_conf = {
  # path to your examples scripts
- 'examples_dirs': ['../examples'],
+ 'examples_dirs': ['../examples', '../data/examples'],
  # path where to save gallery generated examples
- 'gallery_dirs': ['gallery'],
+ 'gallery_dirs': ['gallery', 'sample_data'],
  'filename_pattern': '\.py',
  # Remove the "Download all examples" button from the top level gallery
  'download_all_examples': False,

doc/index.rst

@@ -4,7 +4,7 @@
 
  <div class="banner">
  <img src="_static/readme-banner.png" alt="harmonica">
- <h2> Forward modeling, inversion, and processing gravity and magnetic data </h2>
+ <h2>Forward modeling, inversion, and processing gravity and magnetic data</h2>
  <p>A part of the <a href="https://www.fatiando.org/">Fatiando a Terra</a> project.</p>
  </div>
 
@@ -24,5 +24,12 @@
  :hidden:
  :caption: User Guide
 
+ sample_data/index.rst
+
+.. toctree::
+ :maxdepth: 2
+ :hidden:
+ :caption: Reference documentation
+
  api/index.rst
  references.rst

doc/references.rst

@@ -1,6 +1,8 @@
 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 <https://doi.org/10.7289/V5C8276M>`__
 .. [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 <https://doi.org/10.1007/s00190-016-0948-z>`__
+.. [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 <http://doi.org/10.5880/icgem.2015.1>`__
 .. [Hofmann-WellenhofMoritz2006] Hofmann-Wellenhof, B., & Moritz, H. (2006). Physical Geodesy (2nd, corr. ed. 2006 edition ed.). Wien ; New York: Springer.
 .. [LiGotze2001] Li, X. and H. J. Gotze, 2001, Tutorial: Ellipsoid, geoid, gravity, geodesy, and geophysics, Geophysics, 66(6), p. 1660-1668, doi:`10.1190/1.1487109 <https://doi.org/10.1190/1.1487109>`__

examples/gravity_disturbance.py

@@ -0,0 +1,35 @@
+"""
+Gravity Disturbances
+====================
+
+Gravity disturbances are the differences between the measured gravity and a reference
+(normal) gravity produced by an ellipsoid. The disturbances are what allows
+geoscientists to infer the internal structure of the Earth. We'll use the
+:func:`harmonica.normal_gravity` function to calculate the global gravity disturbance of
+the Earth using our sample gravity data.
+"""
+import matplotlib.pyplot as plt
+import cartopy.crs as ccrs
+import harmonica as hm
+
+# Load the global gravity grid
+data = hm.datasets.fetch_gravity_earth()
+print(data)
+
+# Calculate normal gravity and the disturbance
+gamma = hm.normal_gravity(data.latitude, data.height_over_ell)
+data["disturbance"] = data.gravity - gamma
+
+# Make a plot of data using Cartopy
+plt.figure(figsize=(7, 8))
+ax = plt.axes(projection=ccrs.Orthographic(central_longitude=100))
+pc = data.disturbance.plot.pcolormesh(
+ ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False
+)
+plt.colorbar(
+ pc, label="mGal", orientation="horizontal", aspect=50, pad=0.005, shrink=0.7
+)
+ax.set_title("Gravity of disturbance of the Earth")
+ax.coastlines()
+plt.tight_layout()
+plt.show()

harmonica/__init__.py

@@ -9,6 +9,7 @@
 )
 from .io import load_icgem_gdf
 from .normal_gravity import normal_gravity
+from . import datasets
 
 # Get the version number through versioneer
 __version__ = version.full_version

harmonica/datasets/__init__.py

@@ -0,0 +1,2 @@
+# pylint: disable=missing-docstring
+from .sample_data import fetch_gravity_earth, fetch_topography_earth

harmonica/datasets/registry.txt

@@ -0,0 +1,2 @@
+etopo1-0.5deg.nc.xz bb44b4aaa49778159dc3655ae0c0c23ad88789328829557ba7ea0113f26baa7f
+gravity-earth-0.5deg.nc.xz 84340b771ac91695b0874639a1497b6b81e0f5a8d21ac51cc221d1c48598f5a7

harmonica/datasets/sample_data.py

@@ -0,0 +1,84 @@
+"""
+Functions to load sample datasets used in the Harmonica docs.
+"""
+import os
+import tempfile
+import lzma
+import shutil
+
+import xarray as xr
+import pooch
+
+from ..version import full_version
+
+POOCH = pooch.create(
+ path=["~", ".harmonica", "data"],
+ base_url="https:/fatiando/harmonica/raw/{version}/data/",
+ version=full_version,
+ version_dev="master",
+ env="HARMONICA_DATA_DIR",
+)
+POOCH.load_registry(os.path.join(os.path.dirname(__file__), "registry.txt"))
+
+
+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
+ <http://icgem.gfz-potsdam.de/>`__. 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 height grid
+ (``height_over_ell``). Coordinates are latitude and longitude.
+
+ """
+ return _load_xz_compressed_grid(POOCH.fetch("gravity-earth-0.5deg.nc.xz"))
+
+
+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
+ <http://icgem.gfz-potsdam.de/>`__. 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 topography grid (in meters). Coordinates are latitude and longitude.
+
+ """
+ return _load_xz_compressed_grid(POOCH.fetch("etopo1-0.5deg.nc.xz"))
+
+
+def _load_xz_compressed_grid(fname, **kwargs):
+ """
+ Load a netCDF grid that has been xz compressed. Keyword arguments are passed to
+ :func:`xarray.open_dataset`.
+ """
+ decompressed = tempfile.NamedTemporaryFile(suffix=".nc", delete=False)
+ try:
+ with decompressed:
+ with lzma.open(fname, "rb") as compressed:
+ shutil.copyfileobj(compressed, decompressed)
+ grid = xr.open_dataset(decompressed.name, **kwargs)
+ finally:
+ os.remove(decompressed.name)
+ return grid

harmonica/normal_gravity.py

@@ -31,14 +31,20 @@ def normal_gravity(latitude, height):
  """
  ellipsoid = get_ellipsoid()
  ecc = ellipsoid.linear_eccentricity
- omega2 = ellipsoid.angular_velocity**2
+ omega2 = ellipsoid.angular_velocity ** 2
 
  big_w, q_l, q_0, b_l, cosbeta_l2, sinbeta_l2 = _normal_grav_factors(
- ellipsoid, latitude, height)
- term1 = ellipsoid.geocentric_grav_const / (b_l**2 + ecc**2)
+ ellipsoid, latitude, height
+ )
+ term1 = ellipsoid.geocentric_grav_const / (b_l ** 2 + ecc ** 2)
  term2 = 0.5 * sinbeta_l2 - 1 / 6
- term2 *= ellipsoid.semimajor_axis**2 * ecc * q_l * omega2 / \
- ((b_l**2 + ecc**2) * q_0)
+ term2 *= (
+ ellipsoid.semimajor_axis ** 2
+ * ecc
+ * q_l
+ * omega2
+ / ((b_l ** 2 + ecc ** 2) * q_0)
+ )
  term3 = -cosbeta_l2 * b_l * omega2
  gamma = (term1 + term2 + term3) / big_w
  # Convert gamma from SI to mGal
@@ -50,16 +56,17 @@ def _normal_grav_factors(ellipsoid, latitude, height):
  ecc = ellipsoid.linear_eccentricity
  semiminor = ellipsoid.semiminor_axis
  rl2, zl2 = _normal_grav_r_z(ellipsoid, latitude, height)
- big_d = (rl2 - zl2) / ecc**2
- big_r = (rl2 + zl2) / ecc**2
- cosbeta_l2 = 0.5 * (1 + big_r) - \
- np.sqrt(0.25 * (1 + big_r**2) - 0.5 * big_d)
+ big_d = (rl2 - zl2) / ecc ** 2
+ big_r = (rl2 + zl2) / ecc ** 2
+ cosbeta_l2 = 0.5 * (1 + big_r) - np.sqrt(0.25 * (1 + big_r ** 2) - 0.5 * big_d)
  sinbeta_l2 = 1 - cosbeta_l2
- b_l = np.sqrt(rl2 + zl2 - ecc**2 * cosbeta_l2)
- q_0 = 0.5 * ((1 + 3 * (semiminor / ecc)**2) * np.arctan2(ecc, semiminor)
- - 3 * semiminor / ecc)
- q_l = 3 * (1 + (b_l / ecc)**2) * (1 - b_l / ecc * np.arctan2(ecc, b_l)) - 1
- big_w = np.sqrt((b_l**2 + ecc**2 * sinbeta_l2) / (b_l**2 + ecc**2))
+ b_l = np.sqrt(rl2 + zl2 - ecc ** 2 * cosbeta_l2)
+ q_0 = 0.5 * (
+ (1 + 3 * (semiminor / ecc) ** 2) * np.arctan2(ecc, semiminor)
+ - 3 * semiminor / ecc
+ )
+ q_l = 3 * (1 + (b_l / ecc) ** 2) * (1 - b_l / ecc * np.arctan2(ecc, b_l)) - 1
+ big_w = np.sqrt((b_l ** 2 + ecc ** 2 * sinbeta_l2) / (b_l ** 2 + ecc ** 2))
  return big_w, q_l, q_0, b_l, cosbeta_l2, sinbeta_l2
 
 
@@ -70,6 +77,6 @@ def _normal_grav_r_z(ellipsoid, latitude, height):
  coslat, sinlat = np.cos(lat), np.sin(lat)
  beta = np.arctan2(semiminor * sinlat / coslat, semimajor)
  cosbeta, sinbeta = np.cos(beta), np.sin(beta)
- zl2 = (semiminor * sinbeta + height * sinlat)**2
- rl2 = (semimajor * cosbeta + height * coslat)**2
+ zl2 = (semiminor * sinbeta + height * sinlat) ** 2
+ rl2 = (semimajor * cosbeta + height * coslat) ** 2
  return rl2, zl2