1998年夏季降水距平

1998年夏季降水距平

效果图

导入包

import netCDF4 as nc
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import cartopy.feature as cfeature
from matplotlib.colors import ListedColormap
import cartopy.io.shapereader as shpreader
import matplotlib as mpl
import maskout
import matplotlib.ticker as mticker

设置中文字体

#  设置画图字体的大小
plt.rcParams.update({'font.size': 20})
#  解决中文乱码问题
plt.rcParams['font.sans-serif'] = ['SimHei']
#  解决负号乱码问题
plt.rcParams['axes.unicode_minus'] = False

读取数据

# 读取数据
file = r'pr_wtr.eatm.mon.mean.nc'
data = nc.Dataset(file)
lats = data.variables['lat'][:]
lons = data.variables['lon'][:]
p = data.variables['pr_wtr'][:]

数据处理

# 计算多年平均夏季 64年
p_ave_summer = np.array([])
for i in range(64):
    for j in [5, 6, 7]:
        p_ave_summer = np.append(p_ave_summer, p[i * 12 + j])
p_ave_summer = p_ave_summer.reshape(64 * 3, 73, 144)
p_ave_summer = np.mean(p_ave_summer, axis=0)

# 计算98年夏季平均
p_98_ave = np.mean(p[605:608, :, :], axis=0)

# 距平
p_x = p_98_ave - p_ave_summer

应该有更加简单的方式筛选出每年的678月，笔者还在学习，如果有大佬知道还请不吝赐教。

导入地图

# 导入地图
china = shpreader.Reader('provinces.shp').geometries()
rivers = shpreader.Reader('rivers.shp').geometries()
proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(6, 4.5))
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.add_geometries(china, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
ax.add_geometries(rivers, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=0.75, zorder=5)

画图

cmap = mpl.cm.RdBu
newcolors = cmap(np.linspace(0, 1, 256))
newcmp = ListedColormap(newcolors[50:206])
lons_grid, lats_grid = np.meshgrid(lons, lats)
cf = ax.contourf(lons_grid, lats_grid, p_x,
                 levels=[-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7],
                 transform=proj, extend='both',
                 colors=['#de735c', '#f2a17f', '#f9c4a9', '#fce2d2', '#f8f4f2', '#e3edf3', '#d1e5f0', '#b6d7e8',
                         '#9bc9e0', '#7bb6d6', '#68abd0', '#4695c4', '#3681ba', '#276eb0', '#1a5899', '#134c87',
                         '#0e4179', '#073467'])
cb = fig.colorbar(cf, shrink=0.6, orientation='horizontal', pad=0.08)
cb.ax.set_xlabel('大气可降水量距平（mm）', fontweight='bold')

白化

clip = maskout.shp2clip(cf, ax, 'china0.shp')  # 白化

坐标轴设置

for i in [i * 5 for i in range(11)]:
    plt.axhline(y=i, c='grey', ls='--', lw=0.5)
for i in [i * 10 + 70 for i in range(7)]:
    plt.axvline(x=i, c='grey', ls='--', lw=0.5)
plt.xticks(ticks=[i * 10 + 70 for i in range(7)], labels=[i * 10 + 70 for i in range(7)])
plt.yticks(ticks=[i * 5 for i in range(11)], labels=[i * 5 for i in range(11)])
plt.gca().yaxis.set_major_formatter(mticker.FormatStrFormatter('%.0f°N'))
plt.gca().xaxis.set_major_formatter(mticker.FormatStrFormatter('%.0f°E'))
plt.ylim([0, 55])
plt.xlim([70, 135])
plt.title('98年夏季（JJA）大气可降水量距平图')
plt.show()

总代码

import netCDF4 as nc
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import cartopy.feature as cfeature
from matplotlib.colors import ListedColormap
import cartopy.io.shapereader as shpreader
import matplotlib as mpl
import maskout
import matplotlib.ticker as mticker
import matplotlib.colors as mcolors

#  设置画图字体的大小
plt.rcParams.update({'font.size': 20})
#  解决中文乱码问题
plt.rcParams['font.sans-serif'] = ['SimHei']
#  解决负号乱码问题
plt.rcParams['axes.unicode_minus'] = False

# 读取数据
file = r'pr_wtr.eatm.mon.mean.nc'
data = nc.Dataset(file)
lats = data.variables['lat'][:]
lons = data.variables['lon'][:]
p = data.variables['pr_wtr'][:]

# 计算多年平均夏季 64年
p_ave_summer = np.array([])
for i in range(64):
    for j in [5, 6, 7]:
        p_ave_summer = np.append(p_ave_summer, p[i * 12 + j])
p_ave_summer = p_ave_summer.reshape(64 * 3, 73, 144)
p_ave_summer = np.mean(p_ave_summer, axis=0)

# 计算98年夏季平均
p_98_ave = np.mean(p[605:608, :, :], axis=0)

# 距平
p_x = p_98_ave - p_ave_summer

# 导入地图
china = shpreader.Reader('provinces.shp').geometries()
rivers = shpreader.Reader('rivers.shp').geometries()
proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(6, 4.5))
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.add_geometries(china, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
ax.add_geometries(rivers, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=0.75, zorder=5)
cmap = mpl.cm.RdBu
newcolors = cmap(np.linspace(0, 1, 256))
newcmp = ListedColormap(newcolors[50:206])
lons_grid, lats_grid = np.meshgrid(lons, lats)
cf = ax.contourf(lons_grid, lats_grid, p_x,
                 levels=[-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7],
                 transform=proj, extend='both',
                 colors=['#de735c', '#f2a17f', '#f9c4a9', '#fce2d2', '#f8f4f2', '#e3edf3', '#d1e5f0', '#b6d7e8',
                         '#9bc9e0', '#7bb6d6', '#68abd0', '#4695c4', '#3681ba', '#276eb0', '#1a5899', '#134c87',
                         '#0e4179', '#073467'])
cb = fig.colorbar(cf, shrink=0.6, orientation='horizontal', pad=0.08)
cb.ax.set_xlabel('大气可降水量距平（mm）', fontweight='bold')
clip = maskout.shp2clip(cf, ax, 'china0.shp')  # 白化
for i in [i * 5 for i in range(11)]:
    plt.axhline(y=i, c='grey', ls='--', lw=0.5)
for i in [i * 10 + 70 for i in range(7)]:
    plt.axvline(x=i, c='grey', ls='--', lw=0.5)
# 坐标轴设置
plt.xticks(ticks=[i * 10 + 70 for i in range(7)], labels=[i * 10 + 70 for i in range(7)])
plt.yticks(ticks=[i * 5 for i in range(11)], labels=[i * 5 for i in range(11)])
plt.gca().yaxis.set_major_formatter(mticker.FormatStrFormatter('%.0f°N'))
plt.gca().xaxis.set_major_formatter(mticker.FormatStrFormatter('%.0f°E'))
plt.ylim([0, 55])
plt.xlim([70, 135])
plt.title('98年夏季（JJA）大气可降水量距平图')
plt.show()