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
) is a common operation, e.g. to extract a Region of Interest.
Currently, the pygmt.grdcut
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
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 😄
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) atLines 17 to 18 in bafb8ab
- Implement the GeoDataFrame -> ASCII file conversion that can be read by
gmt grdcut -F
, may need to usetempfile_from_geojson
or some alternative way. - Add a gallery example for this to https://github.com/GenericMappingTools/pygmt/tree/main/examples/gallery/images