GenericMappingTools / GenericMappingTools/pygmt
Subset grid using a geopandas.GeoDataFrame polygon with grdcut -F
- 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
Assessment
This issue has not been assessed yet.