Support polygon-based subsetting in grdcut -F with GeoDataFrame #4135

kshitijaucharmal · 2025-10-01T06:07:13Z

Description of proposed changes

This PR enhances the pygmt.grdcut function to support polygon-based subsetting of grids using the -F option. Users can now provide polygons in three different formats:

A geopandas.GeoDataFrame
A shapely.geometry.Polygon
A GMT ASCII polygon file

Additionally, two new boolean parameters, crop and invert, have been added to control how the polygon is applied:

crop=True crops the output grid region to the bounding box of the polygon.
invert=True inverts the selection, setting all nodes inside the polygon to NaN.

I have added both of these as parameters to grdcut, but they only apply to polygon mode. Please let me know if this is undesirable, and how would I implement this otherwise

A gallery example has been added to demonstrate all three polygon input methods and the effects of the crop and invert options, showing original, cropped, inverted, and cropped+inverted grids. This improves PyGMT’s functionality for extracting regions of interest from raster grids and makes the workflow more Pythonic and flexible.

Fixes #1529 (implements rather)

Preview: https://pygmt-dev--4135.org.readthedocs.build/en/4135/gallery/images/grdpolygon.html

Users can crop grids using polygons without manually creating ASCII files.
Example figure showing original grid and cropped and inverted grids

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 wrapping a new module, open a 'Wrap new GMT module' issue and submit reasonably-sized PRs.
If adding new functionality, add an example to docstrings or tutorials.

kshitijaucharmal · 2025-10-01T11:36:03Z

@yvonnefroehlich I have added both of the +c and +i options as parameters to grdcut, but they only apply to polygon mode. Please let me know if this is undesirable, and how would I implement this otherwise

Also, the label is hacktoberfest-accepted (if you approve of this PR of course) :)

weiji14

Thanks @kshitijaucharmal! Just a few comments for now. If you're able to add some unit tests to test_grdcut.py, that would be fantastic too!

examples/gallery/images/grdpolygon.py

pygmt/src/grdcut.py

- 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

- 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

kshitijaucharmal · 2025-10-02T05:59:02Z

@weiji14 I have resolved all comments and added the unit tests as well, let me know if there are any more corrections! :)

weiji14

Thanks for adding the unit tests! Here are some more comments, and apologies for the failing CI tests, we have some issues with DVC running on forks that we'll probably try to sort out next week.

weiji14 · 2025-10-02T06:20:07Z

pygmt/tests/test_grdcut.py

+# 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),
+    ]
+)


This should be turned into a pytest fixture, see examples a few lines below.

Will convert this into a fixture

weiji14 · 2025-10-02T06:23:11Z

pygmt/tests/test_grdcut.py

+    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


Could you add a more precise assertion(s)? Ideally, you want to check the exact pixel count, and maybe also check that there are NaN values in the four corners or something.

I understand, for the polygon each vertex will need to not be NaN

weiji14 · 2025-10-02T06:29:58Z

pygmt/src/grdcut.py

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


I'm not sure if that specific CRS is necessary, of if a CRS is needed actually. Other input types (ASCII file, shapely.geometry.Polygon) don't even have a CRS.

weiji14 · 2025-10-02T06:30:41Z

pygmt/src/grdcut.py

+        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


Could you add a unit test to confirm that MultiPolygon does work?

kshitijaucharmal · 2025-10-02T07:29:28Z

Thanks for adding the unit tests! Here are some more comments, and apologies for the failing CI tests, we have some issues with DVC running on forks that we'll probably try to sort out next week.

Yeah I was wondering why all the test keep failing
Anyways, will get back to you after resolving the comments

- Reusable polygon fixture for grdcut tests - Improve assertions to check pixel counts and presence of NaNs outside polygons - Add a test for MultiPolygon support

kshitijaucharmal · 2025-10-02T08:01:29Z

I do not quite understand the precise assertions part, is checking for some pixels inside polygon to not be NaN enough? I hope the new cases are correct or at least in the right direction. Thanks !

kshitijaucharmal · 2025-10-16T06:34:31Z

Hey @weiji14 are the CI/CD issues not fixed yet? Please let me know if I can help in any way, thanks !

pygmt/src/grdcut.py

seisman · 2025-10-16T06:46:17Z

Hey @weiji14 are the CI/CD issues not fixed yet? Please let me know if I can help in any way, thanks !

We still have issues with CI, but you can just go ahead to improve this PR. We will run the tests locally as a temporary workaround.

kshitijaucharmal · 2025-10-18T07:08:52Z

@seisman @weiji14 Will fix the issues ASAP

seisman · 2025-11-07T09:19:13Z

Kindly ping @kshitijaucharmal to see if you have time to continue this PR. The CI issue is likely fixed.

kshitijaucharmal · 2025-11-07T09:25:50Z

Hey @seisman sorry for the delay, will be fixing the issues in the PR today as soon as I can

kshitijaucharmal · 2025-11-07T14:02:45Z

Have made the changes requested, please check

seisman

I've made some suggestions, mainly to migrate the polgyon parameter to the new alias system

pygmt/src/grdcut.py

kshitijaucharmal · 2025-11-07T16:47:04Z

Have made the changes, the new alias system is much better honestly :)

seisman · 2025-11-08T05:04:54Z

pygmt/src/grdcut.py

+                if isinstance(polygon, PathLike) or hasattr(
+                    polygon, "__geo_interface__"
+                ):
+                    aliasdict["F"] = str(vinpoly) + aliasdict.get("F", "")
+                else:
+                    polygon_type = type(polygon).__name__
+                    err_msg = "Invalid polygon type"
+                    raise GMTValueError(
+                        err_msg,
+                        polygon_type,
+                        "Must be a PathLike, GeoDataFrame, or Shapely geometry.",
+                    )


This check can be done by the data_kind function, e.g.,

if polyon is not None and data_kind(polygon) not in ("file", "geojson"): raise GMTTypeError(...)

Ideally we should raise the error as early as possible, so please put it at line 140.

seisman · 2025-11-08T05:06:49Z

pygmt/src/grdcut.py

+            if polygon is not None:
+                if isinstance(polygon, PathLike) or hasattr(
+                    polygon, "__geo_interface__"
+                ):
+                    aliasdict["F"] = str(vinpoly) + aliasdict.get("F", "")
+                else:
+                    polygon_type = type(polygon).__name__
+                    err_msg = "Invalid polygon type"
+                    raise GMTValueError(
+                        err_msg,
+                        polygon_type,
+                        "Must be a PathLike, GeoDataFrame, or Shapely geometry.",
+                    )


Please see my comment above. The check can be done earlier.

Suggested change

if polygon is not None:

if isinstance(polygon, PathLike) or hasattr(

polygon, "__geo_interface__"

):

aliasdict["F"] = str(vinpoly) + aliasdict.get("F", "")

else:

polygon_type = type(polygon).__name__

err_msg = "Invalid polygon type"

raise GMTValueError(

err_msg,

polygon_type,

"Must be a PathLike, GeoDataFrame, or Shapely geometry.",

)

if polygon:

aliasdict["F"] = str(vinpoly) + aliasdict.get("F", "")

seisman · 2025-11-08T07:06:31Z

pygmt/src/grdcut.py

    build_arg_list,
    data_kind,
    fmt_docstring,
+    tempfile_from_geojson,


Suggested change

tempfile_from_geojson,

seisman · 2025-11-08T07:07:58Z

pygmt/src/grdcut.py

            outkind = kind
        case _:
            raise GMTTypeError(type(grid))



Suggested change

if polygon and data_kind(polygon) not in ("file", "geojson"):

raise GMTTypeError(

type(polygon),

reason="Must be a PathLike, GeoDataFrame, or Shapely geometry.",

)

kshitijaucharmal added 2 commits October 1, 2025 10:58

Add polygon-based subregion support to grdcut with gallery example

f256ab4

Add crop and invert boolean support for polygon input

bac59da

yvonnefroehlich assigned kshitijaucharmal Oct 1, 2025

yvonnefroehlich added hacktoberfest Hacking away in the month of October enhancement Improving an existing feature labels Oct 1, 2025

yvonnefroehlich added this to the 0.18.0 milestone Oct 1, 2025

weiji14 reviewed Oct 1, 2025

View reviewed changes

kshitijaucharmal added 2 commits October 2, 2025 11:02

grdcut: PathLike for polygon input and improve example

f6d4179

- 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

weiji14 marked this pull request as draft October 2, 2025 06:18

weiji14 reviewed Oct 2, 2025

View reviewed changes

Add polygon fixture and more precise grdcut unit tests

4e45607

- Reusable polygon fixture for grdcut tests - Improve assertions to check pixel counts and presence of NaNs outside polygons - Add a test for MultiPolygon support

kshitijaucharmal marked this pull request as ready for review October 3, 2025 10:00

kshitijaucharmal requested a review from weiji14 October 3, 2025 10:01

kshitijaucharmal marked this pull request as draft October 3, 2025 10:01

kshitijaucharmal marked this pull request as ready for review October 3, 2025 13:44

Merge branch 'main' into grdcut-polygon

e564199

seisman mentioned this pull request Oct 9, 2025

Running "dvc pull" now requires authentication #4147

Open

weiji14 mentioned this pull request Oct 10, 2025

Set CI Tests workflow to run under pr-tests github environment #4148

Merged

6 tasks

Merge branch 'main' into grdcut-polygon

5078cb4

weiji14 marked this pull request as draft October 13, 2025 20:33

weiji14 had a problem deploying to pr-tests October 13, 2025 20:34 — with GitHub Actions Failure

seisman reviewed Oct 16, 2025

View reviewed changes

pygmt/src/grdcut.py Show resolved Hide resolved

pygmt/src/grdcut.py Outdated Show resolved Hide resolved

pygmt/src/grdcut.py Outdated Show resolved Hide resolved

kshitijaucharmal and others added 2 commits November 7, 2025 19:24

Merge branch 'GenericMappingTools:main' into grdcut-polygon

f6d9fa3

removed ExitStack

3bb5ab9

seisman reviewed Nov 7, 2025

View reviewed changes

migrate polygon param to new alias system

c9bee38

seisman reviewed Nov 8, 2025

View reviewed changes

seisman reviewed Nov 9, 2025

View reviewed changes

+    if polygon and data_kind(polygon) not in ("file", "geojson"):
+        raise GMTTypeError(
+            type(polygon),
+            reason="Must be a PathLike, GeoDataFrame, or Shapely geometry.",
+        )

Support polygon-based subsetting in grdcut -F with GeoDataFrame #4135

Are you sure you want to change the base?

Support polygon-based subsetting in grdcut -F with GeoDataFrame #4135

Uh oh!

Conversation

kshitijaucharmal commented Oct 1, 2025 • edited by weiji14 Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

kshitijaucharmal commented Oct 1, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

weiji14 left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

kshitijaucharmal commented Oct 2, 2025

Uh oh!

weiji14 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

kshitijaucharmal commented Oct 2, 2025

Uh oh!

kshitijaucharmal commented Oct 2, 2025

Uh oh!

kshitijaucharmal commented Oct 16, 2025

Uh oh!

Uh oh!

Uh oh!

Uh oh!

seisman commented Oct 16, 2025

Uh oh!

kshitijaucharmal commented Oct 18, 2025

Uh oh!

seisman commented Nov 7, 2025

Uh oh!

kshitijaucharmal commented Nov 7, 2025

Uh oh!

kshitijaucharmal commented Nov 7, 2025

Uh oh!

seisman left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

kshitijaucharmal commented Nov 7, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

seisman Nov 8, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

kshitijaucharmal commented Oct 1, 2025 •

edited by weiji14

Loading

kshitijaucharmal commented Oct 1, 2025 •

edited

Loading

seisman left a comment •

edited

Loading

seisman Nov 8, 2025 •

edited

Loading