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

seisman · 2020-07-22T06:45:09Z

Description of proposed changes

See #489 for the reason of the function redesign. In summary, Earth relief data >05m are stored as single grids, and we can use the which function to get its full path and load it as DataArray. It doesn't work for earth relief data <=05m, which are split into smaller tiles. As suggested in #489, it's easier to use grdcut @earth_relief_xxy -Rregion -Goutput.grd. The grdcut function works for all resolutions, but due to the issue #524, the DataArray returned doesn't support slice operations.

Changes in this PR:

For grid>05m, use the which function. The returned DataArray still supports slice operations, so all tests pass
For grids<=05m, use the grdcut function. The returned DataArray doesn't support slice operation (see GMTDataArrayAccessor doesn't work for temporary files that have been sliced #524). The bug is also documented.
For grids<=05m, users have to pass the region argument to select a subregion.
Remove the unused function _shape_from_resolution
Use two lists tiled_resolutions and non_tiled_resolutions, instead of the old _is_valid_resolution function
Remove 01s and 03s from the test, because these two resolutions are valid
Add a few examples
Add tests for high resolution grids

Fixes #489.

Reminders

Run make format and make check to make sure the code follows the style guide.
Add tests for new features or tests that would have caught the bug that you're fixing.
Add new public functions/methods/classes to doc/api/index.rst.
Write detailed docstrings for all functions/methods.
If adding new functionality, add an example to docstrings or tutorials.

seisman · 2020-07-23T01:41:01Z

@weiji14 I think the PR is ready for review, but the grdcut call is not tested, because it means we have to cache at least a subset of high-resolution grids.

I need some help with:

I added example to this functions, but they're not corrected rendered (https://pygmt-8ovwmg7oz.vercel.app/api/generated/pygmt.datasets.load_earth_relief.html#pygmt.datasets.load_earth_relief)
There is a styling error which I don't understand.

weiji14

@weiji14 I think the PR is ready for review, but the grdcut call is not tested, because it means we have to cache at least a subset of high-resolution grids.

We could add the 05m grid to the cache perhaps? They're 'only' 13MB, and with #530 in place it shouldn't be too hard.

pygmt/datasets/earth_relief.py

Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com>

pygmt/datasets/earth_relief.py

Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com>

weiji14 · 2020-07-23T04:00:12Z

pygmt/datasets/earth_relief.py

+    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)


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.

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.

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.

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.

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.

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?

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.

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.

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.

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.

…tion grid

weiji14 added the enhancement Improving an existing feature label Jul 22, 2020

weiji14 added this to the 0.2.x milestone Jul 22, 2020

seisman force-pushed the load_earth_relief_region branch from 924862e to eab3c58 Compare July 23, 2020 00:47

vercel bot temporarily deployed to Preview July 23, 2020 00:47 Inactive

Redesign the load_earth_relief function

d5f8ce0

seisman force-pushed the load_earth_relief_region branch from 1a91063 to d5f8ce0 Compare July 23, 2020 00:50

vercel bot temporarily deployed to Preview July 23, 2020 00:50 Inactive

The "region" argument is only valid for resolutions>=05m

fceab54

vercel bot temporarily deployed to Preview July 23, 2020 01:12 Inactive

seisman changed the title ~~WIP: Redesign the load_earth_relief function~~ Redesign the load_earth_relief function Jul 23, 2020

vercel bot temporarily deployed to Preview July 23, 2020 01:18 Inactive

seisman added 2 commits July 22, 2020 21:25

Fix the error message

757cadf

Add examples and a test

26125cb

vercel bot temporarily deployed to Preview July 23, 2020 01:40 Inactive

seisman marked this pull request as ready for review July 23, 2020 01:40

weiji14 reviewed Jul 23, 2020

View reviewed changes

pygmt/datasets/earth_relief.py Outdated Show resolved Hide resolved

pygmt/datasets/earth_relief.py Outdated Show resolved Hide resolved

weiji14 reviewed Jul 23, 2020

View reviewed changes

pygmt/datasets/earth_relief.py Outdated Show resolved Hide resolved

Merge branch 'master' into load_earth_relief_region

c5651e1

vercel bot temporarily deployed to Preview July 23, 2020 02:13 Inactive

Update pygmt/datasets/earth_relief.py

7b2abd3

Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com>

vercel bot temporarily deployed to Preview July 23, 2020 02:16 Inactive

seisman added 2 commits July 22, 2020 22:16

Fix the description of "region"

80b8350

Fix the examples in docstrings

8e066f7

vercel bot temporarily deployed to Preview July 23, 2020 02:16 Inactive

Add a test

4b3055c

vercel bot temporarily deployed to Preview July 23, 2020 02:34 Inactive

seisman requested a review from weiji14 July 23, 2020 03:09

weiji14 reviewed Jul 23, 2020

View reviewed changes

pygmt/datasets/earth_relief.py Outdated Show resolved Hide resolved

Improve code style

70498b6

Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com>

vercel bot temporarily deployed to Preview July 23, 2020 03:42 Inactive

weiji14 reviewed Jul 23, 2020

View reviewed changes

Raise an NotImplementedError if subregion is requested for low-resolu…

f50db14

…tion grid

vercel bot temporarily deployed to Preview July 23, 2020 05:31 Inactive

weiji14 approved these changes Jul 23, 2020

View reviewed changes

seisman changed the title ~~Redesign the load_earth_relief function~~ Let load_earth_relief() support all resolutions and optional subregion Jul 23, 2020

seisman merged commit 32c4d4a into master Jul 23, 2020

seisman deleted the load_earth_relief_region branch July 23, 2020 06:00

This was referenced Feb 12, 2021

Low-resolution Earth relief grids should support the "region" argument #870

Closed

Let load_earth_relief() support 'region' argument for all resolutions #873

Merged

weiji14 mentioned this pull request Feb 26, 2023

Remove the notes that high-resolution remote datasets don't support slice operations #2378

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

seisman commented Jul 22, 2020 •

edited

Loading

seisman commented Jul 23, 2020 •

edited

Loading

weiji14 left a comment •

edited

Loading

weiji14 Jul 23, 2020 •

edited

Loading

seisman Jul 23, 2020

seisman Jul 23, 2020

weiji14 Jul 23, 2020

weiji14 Jul 23, 2020 •

edited

Loading

seisman Jul 23, 2020

weiji14 Jul 23, 2020

seisman Jul 23, 2020

seisman Jul 23, 2020 •

edited

Loading

weiji14 Jul 23, 2020

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

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

Conversation

seisman commented Jul 22, 2020 • edited Loading

seisman commented Jul 23, 2020 • edited Loading

weiji14 left a comment • edited Loading

Choose a reason for hiding this comment

weiji14 Jul 23, 2020 • edited Loading

Choose a reason for hiding this comment

seisman Jul 23, 2020

Choose a reason for hiding this comment

seisman Jul 23, 2020

Choose a reason for hiding this comment

weiji14 Jul 23, 2020

Choose a reason for hiding this comment

weiji14 Jul 23, 2020 • edited Loading

Choose a reason for hiding this comment

seisman Jul 23, 2020

Choose a reason for hiding this comment

weiji14 Jul 23, 2020

Choose a reason for hiding this comment

seisman Jul 23, 2020

Choose a reason for hiding this comment

seisman Jul 23, 2020 • edited Loading

Choose a reason for hiding this comment

weiji14 Jul 23, 2020

Choose a reason for hiding this comment

seisman commented Jul 22, 2020 •

edited

Loading

seisman commented Jul 23, 2020 •

edited

Loading

weiji14 left a comment •

edited

Loading

weiji14 Jul 23, 2020 •

edited

Loading

weiji14 Jul 23, 2020 •

edited

Loading

seisman Jul 23, 2020 •

edited

Loading