python 气象绘图 计算多年一月的nino指数和站点降水相关图

python 气象绘图 计算多年一月的nino指数和站点降水相关图

效果图

代码实现

导入包

import cartopy.crs as ccrs
import cartopy.io.shapereader as shapereader
import matplotlib.pyplot as plt
import numpy as np
import math
import pandas as pd
from scipy.interpolate import Rbf
import maskout

定义计算相关系数的函数

def calc_r(x, y):
    x_ave = np.mean(x)
    y_ave = np.mean(y)
    x_fenmu = 0
    y_fenmu = 0
    res = 0
    for j in range(len(x)):
        y_fenmu = y_fenmu + (y[j] - y_ave) ** 2
    y_fenmu = math.sqrt(y_fenmu / len(x))
    for j in range(len(x)):
        x_fenmu = x_fenmu + (x[j] - x_ave) ** 2
    x_fenmu = math.sqrt(x_fenmu / len(x))
    for j in range(len(x)):
        res = res + (x[j] - x_ave) * (y[j] - y_ave) / (x_fenmu * y_fenmu)
    return res / len(x)

设置中文字体

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

读取文本数据（txt格式）

# 读取nino
file1 = r'nino34.txt'
f1 = open(file1).read()
nino = np.array(f1.split()).astype(float).reshape(63, 13)

# 读取站点数据
file2 = r'r1607.txt'
f2 = open(file2).read()
zhandian = np.array(f2.split()).astype(float).reshape(160, 63)

# 读取站点位置
file3 = r'id.txt'
f3 = open(file3).read()
loc = np.array(f3.split()).astype(float).reshape(160, 3)

因为数据没有缺测，所以直接read 然后spilt ，若有缺测可以用 nino[nino==缺测值] = nan

站点与相关系数汇总

# 站点与相关系数汇总
loc = pd.DataFrame(loc, columns=['站点号', 'lat', 'lon'])
r = pd.DataFrame(r, columns=['r'])
data = pd.concat([loc, r], axis=1)
# print(data.shape)

将数据转成DataFrame格式（或许可以直接用pandas读，笔者很菜）

设置经纬度并用scipy包中的函数插值

lon = data['lon']
lat = data['lat']
r = data['r']
olon = np.linspace(70, 140, 100)
olat = np.linspace(0, 60, 100)
olon, olat = np.meshgrid(olon, olat)
# 插值处理
func = Rbf(lon, lat, r, function='linear')
data_new = func(olon, olat)

导入中国地图

fig = plt.figure(figsize=(16, 12))
proj = ccrs.PlateCarree()
ax = fig.add_subplot(1, 1, 1, projection=proj)
china = shapereader.Reader(r'1\provinces.shp').geometries()
china_nine_dotted_line = shapereader.Reader(r'1\china_nine_dotted_line.shp').geometries()
ax.add_geometries(china_nine_dotted_line, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
ax.add_geometries(china, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
ax.set_extent([70, 140, 0, 55])

绘制填色图并表明站点的值

x, y = (olon, olat)
xx, yy = (lon, lat)
levels = np.linspace(-0.4, 0.4, 50)
cf = plt.contourf(x, y, data_new, levels=levels, cmap='RdBu')
sc = ax.scatter(xx - 0.1, yy, c='k', s=5, marker='o')
cbar = plt.colorbar(cf, ticks=np.linspace(-0.4, 0.4, 11))
for i in range(0, len(xx)):
    plt.text(xx[i], yy[i], '{:.2f}'.format(data['r'][i]), va='center', fontsize=8)

添加经纬度线

gl = ax.gridlines(draw_labels=True, linestyle=":", linewidth=0.3, color='k')

白化

clip = maskout.shp2clip(cf, ax, r'1\china0.shp')

展示

plt.title('各站点和nino一月指数的相关系数')
plt.show()

全部代码

import cartopy.crs as ccrs
import cartopy.io.shapereader as shapereader
import matplotlib.pyplot as plt
import numpy as np
import math
import pandas as pd
from scipy.interpolate import Rbf
import maskout


def calc_r(x, y):
    x_ave = np.mean(x)
    y_ave = np.mean(y)
    x_fenmu = 0
    y_fenmu = 0
    res = 0
    for j in range(len(x)):
        y_fenmu = y_fenmu + (y[j] - y_ave) ** 2
    y_fenmu = math.sqrt(y_fenmu / len(x))
    for j in range(len(x)):
        x_fenmu = x_fenmu + (x[j] - x_ave) ** 2
    x_fenmu = math.sqrt(x_fenmu / len(x))
    for j in range(len(x)):
        res = res + (x[j] - x_ave) * (y[j] - y_ave) / (x_fenmu * y_fenmu)
    return res / len(x)


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

# 读取数据
# 读取nino
file1 = r'nino34.txt'
f1 = open(file1).read()
nino = np.array(f1.split()).astype(float).reshape(63, 13)

# 读取站点数据
file2 = r'r1607.txt'
f2 = open(file2).read()
zhandian = np.array(f2.split()).astype(float).reshape(160, 63)

# 读取站点位置
file3 = r'id.txt'
f3 = open(file3).read()
loc = np.array(f3.split()).astype(float).reshape(160, 3)

# 计算相关系数
r = np.array([])
for i in range(160):
    r = np.append(r, calc_r(x=zhandian[i, :], y=nino[:, 1]))

# 站点与相关系数汇总
loc = pd.DataFrame(loc, columns=['站点号', 'lat', 'lon'])
r = pd.DataFrame(r, columns=['r'])
data = pd.concat([loc, r], axis=1)
# print(data.shape)

lon = data['lon']
lat = data['lat']
r = data['r']
olon = np.linspace(70, 140, 100)
olat = np.linspace(0, 60, 100)
olon, olat = np.meshgrid(olon, olat)
# 插值处理
func = Rbf(lon, lat, r, function='linear')
data_new = func(olon, olat)

# 绘制地图
fig = plt.figure(figsize=(16, 12))
proj = ccrs.PlateCarree()
ax = fig.add_subplot(1, 1, 1, projection=proj)
china = shapereader.Reader(r'1\provinces.shp').geometries()
china_nine_dotted_line = shapereader.Reader(r'1\china_nine_dotted_line.shp').geometries()
ax.add_geometries(china_nine_dotted_line, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
ax.add_geometries(china, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
ax.set_extent([70, 140, 0, 55])
x, y = (olon, olat)
xx, yy = (lon, lat)
levels = np.linspace(-0.4, 0.4, 50)
cf = plt.contourf(x, y, data_new, levels=levels, cmap='RdBu')
sc = ax.scatter(xx - 0.1, yy, c='k', s=5, marker='o')
cbar = plt.colorbar(cf, ticks=np.linspace(-0.4, 0.4, 11))
for i in range(0, len(xx)):
    plt.text(xx[i], yy[i], '{:.2f}'.format(data['r'][i]), va='center', fontsize=8)
gl = ax.gridlines(draw_labels=True, linestyle=":", linewidth=0.3, color='k')
clip = maskout.shp2clip(cf, ax, r'1\china0.shp')
plt.title('各站点和nino一月指数的相关系数')
plt.show()

（如有建议或问题，欢迎讨论）