GenericMappingTools / GenericMappingTools/pygmt

Subset grid using a geopandas.GeoDataFrame polygon with grdcut -F

Open
#1,529 2 comments 0 reactions 1 assignee Claimed by @kshitijaucharmal View on GitHub
feature request good first issue help wanted
Dominant language
Python
Stars
874
Forks
255
Avg merge
1d 21h
Merged PRs (30d)
40

Description

**Description of the desired feature**

Spatial subsetting/cropping/masking/clipping of a raster grid using a polygon (stored in an ESRI shapefile or [`geopandas.GeoDataFrame`](https://geopandas.org/docs/reference/api/geopandas.GeoDataFrame.html#geopandas.GeoDataFrame)) is a common operation, e.g. to extract a Region of Interest.

Currently, the [`pygmt.grdcut`](https://www.pygmt.org/v0.4.1/api/generated/pygmt.grdcut.html) function in PyGMT v0.4.1 (added in #492) can do bounding box (via 'region' `-R`), circular (via 'circ_subregion' `-S`) and z-value (via 'z_subregion' `-Z`) based subsetting.

It would be nice to add in the [`grdcut -F`](https://docs.generic-mapping-tools.org/latest/grdcut.html#f) option too which does polygon-based crops using an ASCII file. Furthermore, it would be even nicer to perform the crop using an in-memory `geopandas.GeoDataFrame` or `shapely.geometry` polygon object!

Technically, this is reinventing the wheel, because there are a few Python libraries that already do this:
- https://corteva.github.io/rioxarray/stable/examples/clip_geom.html
- https://rasterio.readthedocs.io/en/latest/topics/masking-by-shapefile.html#masking-a-raster-using-a-shapefile
- https://geosynopsis.github.io/xgeo/docs/html/examples.html#subset-data (not very maintained)

But who says the wheel can't be reinvented :smile:

**Are you willing to help implement and maintain this feature?** A bit low on my priority list, but anyone interested is welcome to try.

Steps to implement:
- [ ] Alias the `-F` parameter in grdcut to "polygon" (or another name) at https://github.com/GenericMappingTools/pygmt/blob/bafb8ab906aa87eb1dab3ce286b9baa91e145099/pygmt/src/grdcut.py#L17-L18
- [ ] Implement the GeoDataFrame -> ASCII file conversion that can be read by `gmt grdcut -F`, may need to use [`tempfile_from_geojson`](https://github.com/GenericMappingTools/pygmt/blob/a0286345933c1eaf2e1fb215809e33f2a792be0b/pygmt/helpers/tempfile.py#L111) or some alternative way.
- [ ] Add a gallery example for this to https://github.com/GenericMappingTools/pygmt/tree/main/examples/gallery/images

Contributor guide

Open the contributing guide

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.