python - 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')

标签： pythoncartopy

Answer 1

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

这是一个有用的链接，用于解决我的问题

Python cartopy 地图，国家以外的剪辑区域（多边形）

在此处输入图像描述

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