终于解决了我的问题。
这是我正在寻找的情节和使用的脚本。
这里有一个有用的链接，用来解决我的问题
Python cartopy map, clip area outside country (polygon)
enter image description here
import numpy as np
from scipy import stats
from netCDF4 import Dataset
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy
import cartopy.io.shapereader as shapereader
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.font_manager import FontProperties
import geopandas
from matplotlib.ticker import (AutoMinorLocator, MultipleLocator)
import matplotlib.ticker as mticker
from shapely.geometry import Polygon
from cartopy.io import shapereader
import cartopy.io.img_tiles as cimgt
import geopandas
def rect_from_bound(xmin, xmax, ymin, ymax):
"""Returns list of (x,y)'s for a rectangle"""
xs = [xmax, xmin, xmin, xmax, xmax]
ys = [ymax, ymax, ymin, ymin, ymax]
return [(x, y) for x, y in zip(xs, ys)]
ref_dat = Dataset("R95ptot_Kmer_chirps-v2.0.1981_2019.days_p05.nc", mode='r')
lat = ref_dat.variables["latitude"][:]
lon = ref_dat.variables["longitude"][:]
time = ref_dat.variables["time"]
pre0 = ref_dat.variables["precip"][:, :, :]
# request data for use by geopandas
resolution = '10m'
category = 'cultural'
name = 'admin_0_countries'
shpfilename = shapereader.natural_earth(resolution, category, name)
df = geopandas.read_file(shpfilename)
# get geometry of a country
poly = [df.loc[df['ADMIN'] == 'Cameroon']['geometry'].values[0]]
#stamen_terrain = cimgt.Stamen('terrain-background')
# projections that involved
st_proj = ccrs.PlateCarree()#stamen_terrain.crs #projection used by Stamen images
ll_proj = ccrs.PlateCarree() #CRS for raw long/lat
# create fig and axes using intended projection
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(1, 1, 1, projection=st_proj)
ax.add_geometries(poly, crs=ll_proj, facecolor='none', edgecolor='black')
pad1 = .1 #padding, degrees unit
exts = [poly[0].bounds[0] - pad1, poly[0].bounds[2] + pad1, poly[0].bounds[1] - pad1, poly[0].bounds[3] + pad1];
#exts = [min(lon), max(lon), min(lat), max(lat)]
ax.set_extent(exts, crs=ccrs.PlateCarree())
ax.set_extent(exts, crs=ll_proj)
# make a mask polygon by polygon's difference operation
# base polygon is a rectangle, another polygon is simplified switzerland
msk = Polygon(rect_from_bound(*exts)).difference( poly[0].simplify(0.01) )
msk_stm = st_proj.project_geometry (msk, ll_proj) # project geometry to the projection used by stamen
# get and plot Stamen images
#ax.add_image(stamen_terrain, 8) # this requests image, and plot
slope = np.zeros(pre0[0, :, :].shape)
intercept = np.zeros(pre0[0, :, :].shape)
r_value = np.zeros(pre0[0, :, :].shape)
p_value = np.zeros(pre0[0, :, :].shape)
std_err = np.zeros(pre0[0, :, :].shape)
for ilat in range(len(lat)):
for jlon in range(len(lon)):
slope[ilat, jlon], intercept[ilat, jlon], r_value[ilat, jlon], p_value[ilat, jlon], std_err[ilat, jlon] = stats.linregress(time, pre0[:, ilat, jlon])
[lons2d, lats2d] = np.meshgrid(lon, lat)
mymap = ax.contourf(lons2d, lats2d, slope,
transform=ccrs.PlateCarree(),
cmap=plt.cm.Spectral_r, extend='both')
#ax.set_extent([min(lon), max(lon), min(lat), max(lat)], crs=ccrs.PlateCarree())
#cbarax = fig.add_axes([0.2, 0.95, 0.6, 0.02])
cbar = plt.colorbar(mymap, orientation='vertical')
cbar.set_label('Precipitation [mm]', rotation=90, fontsize=12, labelpad=1)
cbar.ax.tick_params(labelsize=12, length=0)
pfieldm = np.ma.masked_greater(p_value, 0.05)
ax.contourf(lons2d, lats2d, pfieldm, transform = ccrs.PlateCarree(), hatches=["..."], alpha=0.0)
ax.add_geometries( msk_stm, st_proj, zorder=12, facecolor='white', edgecolor='none', alpha=1.)
ax.outline_patch.set_edgecolor('white')
#ax.outline_patch.set_visible(False)
plt.savefig('figure.pdf', format='pdf', dpi=500)
plt.show()