From f256ab4d8982382602fcd0378589c290597e7f2d Mon Sep 17 00:00:00 2001 From: Kshitij Aucharmal Date: Wed, 1 Oct 2025 10:58:16 +0530 Subject: [PATCH 1/5] Add polygon-based subregion support to grdcut with gallery example --- examples/gallery/images/grdpolygon.py | 63 +++++++++++++++++++++++++++ pygmt/src/grdcut.py | 20 ++++++++- 2 files changed, 81 insertions(+), 2 deletions(-) create mode 100644 examples/gallery/images/grdpolygon.py diff --git a/examples/gallery/images/grdpolygon.py b/examples/gallery/images/grdpolygon.py new file mode 100644 index 00000000000..d702848538b --- /dev/null +++ b/examples/gallery/images/grdpolygon.py @@ -0,0 +1,63 @@ +""" +Clipping grid to a complex polygon +================================== + +The :func:`pygmt.grdcut` function allows you to extract a subregion from a +grid. In this example we use a complex polygon (GeoDataFrame or GMT ASCII file) +to crop the grid to a region of interest. +""" + +# %% +import geopandas as gpd +import pygmt +from shapely.geometry import Polygon + +fig = pygmt.Figure() + +# Define region of interest around Iceland +region = [-28, -10, 62, 68] + +# Load sample grid (3 arc-minutes global relief) in target area +grid = pygmt.datasets.load_earth_relief(resolution="03m", region=region) + +# Plot original grid +fig.basemap( + region=region, + projection="M12c", + frame=["WSne+toriginal grid", "xa5f1", "ya2f1"], +) +fig.grdimage(grid=grid, cmap="oleron") + +# Shift plot origin of the second map by "width of the first map + 0.5 cm" +# in x-direction +fig.shift_origin(xshift="w+0.5c") + +# Create a more complex polygon (irregular shape) around a smaller ROI +complex_poly = Polygon( + [ + (-26, 63), + (-23, 63.5), + (-20, 64), + (-18, 65), + (-19, 66), + (-22, 66.5), + (-25, 66), + (-27, 65), + (-26, 63), + ] +) +gdf = gpd.GeoDataFrame({"geometry": [complex_poly]}, crs="EPSG:4326") + +# Crop the grid using the complex polygon +cropped_grid = pygmt.grdcut(grid=grid, polygon=gdf) + +# Plot cropped grid +fig.basemap( + region=region, + projection="M12c", + frame=["wSne+tpolygon cropped grid", "xa5f1", "ya2f1"], +) +fig.grdimage(grid=cropped_grid, cmap="oleron") +fig.colorbar(frame=["x+lElevation", "y+lm"], position="JMR+o0.5c/0c+w8c") + +fig.show() diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index cf198d9dc74..340cf9ff0fc 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -2,6 +2,7 @@ grdcut - Extract subregion from a grid or image or a slice from a cube. """ +from contextlib import ExitStack from collections.abc import Sequence from typing import Literal @@ -15,13 +16,14 @@ data_kind, fmt_docstring, use_alias, + tempfile_from_geojson, ) __doctest_skip__ = ["grdcut"] @fmt_docstring -@use_alias(N="extend", S="circ_subregion", Z="z_subregion", f="coltypes") +@use_alias(N="extend", S="circ_subregion", Z="z_subregion", f="coltypes", F="polygon") def grdcut( grid: PathLike | xr.DataArray, kind: Literal["grid", "image"] = "grid", @@ -86,6 +88,12 @@ def grdcut( NaNs, append **+N** to strip off such columns before (optionally) considering the range of the core subset for further reduction of the area. + polygon : str, geopandas.GeoDataFrame, or shapely.geometry + Extract a subregion using a polygon. Can be either: + - A GMT ASCII polygon file (`.gmt`) + - A geopandas.GeoDataFrame (must have CRS EPSG:4326) + - A shapely.geometry.Polygon or MultiPolygon + The polygon can have holes or multiple rings. {verbose} {coltypes} @@ -130,12 +138,20 @@ def grdcut( ) aliasdict.merge(kwargs) - with Session() as lib: + with Session() as lib, ExitStack() as stack: with ( lib.virtualfile_in(check_kind="raster", data=grid) as vingrd, lib.virtualfile_out(kind=outkind, fname=outgrid) as voutgrd, ): aliasdict["G"] = voutgrd + + if "F" in kwargs and kwargs["F"] is not None: + polygon_input = kwargs["F"] + if not isinstance(polygon_input, (str, bytes)): + tmpfile = stack.enter_context(tempfile_from_geojson(polygon_input)) + aliasdict["F"] = tmpfile + + lib.call_module( module="grdcut", args=build_arg_list(aliasdict, infile=vingrd) ) From bac59da9488e8962e3b0d6355d08cbcc513c9e97 Mon Sep 17 00:00:00 2001 From: Kshitij Aucharmal Date: Wed, 1 Oct 2025 11:24:50 +0530 Subject: [PATCH 2/5] Add crop and invert boolean support for polygon input --- examples/gallery/images/grdpolygon.py | 51 ++++++++++++++++++--------- pygmt/src/grdcut.py | 18 +++++++--- 2 files changed, 48 insertions(+), 21 deletions(-) diff --git a/examples/gallery/images/grdpolygon.py b/examples/gallery/images/grdpolygon.py index d702848538b..0567123d155 100644 --- a/examples/gallery/images/grdpolygon.py +++ b/examples/gallery/images/grdpolygon.py @@ -4,7 +4,7 @@ The :func:`pygmt.grdcut` function allows you to extract a subregion from a grid. In this example we use a complex polygon (GeoDataFrame or GMT ASCII file) -to crop the grid to a region of interest. +to crop the grid, invert it, or do both. """ # %% @@ -20,18 +20,6 @@ # Load sample grid (3 arc-minutes global relief) in target area grid = pygmt.datasets.load_earth_relief(resolution="03m", region=region) -# Plot original grid -fig.basemap( - region=region, - projection="M12c", - frame=["WSne+toriginal grid", "xa5f1", "ya2f1"], -) -fig.grdimage(grid=grid, cmap="oleron") - -# Shift plot origin of the second map by "width of the first map + 0.5 cm" -# in x-direction -fig.shift_origin(xshift="w+0.5c") - # Create a more complex polygon (irregular shape) around a smaller ROI complex_poly = Polygon( [ @@ -48,16 +36,45 @@ ) gdf = gpd.GeoDataFrame({"geometry": [complex_poly]}, crs="EPSG:4326") -# Crop the grid using the complex polygon -cropped_grid = pygmt.grdcut(grid=grid, polygon=gdf) +# Original grid +fig.basemap( + region=region, + projection="M12c", + frame=["WSne+toriginal grid", "xa5f1", "ya2f1"], +) +fig.grdimage(grid=grid, cmap="oleron") -# Plot cropped grid +# Cropped grid +fig.shift_origin(xshift="w+0.5c") +cropped_grid = pygmt.grdcut(grid=grid, polygon=gdf, crop=True) fig.basemap( region=region, projection="M12c", - frame=["wSne+tpolygon cropped grid", "xa5f1", "ya2f1"], + frame=["WSne+tcropped", "xa5f1", "ya2f1"], ) fig.grdimage(grid=cropped_grid, cmap="oleron") + +# Inverted grid +fig.shift_origin(xshift="w+0.5c") +inverted_grid = pygmt.grdcut(grid=grid, polygon=gdf, invert=True) +fig.basemap( + region=region, + projection="M12c", + frame=["WSne+tinverted", "xa5f1", "ya2f1"], +) +fig.grdimage(grid=inverted_grid, cmap="oleron") + +# Cropped + inverted grid +fig.shift_origin(xshift="w+0.5c") +cropped_inverted_grid = pygmt.grdcut(grid=grid, polygon=gdf, crop=True, invert=True) +fig.basemap( + region=region, + projection="M12c", + frame=["WSne+tcropped and inverted", "xa5f1", "ya2f1"], +) +fig.grdimage(grid=cropped_inverted_grid, cmap="oleron") + +# Shared colorbar fig.colorbar(frame=["x+lElevation", "y+lm"], position="JMR+o0.5c/0c+w8c") fig.show() diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index 340cf9ff0fc..bed82e0286c 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -2,8 +2,8 @@ grdcut - Extract subregion from a grid or image or a slice from a cube. """ -from contextlib import ExitStack from collections.abc import Sequence +from contextlib import ExitStack from typing import Literal import xarray as xr @@ -15,8 +15,8 @@ build_arg_list, data_kind, fmt_docstring, - use_alias, tempfile_from_geojson, + use_alias, ) __doctest_skip__ = ["grdcut"] @@ -32,6 +32,8 @@ def grdcut( region: Sequence[float | str] | str | None = None, verbose: Literal["quiet", "error", "warning", "timing", "info", "compat", "debug"] | bool = False, + crop: bool = False, + invert: bool = False, **kwargs, ) -> xr.DataArray | None: r""" @@ -94,6 +96,11 @@ def grdcut( - A geopandas.GeoDataFrame (must have CRS EPSG:4326) - A shapely.geometry.Polygon or MultiPolygon The polygon can have holes or multiple rings. + Optional modifiers: + - crop : bool + If True, crop the output grid region to the bounding box of the polygon. + - invert : bool + If True, invert the selection, setting all nodes inside the polygon to NaN. {verbose} {coltypes} @@ -147,10 +154,13 @@ def grdcut( if "F" in kwargs and kwargs["F"] is not None: polygon_input = kwargs["F"] + modifiers = ("+c" * crop) + ("+i" * invert) + if not isinstance(polygon_input, (str, bytes)): tmpfile = stack.enter_context(tempfile_from_geojson(polygon_input)) - aliasdict["F"] = tmpfile - + aliasdict["F"] = tmpfile + modifiers + else: + aliasdict["F"] = str(polygon_input) + modifiers lib.call_module( module="grdcut", args=build_arg_list(aliasdict, infile=vingrd) From f6d4179c5d93f8674d08b3d12bd30a69ab00fb0f Mon Sep 17 00:00:00 2001 From: Kshitij Aucharmal Date: Thu, 2 Oct 2025 10:55:08 +0530 Subject: [PATCH 3/5] grdcut: PathLike for polygon input and improve example - Checking for PathLike instead of str, bytes - Inverted logic and comments to make it clearer - Update gallery example: - removed "complex" from polygon desc. - improved docs - Switch CRS to OGC:CRS84 --- examples/gallery/images/grdpolygon.py | 16 ++++++++-------- pygmt/src/grdcut.py | 11 ++++++----- 2 files changed, 14 insertions(+), 13 deletions(-) diff --git a/examples/gallery/images/grdpolygon.py b/examples/gallery/images/grdpolygon.py index 0567123d155..2e2eaefb321 100644 --- a/examples/gallery/images/grdpolygon.py +++ b/examples/gallery/images/grdpolygon.py @@ -1,10 +1,10 @@ """ -Clipping grid to a complex polygon -================================== +Clipping grid to a polygon +========================== The :func:`pygmt.grdcut` function allows you to extract a subregion from a -grid. In this example we use a complex polygon (GeoDataFrame or GMT ASCII file) -to crop the grid, invert it, or do both. +grid. In this example we use a :class:`geopandas.GeoDataFrame` +to crop the grid to the polygon's inner or outer region, invert it, or do both. """ # %% @@ -15,13 +15,13 @@ fig = pygmt.Figure() # Define region of interest around Iceland -region = [-28, -10, 62, 68] +region = [-28, -10, 62, 68] # xmin, xmax, ymin, ymax # Load sample grid (3 arc-minutes global relief) in target area grid = pygmt.datasets.load_earth_relief(resolution="03m", region=region) -# Create a more complex polygon (irregular shape) around a smaller ROI -complex_poly = Polygon( +# Create a more polygon (irregular shape) around a smaller ROI +poly = Polygon( [ (-26, 63), (-23, 63.5), @@ -34,7 +34,7 @@ (-26, 63), ] ) -gdf = gpd.GeoDataFrame({"geometry": [complex_poly]}, crs="EPSG:4326") +gdf = gpd.GeoDataFrame({"geometry": [poly]}, crs="OGC:CRS84") # Original grid fig.basemap( diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index bed82e0286c..1b270f55a13 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -153,14 +153,15 @@ def grdcut( aliasdict["G"] = voutgrd if "F" in kwargs and kwargs["F"] is not None: - polygon_input = kwargs["F"] modifiers = ("+c" * crop) + ("+i" * invert) - if not isinstance(polygon_input, (str, bytes)): - tmpfile = stack.enter_context(tempfile_from_geojson(polygon_input)) - aliasdict["F"] = tmpfile + modifiers + # if file path provided + if isinstance(kwargs["F"], PathLike): + aliasdict["F"] = str(kwargs["F"]) + modifiers + # assuming its geojson else: - aliasdict["F"] = str(polygon_input) + modifiers + tmpfile = stack.enter_context(tempfile_from_geojson(kwargs["F"])) + aliasdict["F"] = tmpfile + modifiers lib.call_module( module="grdcut", args=build_arg_list(aliasdict, infile=vingrd) From bf78b8b80980191c713127ceb45793f708f1e8f0 Mon Sep 17 00:00:00 2001 From: Kshitij Aucharmal Date: Thu, 2 Oct 2025 11:24:53 +0530 Subject: [PATCH 4/5] Add polygon input support tests for grdcut - Validated polygon input types, raising GMTValueError for unsupported types. - Unit tests added for: - Shapely Polygon input - GeoDataFrame polygon input - GMT ASCII polygon file input - Crop (+c) and invert (+i) modifiers - Invalid polygon input type to raise error --- pygmt/src/grdcut.py | 15 +++++++-- pygmt/tests/test_grdcut.py | 68 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 80 insertions(+), 3 deletions(-) diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index 1b270f55a13..3f632ffbc17 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -153,15 +153,24 @@ def grdcut( aliasdict["G"] = voutgrd if "F" in kwargs and kwargs["F"] is not None: + polygon_input = kwargs["F"] modifiers = ("+c" * crop) + ("+i" * invert) # if file path provided - if isinstance(kwargs["F"], PathLike): + if isinstance(polygon_input, PathLike): aliasdict["F"] = str(kwargs["F"]) + modifiers # assuming its geojson - else: - tmpfile = stack.enter_context(tempfile_from_geojson(kwargs["F"])) + elif hasattr(polygon_input, "__geo_interface__"): + tmpfile = stack.enter_context(tempfile_from_geojson(polygon_input)) aliasdict["F"] = tmpfile + modifiers + else: + polygon_type = type(polygon_input).__name__ + err_msg = "Invalid polygon type" + raise GMTValueError( + err_msg, + polygon_type, + "Must be a PathLike, GeoDataFrame, or Shapely geometry.", + ) lib.call_module( module="grdcut", args=build_arg_list(aliasdict, infile=vingrd) diff --git a/pygmt/tests/test_grdcut.py b/pygmt/tests/test_grdcut.py index 1cd54caf4b6..818eee4037d 100644 --- a/pygmt/tests/test_grdcut.py +++ b/pygmt/tests/test_grdcut.py @@ -2,6 +2,7 @@ Test pygmt.grdcut. """ +import geopandas as gpd import numpy as np import pytest import xarray as xr @@ -9,6 +10,18 @@ from pygmt.exceptions import GMTTypeError, GMTValueError from pygmt.helpers import GMTTempFile from pygmt.helpers.testing import load_static_earth_relief +from shapely.geometry import Polygon + +# Reusable polygon geometry +poly = Polygon( + [ + (-52.5, -19.5), + (-51.5, -19.0), + (-50.5, -18.5), + (-51.5, -18.0), + (-52.5, -19.5), + ] +) @pytest.fixture(scope="module", name="grid") @@ -78,3 +91,58 @@ def test_grdcut_invalid_kind(grid, region): """ with pytest.raises(GMTValueError): grdcut(grid, kind="invalid", region=region) + + +def test_grdcut_with_shapely_polygon(grid): + """ + Grdcut should accept a shapely.geometry.Polygon as polygon input. + """ + outgrid = grdcut(grid=grid, polygon=poly) + assert outgrid is not None + assert outgrid.rio.crs is None or str(outgrid.rio.crs) == "" # CRS not enforced + assert outgrid.size > 0 + + +def test_grdcut_with_geodataframe_polygon(grid): + """ + Grdcut should accept a geopandas.GeoDataFrame as polygon input. + """ + gdf = gpd.GeoDataFrame({"geometry": [poly]}, crs="OGC:CRS84") + outgrid = grdcut(grid=grid, polygon=gdf) + assert outgrid is not None + assert outgrid.size > 0 + + +def test_grdcut_with_polygon_file(grid, tmp_path): + """ + Grdcut should accept a GMT ASCII polygon file as polygon input. + """ + gmtfile = tmp_path / "poly.gmt" + gmtfile.write_text( + ">\n-52.5 -19.5\n-51.5 -19.0\n-50.5 -18.5\n-51.5 -18.0\n-52.5 -19.5\n" + ) + + outgrid = grdcut(grid=grid, polygon=gmtfile) + assert outgrid is not None + assert outgrid.size > 0 + + +@pytest.mark.parametrize( + ("crop", "invert"), [(True, False), (False, True), (True, True)] +) +def test_grdcut_polygon_with_crop_and_invert(grid, crop, invert): + """ + Grdcut should support crop (+c) and invert (+i) modifiers with polygon input. + """ + gdf = gpd.GeoDataFrame({"geometry": [poly]}, crs="OGC:CRS84") + outgrid = grdcut(grid=grid, polygon=gdf, crop=crop, invert=invert) + assert outgrid is not None + assert outgrid.size > 0 + + +def test_grdcut_polygon_invalid_input(grid): + """ + Grdcut should raise an error for invalid polygon input type. + """ + with pytest.raises(GMTValueError): + grdcut(grid=grid, polygon=12345) From 4e4560723ae515aa248fadd47da2673228e99a82 Mon Sep 17 00:00:00 2001 From: Kshitij Aucharmal Date: Thu, 2 Oct 2025 13:27:51 +0530 Subject: [PATCH 5/5] Add polygon fixture and more precise grdcut unit tests - Reusable polygon fixture for grdcut tests - Improve assertions to check pixel counts and presence of NaNs outside polygons - Add a test for MultiPolygon support --- pygmt/src/grdcut.py | 2 +- pygmt/tests/test_grdcut.py | 90 +++++++++++++++++++++++++++++--------- 2 files changed, 70 insertions(+), 22 deletions(-) diff --git a/pygmt/src/grdcut.py b/pygmt/src/grdcut.py index 3f632ffbc17..aea856bfabf 100644 --- a/pygmt/src/grdcut.py +++ b/pygmt/src/grdcut.py @@ -93,7 +93,7 @@ def grdcut( polygon : str, geopandas.GeoDataFrame, or shapely.geometry Extract a subregion using a polygon. Can be either: - A GMT ASCII polygon file (`.gmt`) - - A geopandas.GeoDataFrame (must have CRS EPSG:4326) + - A geopandas.GeoDataFrame - A shapely.geometry.Polygon or MultiPolygon The polygon can have holes or multiple rings. Optional modifiers: diff --git a/pygmt/tests/test_grdcut.py b/pygmt/tests/test_grdcut.py index 818eee4037d..3d9632e3ed6 100644 --- a/pygmt/tests/test_grdcut.py +++ b/pygmt/tests/test_grdcut.py @@ -10,18 +10,23 @@ from pygmt.exceptions import GMTTypeError, GMTValueError from pygmt.helpers import GMTTempFile from pygmt.helpers.testing import load_static_earth_relief -from shapely.geometry import Polygon - -# Reusable polygon geometry -poly = Polygon( - [ - (-52.5, -19.5), - (-51.5, -19.0), - (-50.5, -18.5), - (-51.5, -18.0), - (-52.5, -19.5), - ] -) +from shapely.geometry import MultiPolygon, Polygon + + +@pytest.fixture(scope="module", name="polygon") +def fixture_polygon(): + """ + Provide a reusable polygon geometry for grdcut tests. + """ + return Polygon( + [ + (-52.5, -19.5), + (-51.5, -19.0), + (-50.5, -18.5), + (-51.5, -18.0), + (-52.5, -19.5), + ] + ) @pytest.fixture(scope="module", name="grid") @@ -93,25 +98,45 @@ def test_grdcut_invalid_kind(grid, region): grdcut(grid, kind="invalid", region=region) -def test_grdcut_with_shapely_polygon(grid): +def test_grdcut_with_shapely_polygon(grid, polygon): """ - Grdcut should accept a shapely.geometry.Polygon as polygon input. + Grdcut should accept a shapely.geometry.Polygon as polygon input + and produce a grid of expected size with some valid data. """ - outgrid = grdcut(grid=grid, polygon=poly) + outgrid = grdcut(grid=grid, polygon=polygon) assert outgrid is not None - assert outgrid.rio.crs is None or str(outgrid.rio.crs) == "" # CRS not enforced + + # Check size assert outgrid.size > 0 + # There should be some non-NaN pixels inside the polygon + assert np.any(~np.isnan(outgrid.values)) + + # Optionally, corners outside the polygon should be NaN + corner_values = [ + outgrid.values[0, 0], + outgrid.values[0, -1], + outgrid.values[-1, 0], + outgrid.values[-1, -1], + ] + assert any(np.isnan(val) for val in corner_values) + -def test_grdcut_with_geodataframe_polygon(grid): +def test_grdcut_with_geodataframe_polygon(grid, polygon): """ - Grdcut should accept a geopandas.GeoDataFrame as polygon input. + Grdcut should accept a geopandas.GeoDataFrame as polygon input + and produce a grid of expected size with some valid data. """ - gdf = gpd.GeoDataFrame({"geometry": [poly]}, crs="OGC:CRS84") + gdf = gpd.GeoDataFrame({"geometry": [polygon]}, crs="OGC:CRS84") outgrid = grdcut(grid=grid, polygon=gdf) assert outgrid is not None + + # Check size assert outgrid.size > 0 + # Some pixels should remain non-NaN + assert np.any(~np.isnan(outgrid.values)) + def test_grdcut_with_polygon_file(grid, tmp_path): """ @@ -130,15 +155,38 @@ def test_grdcut_with_polygon_file(grid, tmp_path): @pytest.mark.parametrize( ("crop", "invert"), [(True, False), (False, True), (True, True)] ) -def test_grdcut_polygon_with_crop_and_invert(grid, crop, invert): +def test_grdcut_polygon_with_crop_and_invert(grid, polygon, crop, invert): """ Grdcut should support crop (+c) and invert (+i) modifiers with polygon input. """ - gdf = gpd.GeoDataFrame({"geometry": [poly]}, crs="OGC:CRS84") + gdf = gpd.GeoDataFrame({"geometry": [polygon]}, crs="OGC:CRS84") outgrid = grdcut(grid=grid, polygon=gdf, crop=crop, invert=invert) assert outgrid is not None assert outgrid.size > 0 + if invert: + assert np.isnan(outgrid.data).any() + + assert np.isfinite(outgrid.data).any() + + assert np.count_nonzero(np.isfinite(outgrid.data)) < grid.size + + +def test_grdcut_with_multipolygon(grid): + """ + Grdcut should accept a shapely.geometry.MultiPolygon as polygon input. + """ + # Create two simple polygons and combine into a MultiPolygon + poly1 = Polygon([(-52.5, -19.5), (-51.5, -19.0), (-50.5, -18.5), (-52.5, -19.5)]) + poly2 = Polygon([(-51.0, -18.5), (-50.5, -18.0), (-50.0, -18.5), (-51.0, -18.5)]) + multipoly = MultiPolygon([poly1, poly2]) + + outgrid = grdcut(grid=grid, polygon=multipoly) + assert outgrid is not None + # Ensure some grid pixels exist + assert outgrid.size > 0 + assert np.any(np.isnan(outgrid.values)) + def test_grdcut_polygon_invalid_input(grid): """