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.

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

Add gallery example for grdclip #1396

Merged
merged 13 commits into from
Aug 4, 2021
35 changes: 35 additions & 0 deletions examples/gallery/images/grdclip.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
"""
Clip grid values
----------------
michaelgrund marked this conversation as resolved.
Show resolved Hide resolved
The :meth:`pygmt.grdclip` method allows to clip defined ranges of grid values.
In the example shown below we set all elevation values (grid points) smaller
than 0 m (in general the bathymetric part of the grid) to a common value of
-2000 m via the ``below`` parameter.
"""

import pygmt

fig = pygmt.Figure()

# Load sample grid (1 arc minute global relief) and use area around the Hawaiian Islands
grid = pygmt.datasets.load_earth_relief(resolution="01m", region=[-162, -153, 18, 23])
Copy link
Member

Choose a reason for hiding this comment

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

The 1 arc minute global relief data is NOT cached in the CI. We should either cache the 01m data or use a lower resolution earth relief.

Copy link
Member Author

@michaelgrund michaelgrund Jul 29, 2021

Choose a reason for hiding this comment

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

I would prefer to use a high resolution as the 1 arc minute global relief, otherwise the map looks bad. Will add the 01m data to cache in a separate PR (#1398).

Copy link
Member

Choose a reason for hiding this comment

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

Or can we use a different location? I like the Hawaii example, but caching the 1 arc minute data (258MB) isn't very feasible... How about another island like Iceland (or another medium-ish sized one).

Copy link
Member

Choose a reason for hiding this comment

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

Or can we use a different location? I like the Hawaii example, but caching the 1 arc minute data (258MB) isn't very feasible... How about another island like Iceland (or another medium-ish sized one).

Sounds a better solution.

Copy link
Member Author

@michaelgrund michaelgrund Jul 31, 2021

Choose a reason for hiding this comment

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

I changed it to the 3 arc minute grid (27 MB) for the area around Iceland. Is that fine or still to large for caching?

Copy link
Member

Choose a reason for hiding this comment

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

The current cache size is 30.5MB (https://github.com/GenericMappingTools/pygmt/actions/runs/1064597430), so adding 27MB would make it 50.5MB. I'm still a bit hesitant to add it... How about adding just the tile that is downloaded. Running gmt grdcut @earth_relief_03m -R-28/-10/62/68 -Gtemp.nc gives me the tile name "N00W090.earth_relief_03m_p" in .gmt/server/earth/earth_relief/earth_relief_03m_p so maybe try just caching that?

Copy link
Member Author

@michaelgrund michaelgrund Aug 2, 2021

Choose a reason for hiding this comment

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

The current cache size is 30.5MB (https://github.com/GenericMappingTools/pygmt/actions/runs/1064597430), so adding 27MB would make it 50.5MB. I'm still a bit hesitant to add it... How about adding just the tile that is downloaded. Running gmt grdcut @earth_relief_03m -R-28/-10/62/68 -Gtemp.nc gives me the tile name "N00W090.earth_relief_03m_p" in .gmt/server/earth/earth_relief/earth_relief_03m_p so maybe try just caching that?

Does adding "N00W090.earth_relief_03m_p" allow to automatically use only this tile for plotting or is there anything else to consider @weiji14 ?

Copy link
Member

@core-man core-man Aug 2, 2021

Choose a reason for hiding this comment

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

See another similar case:

"@N35E135.earth_relief_03s_g.nc",

So, I think we could only cache such a single tile. See some notes for this case in #985 (comment).

Copy link
Member Author

Choose a reason for hiding this comment

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

See another similar case:

"@N35E135.earth_relief_03s_g.nc",

I think we could only cache such a single tile.

I already recognized this one 😉 . My question is: Is that tile automatically used by pygmt.datasets.load_earth_relief() when the corresponding region is defined?

Copy link
Member

Choose a reason for hiding this comment

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

Is that tile automatically used by pygmt.datasets.load_earth_relief() when the corresponding region is defined?

Yes.


# Plot original grid
fig.basemap(
region=[-162, -153, 18, 23], projection="M12c", frame=["f", '+t"original grid"']
)
fig.grdimage(grid=grid, cmap="oleron")

fig.shift_origin(yshift="-9c")
michaelgrund marked this conversation as resolved.
Show resolved Hide resolved

# Set all grid points < 0 m to a value of -2000 m.
grid = pygmt.grdclip(grid, below=[0, -2000])

# Plot clipped grid
fig.basemap(
region=[-162, -153, 18, 23], projection="M12c", frame=["f", '+t"clipped grid"']
)
fig.grdimage(grid=grid)
fig.colorbar(frame=["x+lElevation", "y+lm"])

fig.show()