短期气候预测实习3

xiaoxiao啊

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


# 导入所需库
import numpy as np
import pandas as pd
import xarray as xr
from xskillscore import pearson_r, pearson_r_p_value
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


# 数据处理，根据公式计算EU指数
f1 = xr.open_dataset('hgt.mon.mean.nc')
z = f1.hgt.loc[f1.time.dt.month.isin([1]), 500].loc['1948':'2012']
z1 = z.loc[:, 55, 20].drop(['lat', 'lon'])
z2 = z.loc[:, 55, 75].drop(['lat', 'lon'])
z3 = z.loc[:, 40, 145].drop(['lat', 'lon'])
eu = -z1 / 4 + z2 / 2 - z3 / 4
eu_norm = (eu - eu.mean('time')) / eu.std('time')  # 标准化后的EU指数


# 绘画EU指数
def eu_plot():
    fig = plt.figure(figsize=(6, 3), facecolor='w')
    ax = fig.add_subplot(111)
    ax.tick_params(which='major', width=0.5, length=5)
    bars = ax.bar(range(1948, 2013), eu_norm, facecolor='#ff9999', edgecolor='r', width=1)
    for bar, height in zip(bars, eu_norm):
        if height < 0:
            bar.set(facecolor='#9999ff', edgecolor='blue')
    for i in ['bottom', 'top', 'left', 'right']:
        ax.spines【i】.set_linewidth(0.5)   # 记得改成英文的中括号，气象家园无法显示英文中括号加i
    ax.set_xlim(1947, 2013)
    ax.set_xlabel('year')
    ax.set_title('EU Index(1948-2012)', loc='left')
eu_plot()


# 将一维的EU指数转化为三维的
def to_3dims(*, var):
    var = np.array(var)
    var = np.expand_dims(var, axis=1)
    var = np.repeat(var, 73 * 144, axis=1)
    var = np.reshape(var, [65, 73, 144])
    var = xr.DataArray(var, coords=[z.time, z.lat, z.lon], dims=['time', 'lat', 'lon'])
    return var
eu_norm_3dims = to_3dims(var=eu_norm)  


# EU遥相关指数与同期环流场（500hPa高度场）的相关系数
r1 = pearson_r(eu_norm_3dims, z, dim='time')
p1 = pearson_r_p_value(eu_norm_3dims, z, dim='time')

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

# 绘图
def z500_plot():
    levels1 = np.arange(-1, 1.2, 0.2)
    levels2 = np.arange(-0.8, 1, 0.2)
    fig = plt.figure(figsize=(10, 6), facecolor='w')
    ax = fig.add_subplot(111, projection=ccrs.PlateCarree(180))
    ax.set_title('Corr EU & Z500 in Jan(1948-2012)', 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(f1.lon, f1.lat, r1, cmap=cmaps.BlueDarkRed18, transform=ccrs.PlateCarree(), levels=levels1)
    cb = fig.colorbar(cf, orientation='horizontal', shrink=0.6, aspect=25, pad=0.12, ticks=levels2)
    cb.outline.set_linewidth(0.1)
    cb.ax.tick_params(length=0)
    dot_area = np.where(p1 < 0.05)
    dot = ax.scatter(mesh_lon[dot_area], mesh_lat[dot_area], color='k', s=1.2, linewidths=0, transform=ccrs.PlateCarree())
    ax.set_aspect(1.2)
z500_plot()


df = pd.read_csv('station.csv', encoding='gbk', names=['station', 'num', 'lat', 'lon'], header=0)  # 读取160站站点数据
f2 = np.loadtxt('t1601.txt')  # 160站气温数据
t = f2.reshape(63, 160)[:-1]
t = xr.DataArray(t, coords=[z.time[3:], df.station], dims=['time', 'station'])


# 将一维的EU指数转化为二维的
def to_2dims_ch(*, var):
    var = np.array(var)
    var = np.expand_dims(var, axis=1)
    var = np.repeat(var, 160, axis=1)
    var = np.reshape(var, [62, 160])
    var = xr.DataArray(var, coords=[z.time[3:], df.station], dims=['time', 'station'])
    return var
eu_norm_2dims = to_2dims_ch(var=eu_norm[3:])


# EU遥相关指数与同期我国气温的相关系数
r2 = pearson_r(eu_norm_2dims, t, dim='time')
p2 = pearson_r_p_value(eu_norm_2dims, t, dim='time')


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)


# 站点数据插值为格点数据，三次样条插值
func = Rbf(df.lon, df.lat, r2, function='cubic')
grid_r2 = func(meshgrid_lon, meshgrid_lat)
func = Rbf(df.lon, df.lat, p2, function='cubic')
grid_p2 = func(meshgrid_lon, meshgrid_lat)


# 绘图
def temp_plot():
    levels1 = np.arange(-0.5, 0.6, 0.1)
    levels2 = np.arange(-0.4, 0.5, 0.1)
    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('Corr EU & Temp in Jan(1951-2012)', 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, grid_r2, cmap=cmaps.BlueDarkRed18, transform=ccrs.PlateCarree(), levels=levels1)
    dot_area = np.where(grid_p2 < 0.05)
    dot = ax1.scatter(meshgrid_lon[dot_area], meshgrid_lat[dot_area], color='k', s=1.2, linewidths=0, transform=ccrs.PlateCarree())
    china = get_map('中国')
    south_sea = get_map('南海')
    draw_map(china, color='dimgrey', linewidth=0.5)
    draw_map(south_sea, color='k', linewidth=0.6)
    clip_contours_by_map(cf, china)
    clip_pcolormesh_by_map(dot, china)
    ax1.add_geometries(Reader('河流.shp').geometries(), crs=ccrs.PlateCarree(), edgecolor='dimgrey', facecolor='none', lw=0.5)
    ax2 = fig.add_axes([0.768, 0.2005, 0.15, 0.15], projection=ccrs.PlateCarree())
    ax2.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5, edgecolor='dimgrey')
    ax2.set_extent([105, 124, 0, 23], ccrs.PlateCarree())
    ax2.spines['geo'].set_linewidth(0.5)
    draw_map(china, color='dimgrey', linewidth=0.5)
    draw_map(south_sea, color='k', linewidth=0.6)

    ax3 = fig.add_axes([0.15, 0.06, 0.7, 0.04])
    cb = fig.colorbar(cf, cax=ax3, orientation='horizontal', ticks=levels2)
    cb.outline.set_linewidth(0.1)
    cb.ax.tick_params(length=0)
temp_plot()

ZP盛夏光年

楼主，请问填加南海小地图怎么搞

来自安卓APP客户端

xiaoxiao啊

ZP盛夏光年 发表于 2022-6-17 11:33
楼主，请问填加南海小地图怎么搞

使用add_axes()添加子图

ZP盛夏光年

xiaoxiao啊 发表于 2022-6-17 16:45
使用add_axes()添加子图

好的谢谢

xiaoxiao啊

ZP盛夏光年 发表于 2022-6-18 09:24
好的谢谢

我现在把绘图代码也补上了

ZP盛夏光年

xiaoxiao啊 发表于 2022-6-18 15:23
我现在把绘图代码也补上了

感谢楼主，来学习啦

ZP盛夏光年

楼主，求求cnmaps包

来自安卓APP客户端

xiaoxiao啊

ZP盛夏光年 发表于 2022-6-26 17:22
楼主，求求cnmaps包

可以pip安装

ZP盛夏光年

xiaoxiao啊 发表于 2022-6-26 23:04
可以pip安装

我的pip安装报错了，可以分析一下吗

xiaoxiao啊

ZP盛夏光年 发表于 2022-6-27 22:16
我的pip安装报错了，可以分析一下吗

如果安装不上这个包的话那我也爱莫能助啦，因为我是能安装并使用的。如果用的是anaconda的话建议创建一个新的虚拟环境再安装试试，ps：不建议安装最新版本的