Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Let load_earth_relief() support all resolutions and optional subregion #542

Merged
merged 11 commits into from
Jul 23, 2020
Merged

Let load_earth_relief() support all resolutions and optional subregion #542

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
d5f8ce0
Redesign the load_earth_relief function
seisman Jul 10, 2020
fceab54
The "region" argument is only valid for resolutions>=05m
seisman Jul 23, 2020
757cadf
Fix the error message
seisman Jul 23, 2020
26125cb
Add examples and a test
seisman Jul 23, 2020
c5651e1
Merge branch 'master' into load_earth_relief_region
seisman Jul 23, 2020
7b2abd3
Update pygmt/datasets/earth_relief.py
seisman Jul 23, 2020
80b8350
Fix the description of "region"
seisman Jul 23, 2020
8e066f7
Fix the examples in docstrings
seisman Jul 23, 2020
4b3055c
Add a test
seisman Jul 23, 2020
70498b6
Improve code style
seisman Jul 23, 2020
f50db14
Raise an NotImplementedError if subregion is requested for low-resolu…
seisman Jul 23, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
169 changes: 73 additions & 96 deletions pygmt/datasets/earth_relief.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,32 +1,39 @@
"""
Functions to download the Earth relief datasets from the GMT data server.
The grids are available in various resolutions.
Function to download the Earth relief datasets from the GMT data server,
and load as DataArray. The grids are available in various resolutions.
"""
import xarray as xr

from .. import which
from .. import grdcut, which
from ..exceptions import GMTInvalidInput
from ..helpers import kwargs_to_strings


def load_earth_relief(resolution="01d", registration=None):
@kwargs_to_strings(region="sequence")
def load_earth_relief(resolution="01d", region=None, 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.
The grids are downloaded to a user data directory
(usually ``~/.gmt/server/earth/earth_relief/``) the first time you invoke
this function. Afterwards, it will load the grid from the data directory.
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.
``'@earth_relief_rru[_reg]'`` to any grid plotting/processing function.
Refer to :gmt-docs:`datasets/remote-data.html` for more details.

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'``.
``'03m'``, ``'02m'``, ``'01m'``, ``'30s'``, ``'15s'``, ``'03s'``,
or ``'01s'``.

region : str or list
The subregion of the grid to load. Only works
seisman marked this conversation as resolved.
Show resolved Hide resolved

registration : str
Grid registration type. Either ``pixel`` for pixel registration or
Expand All @@ -40,8 +47,46 @@ def load_earth_relief(resolution="01d", registration=None):
The Earth relief grid. Coordinates are latitude and longitude in
degrees. Relief is in meters.

Notes
-----
The DataArray doesn's support slice operation, for Earth relief data with
resolutions higher than "05m", which are stored as smaller tiles.

Examples
--------

# load the default grid (1 arc degree resolution and pixel registration)
>>> grid = load_earth_relief()

# load the 30m grid with "gridline" registration
>>> grid = load_earth_relief("30m", registration="gridline")

# load high-resolution grid for a specific region
>>> grid = load_earth_relief("01m", region=[120, 140, 30, 50])
seisman marked this conversation as resolved.
Show resolved Hide resolved

"""
_is_valid_resolution(resolution)

# earth relief data stored as single grids for low resolutions
non_tiled_resolutions = [
"01d",
"30m",
"20m",
"15m",
"10m",
"06m",
]
# earth relief data stored as tiles for high resolutions
tiled_resolutions = [
"05m",
"04m",
"03m",
"02m",
"01m",
"30s",
"15s",
"03s",
"01s",
]
seisman marked this conversation as resolved.
Show resolved Hide resolved

if registration in ("pixel", "gridline", None):
# If None, let GMT decide on Pixel/Gridline type
Expand All @@ -54,8 +99,23 @@ def load_earth_relief(resolution="01d", registration=None):
"gridline-registered grid is available."
)

fname = which(f"@earth_relief_{resolution}{reg}", download="a")
grid = xr.open_dataarray(fname)
# different ways to load tiled and non-tiled earth relief data
if resolution in non_tiled_resolutions:
fname = which(f"@earth_relief_{resolution}{reg}", download="a")
with xr.open_dataarray(fname) as dataarray:
grid = dataarray.load()
_ = grid.gmt # load GMTDataArray accessor information
elif resolution in tiled_resolutions:
# Titled grid can't be sliced.
# See https://github.com/GenericMappingTools/pygmt/issues/524
if region is None:
raise GMTInvalidInput(
"region is required for high resolution (<='05m') Earth relief grid"
)
grid = grdcut(f"@earth_relief_{resolution}{reg}", region=region)
Copy link
Member

@weiji14 weiji14 Jul 23, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, one more thing. What if a user wants to get a non_tiled_resolution grid with a region? The current code doesn't seem to handle this case.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The initial design was to allow users to specify a subregion for all resolutions. However, due to #524, the returned grids can't support slice operation, which makes three or more tests fail. So I have to do different ways for tiled and non-tiled resolutions, to keep all tests pass.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If a user really want it, they need to call load_earth_relief, then grdcut, or using xarray slice operation.

I'm also OK to add "region" to all resolutions, and declare that "load_earth_relief() grids don't support slice operation" as a bug.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right, that's tricky. But, an unsuspecting user might pass in a region argument for resolutions 06m and above, and it would just be silently ignored (and a full grid is returned)! We should at least raise a warning or NotImplementedError for this case.

Copy link
Member

@weiji14 weiji14 Jul 23, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, slow reply. A workaround might be to just set the registration and gtype directly after calling grdcut, since we know it from the user input to load_earth_relief.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

that 05m(and higher resolution) grids which are grdcut won't work with grdimage or grdview, but we can fix this before doing a proper release.

I don't understand it. Why do you say it doesn't work?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand it. Why do you say it doesn't work?

Sorry, it will work actually, I forgot for a second that we applied the workaround in #539.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, if we remove/pop the temp.nc file from grid.encoding["source"] after running grdcut, we can then set the registration and gtype without getting an error. Do you think it's wise to apply this hack?

You mean applying the hack only in load_earth_relief() function? Sounds a good idea.

Copy link
Member Author

@seisman seisman Jul 23, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Implemented locally, but I'm getting this error:

_________________________________________ [doctest] pygmt.clib.conversion.dataarray_to_matrix _________________________________________
060     >>> matrix.flags.c_contiguous
061     True
062     >>> # Using a slice of the grid, the matrix will be copied to guarantee
063     >>> # that it's C-contiguous in memory. The increment should be unchanged.
064     >>> matrix, region, inc = dataarray_to_matrix(grid[10:41,30:101])
065     >>> matrix.flags.c_contiguous
066     True
067     >>> print(matrix.shape)
068     (31, 71)
069     >>> print(region)
Expected:
    [-150.0, -79.0, -80.0, -49.0]
Got:
    [-149.5, -79.5, -79.5, -49.5]

This is because grid.encoding["source"] is removed, then the sliced new grid has the default registration and gtype.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't see an easy fix for #524. Perhaps the best compromise is

raise a warning or NotImplementedError for this case.

RIght, let's just raise a warning here then, having layers and layers of workarounds isn't good anyway.

else:
raise GMTInvalidInput(f"Invalid Earth relief resolution '{resolution}'")
seisman marked this conversation as resolved.
Show resolved Hide resolved

# Add some metadata to the grid
grid.name = "elevation"
grid.attrs["long_name"] = "elevation relative to the geoid"
Expand All @@ -69,86 +129,3 @@ def load_earth_relief(resolution="01d", registration=None):
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)
8 changes: 7 additions & 1 deletion pygmt/tests/test_datasets.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ def test_usgs_quakes():

def test_earth_relief_fails():
"Make sure earth relief fails for invalid resolutions"
resolutions = "1m 1d bla 60d 01s 03s 001m 03".split()
resolutions = "1m 1d bla 60d 001m 03".split()
resolutions.append(60)
for resolution in resolutions:
with pytest.raises(GMTInvalidInput):
Expand Down Expand Up @@ -92,3 +92,9 @@ def test_earth_relief_incorrect_registration():
"Test loading earth relief with incorrect registration type"
with pytest.raises(GMTInvalidInput):
load_earth_relief(registration="improper_type")


def test_earth_relief_none_region():
"Test loading high-resolution earth relief without passing 'region'"
with pytest.raises(GMTInvalidInput):
load_earth_relief("05m")