Python（Cartopy）在特定国家/地区内绘制阴影图形[已解决]

Question

我正试图绘制喀麦隆上空的降水图，只给喀麦隆边界的内侧加上阴影。
所有其他国家都将被屏蔽。
下面是我正在使用的脚本的情节和开头：

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
    

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"][:, :, :]    
        
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) 

fig, ax = plt.subplots(1, 1, sharex=True, sharey=True, figsize=(8.5, 6.98),
subplot_kw={'projection':ccrs.PlateCarree()})

mymap = ax.contourf(lons2d, lats2d, slope, 
                              transform=ccrs.PlateCarree(), 
                              cmap=plt.cm.Spectral_r, extend='both')

Answer 1

最后我解决了我的问题。 这里是我要寻找的情节和使用的脚本

这里有一个有用的链接，用来解决我的问题

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