"""
Functions to download the Earth relief datasets from the GMT data server.
The grids are available in various resolutions.
"""
import xarray as xr
from .. import which
from ..exceptions import GMTInvalidInput
[docs]def load_earth_relief(resolution="01d", registration=None):
"""
Load Earth relief grids (topography and bathymetry) in various resolutions.
The grids are downloaded to a user data directory (usually ``~/.gmt/``) the
first time you invoke this function. Afterwards, it will load the data from
the cache. So you'll need an internet connection the first time around.
These grids can also be accessed by passing in the file name
``'@earth_relief_XXm'`` or ``'@earth_relief_XXs'`` to any grid
plotting/processing function.
Parameters
----------
resolution : str
The grid resolution. The suffix ``d``, ``m`` and ``s`` stand for
arc-degree, arc-minute and arc-second. It can be ``'01d'``, ``'30m'``,
``'20m'``, ``'15m'``, ``'10m'``, ``'06m'``, ``'05m'``, ``'04m'``,
``'03m'``, ``'02m'``, ``'01m'``, ``'30s'`` or ``'15s'``.
registration : str
Grid registration type. Either ``pixel`` for pixel registration or
``gridline`` for gridline registration. Default is ``None``, where
a pixel-registered grid is returned unless only the
gridline-registered grid is available.
Returns
-------
grid : xarray.DataArray
The Earth relief grid. Coordinates are latitude and longitude in
degrees. Relief is in meters.
"""
_is_valid_resolution(resolution)
if registration in ("pixel", "gridline", None):
# If None, let GMT decide on Pixel/Gridline type
reg = f"_{registration[0]}" if registration else ""
else:
raise GMTInvalidInput(
f"Invalid grid registration: {registration}, should be either "
"'pixel', 'gridline' or None. Default is None, where a "
"pixel-registered grid is returned unless only the "
"gridline-registered grid is available."
)
fname = which(f"@earth_relief_{resolution}{reg}", download="a")
grid = xr.open_dataarray(fname)
# Add some metadata to the grid
grid.name = "elevation"
grid.attrs["long_name"] = "elevation relative to the geoid"
grid.attrs["units"] = "meters"
grid.attrs["vertical_datum"] = "EMG96"
grid.attrs["horizontal_datum"] = "WGS84"
# Remove the actual range because it gets outdated when indexing the grid,
# which causes problems when exporting it to netCDF for usage on the
# command-line.
grid.attrs.pop("actual_range")
for coord in grid.coords:
grid[coord].attrs.pop("actual_range")
return grid
def _is_valid_resolution(resolution):
"""
Check if a resolution is valid for the global Earth relief grid.
Parameters
----------
resolution : str
Same as the input for load_earth_relief
Raises
------
GMTInvalidInput
If given resolution is not valid.
Examples
--------
>>> _is_valid_resolution("01d")
>>> _is_valid_resolution("60m")
>>> _is_valid_resolution("5m")
Traceback (most recent call last):
...
pygmt.exceptions.GMTInvalidInput: Invalid Earth relief resolution '5m'.
>>> _is_valid_resolution("15s")
>>> _is_valid_resolution("01s")
Traceback (most recent call last):
...
pygmt.exceptions.GMTInvalidInput: Invalid Earth relief resolution '01s'.
"""
valid_resolutions = ["01d"]
valid_resolutions.extend(
[f"{res:02d}m" for res in [60, 30, 20, 15, 10, 6, 5, 4, 3, 2, 1]]
)
valid_resolutions.extend([f"{res:02d}s" for res in [30, 15]])
if resolution not in valid_resolutions:
raise GMTInvalidInput(
"Invalid Earth relief resolution '{}'.".format(resolution)
)
def _shape_from_resolution(resolution):
"""
Calculate the shape of the global Earth relief grid given a resolution.
Parameters
----------
resolution : str
Same as the input for load_earth_relief
Returns
-------
shape : (nlat, nlon)
The calculated shape.
Examples
--------
>>> _shape_from_resolution('60m')
(181, 361)
>>> _shape_from_resolution('30m')
(361, 721)
>>> _shape_from_resolution('10m')
(1081, 2161)
>>> _shape_from_resolution('30s')
(21601, 43201)
>>> _shape_from_resolution('15s')
(43201, 86401)
"""
_is_valid_resolution(resolution)
unit = resolution[2]
if unit == "d":
seconds = int(resolution[:2]) * 60 * 60
elif unit == "m":
seconds = int(resolution[:2]) * 60
elif unit == "s":
seconds = int(resolution[:2])
nlat = 180 * 60 * 60 // seconds + 1
nlon = 360 * 60 * 60 // seconds + 1
return (nlat, nlon)