Hello,
Is “inside one of the polygons in the polygonfile” option available for pygmt.select … I was expecting a “polygon=polygon” option, but I’m not seeing one that would match GMT’s “-F” option.
Thanks!
The -F option is not yet wrapped for select. In the meantime, you could use F=polygonfile where polygonfile is the name of the polygon file. The option to use single characters for GMT options in PyGMT is undocumented because it will be disallowed in a future release, so no gaurantees of stability here. If you do chose to use it in your code, you will be warned when a parameter becomes available for the -F option.
2 Likes
Here’s a work-around to extract a grid from within a shapefile polygon (‘shape_file’): You’ll need to have the initial grid as an Xarray dataarray (da_all), and “import rioxarray as rio” and “from shapely.geometry import mapping” (as well as Pandas, Geopandas, Xarray and Matplotlib)
# [1] Read polygon shapefile to a Pandas geodataframe
gdf = gpd.read_file(filename='shape_file')
# [2] Ensure the Xarray dataarray has the same CRS as the geodataframe
da_all = da_all.rio.write_crs(gdf.crs)
# [3] Clip grid to polygon
da_insidePoly = da_all.rio.clip(gdf.geometry.apply(mapping), gdf.crs, drop=True, invert=False)
# [4] QC results
fig = plt.figure(figsize=(16,5))
gs = fig.add_gridspec(1, 2)
### ------------------Map 1
ax = fig.add_subplot(gs[0, 0])
da_all.plot.pcolormesh('x', 'y', cmap='nipy_spectral', ax=ax)
ax.set_title('All area')
ax.set_aspect('equal')
### ------------------Map 2
ax = fig.add_subplot(gs[0, 1])
da_insidePoly.plot.pcolormesh('x', 'y', cmap='nipy_spectral', ax=ax)
ax.set_title('Inside polygon')
ax.set_aspect('equal')
plt.show()