短期气候预测实习4

xiaoxiao啊

基于jupyter notebook
数据获取方式如下
链接：https://pan.baidu.com/s/17TSNCN-XXV5zrMXhIG9iTA
提取码：0un5


# 导入所需库
import numpy as np
import pandas as pd
import xarray as xr
from scipy.stats.mstats import ttest_ind
from scipy.interpolate import Rbf
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cnmaps import get_map, draw_map, clip_contours_by_map, clip_pcolormesh_by_map
import cmaps

import warnings
warnings.filterwarnings('ignore')


# 更改字体
plt.rcParams['font.sans-serif'] = ['Times New Roman']
plt.rcParams['axes.unicode_minus'] = False


# 读取站点数据
df = pd.read_csv('station.csv', encoding='gbk', names=['station', 'num', 'lat', 'lon'], header=0)
df.head()  # 查看前五行


# 1951-2010(60 years)
f6 = np.loadtxt('r1606.txt').reshape(71, 160)[:-11]
f7 = np.loadtxt('r1607.txt').reshape(71, 160)[:-11]
f8 = np.loadtxt('r1608.txt').reshape(70, 160)[:-10]

# 1951-2010夏季逐年降水站点数据
pr = np.stack([f6, f7, f8], axis=0).mean(0)
pr = xr.DataArray(pr, coords=[range(1951, 2011), df.station], dims=['year', 'station'])


# 站点经纬度
lon = df.lon
lat = df.lat

# 格点经纬度
grid_lon = np.linspace(73, 135, 120)
grid_lat = np.linspace(15, 55, 80)
meshgrid_lon, meshgrid_lat = np.meshgrid(grid_lon, grid_lat)


f2 = np.loadtxt('ddi.txt', dtype='int32')
_, pr_type = np.where(f2 == 1)
year = np.arange(1951, 2011)
yr1 = year[pr_type == 0]
yr2 = year[pr_type == 1]
yr3 = year[pr_type == 2]


pr_mean = pr.loc['1971':'2000'].mean('year')
RP = (pr - pr_mean) / pr_mean * 100
RP1 = RP.loc[RP.year.isin(yr1)]
RP2 = RP.loc[RP.year.isin(yr2)]
RP3 = RP.loc[RP.year.isin(yr3)]
RP_ano1 = RP1.mean('year') - RP.mean('year')
RP_ano2 = RP2.mean('year') - RP.mean('year')
RP_ano3 = RP3.mean('year') - RP.mean('year')


_, RP_p1 = ttest_ind(RP1, RP, axis=0, equal_var=False)
_, RP_p2 = ttest_ind(RP2, RP, axis=0, equal_var=False)
_, RP_p3 = ttest_ind(RP3, RP, axis=0, equal_var=False)


def my_interpolate(*, var):
    func = Rbf(df.lon, df.lat, var, function='cubic')
    grid_var = func(meshgrid_lon, meshgrid_lat)
    return grid_var


grid_RP_ano1 = my_interpolate(var=RP_ano1)
grid_RP_ano2 = my_interpolate(var=RP_ano2)
grid_RP_ano3 = my_interpolate(var=RP_ano3)
grid_RP_p1 = my_interpolate(var=RP_p1)
grid_RP_p2 = my_interpolate(var=RP_p2)
grid_RP_p3 = my_interpolate(var=RP_p3)


def pr_plot(*, ano, p, num):
    levels = list(np.arange(-25, 30, 5))
    fig = plt.figure(figsize=(6, 5.5), facecolor='w')
    ax1 = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=ccrs.PlateCarree(105))
    ax1.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5, edgecolor='dimgrey')
    ax1.add_feature(cfeature.LAKES.with_scale('110m'), linewidth=0.5, edgecolor='dimgrey', facecolor='none')
    ax1.set_extent([70, 140, 15, 55], ccrs.PlateCarree())
    ax1.set_title(f'RP Anomaly of type {num}', loc='left')
    ax1.set_xticks(np.arange(70, 145, 10), crs=ccrs.PlateCarree())
    ax1.set_yticks(np.arange(15, 60, 10), crs=ccrs.PlateCarree())
    ax1.xaxis.set_major_formatter(LongitudeFormatter())
    ax1.yaxis.set_major_formatter(LatitudeFormatter())
    ax1.tick_params(which='major', width=0.5, length=5)
    ax1.spines['geo'].set_linewidth(0.5)
    ax1.set_aspect(1.2)
    cf = ax1.contourf(grid_lon, grid_lat, ano, cmap='BrBG', extend='both', transform=ccrs.PlateCarree(), levels=levels)
    dot_area = np.where(p < 0.05)
    dot = ax1.scatter(meshgrid_lon[dot_area], meshgrid_lat[dot_area], color='k', s=1.2, linewidths=0, transform=ccrs.PlateCarree())
    draw_map(get_map('中国'), color='dimgrey', linewidth=0.5)
    draw_map(get_map('南海'), color='k', linewidth=0.6)
    clip_contours_by_map(cf, get_map('中国'))
    clip_pcolormesh_by_map(dot, get_map('中国'))
    ax1.add_geometries(Reader('河流.shp').geometries(), crs=ccrs.PlateCarree(), edgecolor='dimgrey', facecolor='none', lw=0.5)

    ax2 = fig.add_axes([0.7679, 0.201, 0.15, 0.15], projection=ccrs.PlateCarree())
    ax2.set_extent([105, 124, 0, 23], ccrs.PlateCarree())
    ax2.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5, edgecolor='dimgrey')
    ax2.spines['geo'].set_linewidth(0.5)
    draw_map(get_map('南海'), color='k', linewidth=0.5)
    draw_map(get_map('中国'), color='dimgrey', linewidth=0.6)

    ax3 = fig.add_axes([0.15, 0.075, 0.7, 0.02])
    cb = fig.colorbar(cf, cax=ax3, orientation='horizontal', shrink=0.9, aspect=30, ticks=levels, pad=0.12)
    cb.outline.set_linewidth(0.1)
    cb.ax.tick_params(length=0)

pr_plot(ano=grid_RP_ano1, p=grid_RP_p1, num=1)
pr_plot(ano=grid_RP_ano2, p=grid_RP_p2, num=2)
pr_plot(ano=grid_RP_ano3, p=grid_RP_p3, num=3)


f3 = xr.open_dataset('hgt.mon.mean.nc')
z = f3.hgt.loc[f3.time.dt.month.isin([1]), 500].loc['1951':'2010']
z1 = z.loc[z.time.dt.year.isin(yr1)]
z2 = z.loc[z.time.dt.year.isin(yr2)]
z3 = z.loc[z.time.dt.year.isin(yr3)]


z_ano1 = z1.mean('time') - z.mean('time')
z_ano2 = z2.mean('time') - z.mean('time')
z_ano3 = z3.mean('time') - z.mean('time')
_, z_p1 = ttest_ind(z1, z, equal_var=False)
_, z_p2 = ttest_ind(z2, z, equal_var=False)
_, z_p3 = ttest_ind(z3, z, equal_var=False)

mesh_lon, mesh_lat = np.meshgrid(f3.lon, f3.lat)


def z500_plot(*, ano, p, num):
    levels = list(np.arange(-50, 55, 10))
    fig = plt.figure(figsize=(10, 5), facecolor='w')
    ax = fig.add_subplot(111, projection=ccrs.PlateCarree(180))
    ax.set_aspect(1.2)
    ax.set_title(f'z500 Anomaly of type {num}', loc='left')
    ax.add_feature(cfeature.COASTLINE.with_scale('110m'), linewidth=0.5, edgecolor='dimgrey')
    ax.set_xticks(np.arange(0, 360, 60), crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(-90, 120, 30), crs=ccrs.PlateCarree())
    ax.xaxis.set_major_formatter(LongitudeFormatter())
    ax.yaxis.set_major_formatter(LatitudeFormatter())
    ax.tick_params(which='major', width=0.5, length=5)
    ax.spines['geo'].set_linewidth(0.5)
    cf = ax.contourf(f3.lon, f3.lat, ano, cmap=cmaps.BlueDarkRed18, levels=levels, extend='both', transform=ccrs.PlateCarree())
    cb = fig.colorbar(cf, orientation='horizontal', shrink=0.55, aspect=35, pad=0.12, ticks=levels)
    cb.outline.set_linewidth(0.1)
    cb.ax.tick_params(length=0)
    dot_area = np.where(p < 0.05)
    dot = ax.scatter(mesh_lon[dot_area], mesh_lat[dot_area], color='k', s=1.2, linewidths=0, transform=ccrs.PlateCarree())

z500_plot(ano=z_ano1, p=z_p1, num=1)
z500_plot(ano=z_ano2, p=z_p2, num=2)
z500_plot(ano=z_ano3, p=z_p3, num=3)

1271756664

终于找到源版了，感谢分享。

myj0219

谢谢分享！！！！！！！！！

kykaiii

谢谢谢谢！！！