Support polygon-based subsetting in grdcut -F with GeoDataFrame #4135
Conversation
yvonnefroehlich
added
hacktoberfest
|
@yvonnefroehlich I have added both of the Also, the label is |
There was a problem hiding this comment.
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!
- 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
|
@weiji14 I have resolved all comments and added the unit tests as well, let me know if there are any more corrections! :) |
pygmt/tests/test_grdcut.py
Outdated
| # 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), | ||
| ] | ||
| ) |
There was a problem hiding this comment.
This should be turned into a pytest fixture, see examples a few lines below.
There was a problem hiding this comment.
Will convert this into a fixture
| 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 |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
I understand, for the polygon each vertex will need to not be NaN
pygmt/src/grdcut.py
Outdated
| 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) |
There was a problem hiding this comment.
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.
| 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 |
There was a problem hiding this comment.
Could you add a unit test to confirm that MultiPolygon does work?
Yeah I was wondering why all the test keep failing |
- Reusable polygon fixture for grdcut tests - Improve assertions to check pixel counts and presence of NaNs outside polygons - Add a test for MultiPolygon support
|
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 ! |
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:
Additionally, two new boolean parameters,
cropandinvert, 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.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
make formatandmake checkto make sure the code follows the style guide.doc/api/index.rst.