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

[weiji14]

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

  pygmt/pygmt/src/grdcut.py

  Lines 17 to 18 in bafb8ab

  	@use_alias(
  	G="outgrid",

• Implement the GeoDataFrame -> ASCII file conversion that can be read by gmt grdcut -F, may need to use tempfile_from_geojson or some alternative way.
• Add a gallery example for this to https://github.com/GenericMappingTools/pygmt/tree/main/examples/gallery/images