From 095c795aabd94966b40d657bfea5b0f0064db24f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yvonne=20Fr=C3=B6hlich?= Date: Thu, 9 Apr 2026 18:33:14 +0200 Subject: [PATCH 01/12] Add first draft for gallery example or tutorial on grdmask --- examples/gallery/images/grdmask.py | 60 ++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100755 examples/gallery/images/grdmask.py diff --git a/examples/gallery/images/grdmask.py b/examples/gallery/images/grdmask.py new file mode 100755 index 00000000000..bdf55753b7f --- /dev/null +++ b/examples/gallery/images/grdmask.py @@ -0,0 +1,60 @@ +""" +Create grid mask from polygons +============================== +:func:`pygmt.grdmask`. + +:func:`pygmt.grdlandmask` and gallery example https://www.pygmt.org/latest/gallery/images/grdlandmask.html. +""" + +# %% +import numpy as np +import pygmt + +# Define a study region +region = [120, 135, 25, 40] + +# Define two polygons, here triangels, use nan to separate the polygons +polygon = np.array( + [ + [125, 30], + [130, 30], + [130, 35], + [125, 30], + [np.nan, np.nan], + [122, 37], + [128, 37], + [128, 39], + [122, 37], + ], +) + +# Download elevation grid +grid = pygmt.datasets.load_earth_relief(region=region, resolution="30s") + +# Create a grid mask based on the two polygons defined above, set all values +# outside the polygons to NaN +mask = pygmt.grdmask(region=region, data=polygon, spacing="30s", outside="NaN") + +# Apply the grid mask to the downloaded elevation grid by multiplying the two grids +grid_mask = grid * mask + + +fig = pygmt.Figure() +pygmt.makecpt(cmap="oleron", series=[-2000, 2000]) + +# Plot the elevation grid +fig.basemap(region=region, projection="M12c", frame=True) +fig.grdimage(grid=grid, cmap=True) +fig.basemap(frame="g1") +fig.plot(data=polygon, pen="2p,darkorange") + +fig.shift_origin(xshift="+w+2c") + +# Plot the masked elevation grid +fig.basemap(region=region, projection="M12c", frame=True) +fig.grdimage(grid=grid_mask, cmap=True) +fig.basemap(frame="g1") +fig.plot(data=polygon, pen="2p,darkorange") + +fig.colorbar(frame=True) +fig.show() From 5456a7fc2fe4d049afa0d596e44add554b538d11 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yvonne=20Fr=C3=B6hlich?= Date: Thu, 9 Apr 2026 20:20:44 +0200 Subject: [PATCH 02/12] Adjust regions --- examples/gallery/images/grdmask.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/examples/gallery/images/grdmask.py b/examples/gallery/images/grdmask.py index bdf55753b7f..32f2abe6fff 100755 --- a/examples/gallery/images/grdmask.py +++ b/examples/gallery/images/grdmask.py @@ -11,20 +11,22 @@ import pygmt # Define a study region -region = [120, 135, 25, 40] +region = [125, 135, 25, 36] -# Define two polygons, here triangels, use nan to separate the polygons +# Define two closed polygons, here a quare and a triangle. +# Use nan to separate the polygons polygon = np.array( [ - [125, 30], - [130, 30], - [130, 35], - [125, 30], + [129, 31], + [134, 31], + [134, 35], + [129, 35], + [129, 31], [np.nan, np.nan], - [122, 37], - [128, 37], - [128, 39], - [122, 37], + [126, 26], + [131, 26], + [131, 30], + [126, 26], ], ) From 31d1b74c505520515ef96eb8d2b7f1bfd873ea50 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yvonne=20Fr=C3=B6hlich?= Date: Thu, 9 Apr 2026 20:52:21 +0200 Subject: [PATCH 03/12] Try extracting a circle --- examples/gallery/images/grdmask.py | 52 ++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/examples/gallery/images/grdmask.py b/examples/gallery/images/grdmask.py index 32f2abe6fff..c65fffb6408 100755 --- a/examples/gallery/images/grdmask.py +++ b/examples/gallery/images/grdmask.py @@ -7,8 +7,14 @@ """ # %% +import geopandas import numpy as np import pygmt +from shapely.geometry import Point + +# %% +# Polygons based on NumPy arrays +# ------------------------------ # Define a study region region = [125, 135, 25, 36] @@ -60,3 +66,49 @@ fig.colorbar(frame=True) fig.show() + + +# %% +# Circle based on GeoPandas polygon geometry +# ------------------------------------------ + +# Define a study region +region = [125, 135, 25, 36] + +# Create a point and buffer it +point = geopandas.GeoSeries([Point(126.5, 33.5)]) +circle = point.buffer(0.6) # 10 is the radius + +# Download elevation grid +grid = pygmt.datasets.load_earth_relief(region=region, resolution="30s") + +# Create a grid mask based on the two polygons defined above, set all values +# outside the polygons to NaN +mask = pygmt.grdmask(region=region, data=circle, spacing="30s", outside="NaN") +mask_lonlat = mask.rename(new_name_or_name_dict={"x": "lon", "y": "lat"}) + +# Apply the grid mask to the downloaded elevation grid by multiplying the two grids +grid_mask = grid * mask_lonlat + + +fig = pygmt.Figure() +pygmt.makecpt(cmap="oleron", series=[-2000, 2000]) + +# Plot the elevation grid +fig.basemap(region=region, projection="M12c", frame=True) +fig.grdimage(grid=grid, cmap=True) +fig.basemap(frame="g1") +fig.plot(data=circle, pen="2p,darkorange") + +fig.shift_origin(xshift="+w+2c") + +# Plot the masked elevation grid +fig.basemap(region=region, projection="M12c", frame=True) +fig.grdimage(grid=grid_mask, cmap=True) +fig.basemap(frame="g1") +fig.plot(data=circle, pen="2p,darkorange") + +fig.colorbar(frame=True) +fig.show() + +# sphinx_gallery_thumbnail_number = 1 From 4584da72b5f8c79e1e1c1b432b96539b727c6fae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yvonne=20Fr=C3=B6hlich?= Date: Thu, 9 Apr 2026 21:00:16 +0200 Subject: [PATCH 04/12] Remove execution permission --- examples/gallery/images/grdmask.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100755 => 100644 examples/gallery/images/grdmask.py diff --git a/examples/gallery/images/grdmask.py b/examples/gallery/images/grdmask.py old mode 100755 new mode 100644 From 2b8a6321d75ec4cc468b8fe5136ea48c3da56a8b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yvonne=20Fr=C3=B6hlich?= Date: Fri, 10 Apr 2026 11:20:15 +0200 Subject: [PATCH 05/12] Select US staat --- examples/gallery/images/grdmask.py | 46 ++++++++++++++++++++++++------ 1 file changed, 37 insertions(+), 9 deletions(-) diff --git a/examples/gallery/images/grdmask.py b/examples/gallery/images/grdmask.py index c65fffb6408..9a7c9fc5340 100644 --- a/examples/gallery/images/grdmask.py +++ b/examples/gallery/images/grdmask.py @@ -68,26 +68,56 @@ fig.show() +# %% +# US staat based on GeoPandas polygon geometry +# -------------------------------------------- + +region = [-126, -66, 25, 49] + +provider = "https://naciscdn.org/naturalearth" +states = geopandas.read_file( + f"{provider}/50m/cultural/ne_50m_admin_1_states_provinces.zip" +) +wyoming = states[states["name"] == "Missouri"] + +grid = pygmt.datasets.load_earth_relief(region=region, resolution="01m") +mask = pygmt.grdmask(region=region, data=wyoming, spacing="01m", inside="NaN") +mask_lonlat = mask.rename(new_name_or_name_dict={"x": "lon", "y": "lat"}) +grid_mask = grid * mask_lonlat + + +fig = pygmt.Figure() +pygmt.makecpt(cmap="oleron", series=[-2000, 2000]) + +# Plot the elevation grid +fig.basemap("L-96/35/33/41/12c", region=region, frame=True) +fig.grdimage(grid=grid, cmap=True) +fig.plot(data=wyoming, pen="1p,darkorange") + +fig.shift_origin(xshift="+w+1c") + +# Plot the masked elevation grid +fig.basemap("L-96/35/33/41/12c", region=region, frame=True) +fig.grdimage(grid=grid_mask, cmap=True) +fig.plot(data=wyoming, pen="1p,darkorange") + +fig.colorbar(frame=True) +fig.show() + + # %% # Circle based on GeoPandas polygon geometry # ------------------------------------------ -# Define a study region region = [125, 135, 25, 36] # Create a point and buffer it point = geopandas.GeoSeries([Point(126.5, 33.5)]) circle = point.buffer(0.6) # 10 is the radius -# Download elevation grid grid = pygmt.datasets.load_earth_relief(region=region, resolution="30s") - -# Create a grid mask based on the two polygons defined above, set all values -# outside the polygons to NaN mask = pygmt.grdmask(region=region, data=circle, spacing="30s", outside="NaN") mask_lonlat = mask.rename(new_name_or_name_dict={"x": "lon", "y": "lat"}) - -# Apply the grid mask to the downloaded elevation grid by multiplying the two grids grid_mask = grid * mask_lonlat @@ -97,7 +127,6 @@ # Plot the elevation grid fig.basemap(region=region, projection="M12c", frame=True) fig.grdimage(grid=grid, cmap=True) -fig.basemap(frame="g1") fig.plot(data=circle, pen="2p,darkorange") fig.shift_origin(xshift="+w+2c") @@ -105,7 +134,6 @@ # Plot the masked elevation grid fig.basemap(region=region, projection="M12c", frame=True) fig.grdimage(grid=grid_mask, cmap=True) -fig.basemap(frame="g1") fig.plot(data=circle, pen="2p,darkorange") fig.colorbar(frame=True) From b936f5d5d661fe4b79c9e205fb079920391390d7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yvonne=20Fr=C3=B6hlich?= Date: Fri, 10 Apr 2026 12:58:08 +0200 Subject: [PATCH 06/12] Zoom in mask grid --- examples/gallery/images/grdmask.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/gallery/images/grdmask.py b/examples/gallery/images/grdmask.py index 9a7c9fc5340..779cccafc1b 100644 --- a/examples/gallery/images/grdmask.py +++ b/examples/gallery/images/grdmask.py @@ -132,7 +132,7 @@ fig.shift_origin(xshift="+w+2c") # Plot the masked elevation grid -fig.basemap(region=region, projection="M12c", frame=True) +fig.basemap(region=[125.5, 127.5, 32.5, 34.5], projection="M12c", frame=True) fig.grdimage(grid=grid_mask, cmap=True) fig.plot(data=circle, pen="2p,darkorange") From 2c1c2644075e10aa412744e44d34bdd2f2373fd2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yvonne=20Fr=C3=B6hlich?= Date: Fri, 10 Apr 2026 13:04:40 +0200 Subject: [PATCH 07/12] Try inside and outside --- examples/gallery/images/grdmask.py | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/examples/gallery/images/grdmask.py b/examples/gallery/images/grdmask.py index 779cccafc1b..05fd7337bc9 100644 --- a/examples/gallery/images/grdmask.py +++ b/examples/gallery/images/grdmask.py @@ -39,12 +39,15 @@ # Download elevation grid grid = pygmt.datasets.load_earth_relief(region=region, resolution="30s") -# Create a grid mask based on the two polygons defined above, set all values -# outside the polygons to NaN -mask = pygmt.grdmask(region=region, data=polygon, spacing="30s", outside="NaN") +# Create a grid mask based on the two polygons defined above +# Set all values outside the polygons to NaN +mask_out = pygmt.grdmask(region=region, data=polygon, spacing="30s", outside="NaN") +# Set all values inside the polygons to NaN +mask_in = pygmt.grdmask(region=region, data=polygon, spacing="30s", inside="NaN") # Apply the grid mask to the downloaded elevation grid by multiplying the two grids -grid_mask = grid * mask +grid_mask_out = grid * mask_out +grid_mask_in = grid * mask_in fig = pygmt.Figure() @@ -53,15 +56,20 @@ # Plot the elevation grid fig.basemap(region=region, projection="M12c", frame=True) fig.grdimage(grid=grid, cmap=True) -fig.basemap(frame="g1") fig.plot(data=polygon, pen="2p,darkorange") fig.shift_origin(xshift="+w+2c") -# Plot the masked elevation grid +# Plot the masked elevation grid outside fig.basemap(region=region, projection="M12c", frame=True) -fig.grdimage(grid=grid_mask, cmap=True) -fig.basemap(frame="g1") +fig.grdimage(grid=grid_mask_out, cmap=True) +fig.plot(data=polygon, pen="2p,darkorange") + +fig.shift_origin(xshift="+w+2c") + +# Plot the masked elevation grid inside +fig.basemap(region=region, projection="M12c", frame=True) +fig.grdimage(grid=grid_mask_in, cmap=True) fig.plot(data=polygon, pen="2p,darkorange") fig.colorbar(frame=True) From a62befc48f5e9c345b89afcb56239d5a087c1536 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yvonne=20Fr=C3=B6hlich?= Date: Fri, 10 Apr 2026 13:39:39 +0200 Subject: [PATCH 08/12] Adjust area for US staat, fix typo, improve code --- examples/gallery/images/grdmask.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/examples/gallery/images/grdmask.py b/examples/gallery/images/grdmask.py index 05fd7337bc9..5ecaf46b5f4 100644 --- a/examples/gallery/images/grdmask.py +++ b/examples/gallery/images/grdmask.py @@ -86,10 +86,10 @@ states = geopandas.read_file( f"{provider}/50m/cultural/ne_50m_admin_1_states_provinces.zip" ) -wyoming = states[states["name"] == "Missouri"] +missouri = states[states["name"] == "Missouri"] grid = pygmt.datasets.load_earth_relief(region=region, resolution="01m") -mask = pygmt.grdmask(region=region, data=wyoming, spacing="01m", inside="NaN") +mask = pygmt.grdmask(region=region, data=missouri, spacing="01m", outside="NaN") mask_lonlat = mask.rename(new_name_or_name_dict={"x": "lon", "y": "lat"}) grid_mask = grid * mask_lonlat @@ -98,16 +98,17 @@ pygmt.makecpt(cmap="oleron", series=[-2000, 2000]) # Plot the elevation grid -fig.basemap("L-96/35/33/41/12c", region=region, frame=True) +fig.basemap(projection="L-96/35/33/41/12c", region=region, frame=True) fig.grdimage(grid=grid, cmap=True) -fig.plot(data=wyoming, pen="1p,darkorange") +fig.plot(data=missouri, pen="1p,darkorange") fig.shift_origin(xshift="+w+1c") # Plot the masked elevation grid -fig.basemap("L-96/35/33/41/12c", region=region, frame=True) +#fig.basemap(projection="L-96/35/33/41/12c", region=region, frame=True) +fig.basemap(projection="M10c", region=[-96.5, -88.5, 35.8, 41], frame=True) fig.grdimage(grid=grid_mask, cmap=True) -fig.plot(data=wyoming, pen="1p,darkorange") +fig.plot(data=missouri, pen="1p,darkorange") fig.colorbar(frame=True) fig.show() From 9cb4b75912a170748a286297a680d70fd6f6549a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yvonne=20Fr=C3=B6hlich?= Date: Fri, 10 Apr 2026 13:42:39 +0200 Subject: [PATCH 09/12] Fix style --- examples/gallery/images/grdmask.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/gallery/images/grdmask.py b/examples/gallery/images/grdmask.py index 5ecaf46b5f4..63cf542ec92 100644 --- a/examples/gallery/images/grdmask.py +++ b/examples/gallery/images/grdmask.py @@ -105,7 +105,7 @@ fig.shift_origin(xshift="+w+1c") # Plot the masked elevation grid -#fig.basemap(projection="L-96/35/33/41/12c", region=region, frame=True) +# fig.basemap(projection="L-96/35/33/41/12c", region=region, frame=True) fig.basemap(projection="M10c", region=[-96.5, -88.5, 35.8, 41], frame=True) fig.grdimage(grid=grid_mask, cmap=True) fig.plot(data=missouri, pen="1p,darkorange") From c0ef6115371367418128fb3886da3c791b8821d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yvonne=20Fr=C3=B6hlich?= Date: Fri, 10 Apr 2026 14:04:46 +0200 Subject: [PATCH 10/12] Set outside parameter correctly --- examples/gallery/images/grdmask.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/examples/gallery/images/grdmask.py b/examples/gallery/images/grdmask.py index 63cf542ec92..c9ae63dc323 100644 --- a/examples/gallery/images/grdmask.py +++ b/examples/gallery/images/grdmask.py @@ -40,10 +40,13 @@ grid = pygmt.datasets.load_earth_relief(region=region, resolution="30s") # Create a grid mask based on the two polygons defined above -# Set all values outside the polygons to NaN +# Set all grid nodes outside the polygons to NaN mask_out = pygmt.grdmask(region=region, data=polygon, spacing="30s", outside="NaN") -# Set all values inside the polygons to NaN -mask_in = pygmt.grdmask(region=region, data=polygon, spacing="30s", inside="NaN") +# Set all grid nodes inside the polygons to NaN +# Set the outside parameter to a value larger 0 to keep the nodes outside unchanged +mask_in = pygmt.grdmask( + region=region, data=polygon, spacing="30s", inside="NaN", outside=1 +) # Apply the grid mask to the downloaded elevation grid by multiplying the two grids grid_mask_out = grid * mask_out From 6326c199e9961d9e4f8c30b321ee6dabd0cca817 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yvonne=20Fr=C3=B6hlich?= Date: Fri, 10 Apr 2026 14:18:45 +0200 Subject: [PATCH 11/12] Fix typo --- examples/gallery/images/grdmask.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/gallery/images/grdmask.py b/examples/gallery/images/grdmask.py index c9ae63dc323..f89b825bbb7 100644 --- a/examples/gallery/images/grdmask.py +++ b/examples/gallery/images/grdmask.py @@ -80,8 +80,8 @@ # %% -# US staat based on GeoPandas polygon geometry -# -------------------------------------------- +# US state Missouri based on GeoPandas polygon geometry +# ----------------------------------------------------- region = [-126, -66, 25, 49] From a4d8f5ebd2d493832e79d4fa22b3e6dcc61f63aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yvonne=20Fr=C3=B6hlich?= Date: Fri, 10 Apr 2026 18:22:47 +0200 Subject: [PATCH 12/12] Adjust title - circles are no polygons --- examples/gallery/images/grdmask.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/gallery/images/grdmask.py b/examples/gallery/images/grdmask.py index f89b825bbb7..dbd955781ca 100644 --- a/examples/gallery/images/grdmask.py +++ b/examples/gallery/images/grdmask.py @@ -1,6 +1,6 @@ """ -Create grid mask from polygons -============================== +Create grid masks from geospatial shapes +======================================== :func:`pygmt.grdmask`. :func:`pygmt.grdlandmask` and gallery example https://www.pygmt.org/latest/gallery/images/grdlandmask.html.