-
Notifications
You must be signed in to change notification settings - Fork 69
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Add global topography data and update gravity data (#23)
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.
- Loading branch information
Showing
22 changed files
with
290 additions
and
15 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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: |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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() |
Binary file not shown.
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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>`__ |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,2 @@ | ||
# pylint: disable=missing-docstring | ||
from .sample_data import fetch_gravity_earth, fetch_topography_earth |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,2 @@ | ||
etopo1-0.5deg.nc.xz bb44b4aaa49778159dc3655ae0c0c23ad88789328829557ba7ea0113f26baa7f | ||
gravity-earth-0.5deg.nc.xz 84340b771ac91695b0874639a1497b6b81e0f5a8d21ac51cc221d1c48598f5a7 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,87 @@ | ||
""" | ||
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) | ||
# Call load to make sure the data are loaded into memory and not linked to file | ||
grid = xr.open_dataset(decompressed.name, **kwargs).load() | ||
# Close any files associated with this dataset to make sure can delete them | ||
grid.close() | ||
finally: | ||
os.remove(decompressed.name) | ||
return grid |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,24 @@ | ||
""" | ||
Test the sample data loading functions. | ||
""" | ||
import numpy.testing as npt | ||
|
||
from ..datasets.sample_data import fetch_gravity_earth, fetch_topography_earth | ||
|
||
|
||
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) | ||
|
||
|
||
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(), 5622) | ||
npt.assert_allclose(grid.topography.min(), -8397) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters