基于Python的NCEP再分析数据的中国区域白化（含南海小地图）

'''
Copyright  2020.
@author:Developed by Haolisheng, 2020.12.18,tianjin
'''
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import xarray as xr
import maskout31
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from matplotlib import rcParams
config = {"font.family":'Times New Roman', "font.size": 12, "mathtext.fontset":'stix'}
rcParams.update(config)
region=[70, 140, 15, 55]
df= xr.open_dataset(r'F:/Rpython/lp27/data/vwnd.mon.ltm.nc')
lat0 = df['lat']  
lon0 = df['lon'] 
z = df['vwnd'][6,2,::]     #取8月、500hPa等，时次,层,经度,纬度
#scale='10m'
fig=plt.figure(dpi = 600, figsize = (12, 6))
proj=ccrs.PlateCarree()
ax = plt.axes(projection = proj)
ax.set_extent(region, crs = proj)
#ax.coastlines(scale, linewidth = 0.5)
ax.set_xticks(np.arange(region[0], region[1] + 1, 10), crs = proj)
ax.set_yticks(np.arange(region[-2], region[-1] + 1, 10), crs = proj)
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=False))
ax.yaxis.set_major_formatter(LatitudeFormatter())
# Add ocean, land, rivers and lakes
#ax.add_feature(cfeature.COASTLINE,lw=0.3) 
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))
ax.add_geometries(Reader(r'F:/Rpython/lp27/data/river1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='b',linewidth=0.2)
ax.add_geometries(Reader(r'F:/Rpython/lp27/data/china1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=0.2)
ax.add_geometries(Reader(r'F:/Rpython/lp27/data/china2.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=0.1)
ax.add_geometries(Reader(r'F:/Rpython/lp27/data/ne_10m_land.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=0.2)
ax.add_geometries(Reader(r'F:/Rpython/lp27/data/ne_50m_lakes.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=0.2)
lon_range1=lon0[(lon0>60)&(lon0<150)]
lat_range1=lat0[(lat0>0)&(lat0<60)]
z1=z.sel(lon=lon_range1,lat=lat_range1)
lon1,lat1=np.meshgrid(lon_range1,lat_range1)
z1=z1.data
c11=ax.contourf(lon1,lat1,z1,np.arange(-6,6.01,0.5),extend='both',transform=ccrs.PlateCarree(),cmap='gist_rainbow')
clip=maskout31.shp2clip(c11,ax,'F:/Rpython/lp27/data/china0')
cbar=plt.colorbar(c11,shrink=0.99,aspect=20,fraction=.03,pad=0.01)   #aspect控制bar宽度，fraction控制大小比例,pad控制与图的距离
cbar.set_ticks(np.arange(-5,5.001,1))            #设置colorbar范围和刻度标记间隔
cbar.ax.tick_params(labelsize=12, direction='in', right=False)
ax.set_title('NCEP再分析1981-2010年7月平均850hPa经向风速. 单位：m.s$^{-1}$',{'family':'simhei','size':14,'color':'k'}) 
ax0 = plt.gca()   #获取边框
ax0.outline_patch.set_linewidth(0.5)    #修改边框粗细         
#添加南海，实际上就是新建一个子图覆盖在之前子图的右下角
ax2 = fig.add_axes([0.7365, 0.11, 0.2, 0.2],projection=proj)
ax2.set_extent([105.8, 122,0,25])
ax2.add_feature(cfeature.OCEAN.with_scale('50m'))
ax2.add_geometries(Reader(r'F:/Rpython/lp27/data/china1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=0.2)
ax2.add_geometries(Reader(r'F:/Rpython/lp27/data/ne_10m_land.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=0.2)
lon_range2=lon0[(lon0>=105)&(lon0<=125)]
lat_range2=lat0[(lat0>=0)&(lat0<=25)]
z2=z.sel(lon=lon_range2,lat=lat_range2)
lon2,lat2=np.meshgrid(lon_range2,lat_range2)
z2=z2.data
c21=ax2.contourf(lon2,lat2,z2,np.arange(-6,6.01,0.5),extend='both',transform=ccrs.PlateCarree(),cmap='gist_rainbow')
clip=maskout31.shp2clip(c21,ax2,'F:/Rpython/lp27/data/china0')
ax0 = plt.gca()   #获取边框
ax0.outline_patch.set_linewidth(0.5)    #修改边框粗细
#保存图片
plt.savefig(r"F:/Rpython/lp27/plot70.4.png",dpi=600)
plt.show()

ax.set_xticks(np.arange(region[0], region[1]+10, 10))
ax.set_yticks(np.arange(region[2], region[3]+10, 10))
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.tick_params(axis='both',labelsize=12,direction='in',right=False,top=False)
ax.grid(linewidth=0.4, color='r', alpha=0.45, linestyle='-')