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 = 0y_fenmu = 0res = 0for j in range(len(x)):y_fenmu = y_fenmu + (y[j] - y_ave) ** 2y_fenmu = math.sqrt(y_fenmu / len(x))for j in range(len(x)):x_fenmu = x_fenmu + (x[j] - x_ave) ** 2x_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 maskoutdef calc_r(x, y):x_ave = np.mean(x)y_ave = np.mean(y)x_fenmu = 0y_fenmu = 0res = 0for j in range(len(x)):y_fenmu = y_fenmu + (y[j] - y_ave) ** 2y_fenmu = math.sqrt(y_fenmu / len(x))for j in range(len(x)):x_fenmu = x_fenmu + (x[j] - x_ave) ** 2x_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()

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

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

python气象绘图技巧之箱线图
python气象绘图技巧–seaborn catplot python气象绘图技巧--seaborn 前言一.构建dataframe和seaborn分类? 二.使用步骤 1.引入库 2.数据分析总 ...
python气象绘图速成_Python气象绘图Day-By-Day
登录后查看更多精彩内容~ 您需要登录才可以下载或查看,没有帐号?立即注册 x 本帖最后由 edwardli 于 2017-6-13 10:02 编辑工作繁事多,先上结果供交流,回头不断细化. 我 ...
【Python气象绘图临摹】处理数据（上）:读入输出nc数据、截取夏季/冬季数据、ButterWorth带通滤波、计算方差
前言 2022.9学习绘图利用python进行气象绘图,本文为学习绘制期间记录笔记,分为上.下两部分:处理数据和图像绘制.处理数据流程:读入olr资料,截取夏季/冬季数据,进行10-30dButte ...
【Python气象绘图临摹】图像绘制（下）:地理子图GeoAxes、xy轴设置、应用ncl色阶colormap、各标题、海岸线、添加文本、添加矩形框
文章目录前言 plt.fig.ax.三者绘制区别: GeoAxes地图投影:绘图投影和数据投影 x轴.y轴设置: 多个子图之间的间距调节: 图上添加文本.矩形框: python中ncl色阶color ...
Python气象绘图笔记——常用气象绘图函数脚本封装与使用记录
由于工作需要,将对我常用的python绘图脚本进行封装,为了防止代码丢失.忘记使用流程等,写个博客记录下. 要加载的包 import os import matplotlib.ticker as mt ...
python气象绘图速成_基于Python气象数据处理与可视化分析
基于 Python 气象数据处理与可视化分析张鑫 ; 曹蕾 ; 韩基良 [期刊名称] <气象灾害防御> [年 ( 卷 ), 期] 2020(027)001 [摘要] 全国综合气象信息共享 ...
python气象数据处理——计算台风方位角平均物理量
研究台风的同学们应该都接触过需要计算以台风为中心的方位角平均物理量,这就需要将笛卡尔坐标系中的数据插值到极坐标系,再对各个方位角的数据进行平均.最近在学习Python,相比NCL和FORTRAN,Py ...
python气象绘图相关链接
一页4图,绘制500风场及高度场.700的风场和相对湿度.850hpa的风场及高度场.海平面气压场.(3条消息) python学习记录-- 利用再分析数据绘制天气图_wdd_1992的博客-CSDN博 ...
python气象绘图_Python气象绘图教程（三）
更多的关于基础折线图技巧前面已经讲了很多关于折线图的常用参数,但是像颜色关键词在黑白文献中应该如何修改呢?plot()提供了一个marker=' '参数,其具体变化如下: plt.plot(x,te ...

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

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

效果图

代码实现

导入包

定义计算相关系数的函数

设置中文字体

读取文本数据（txt格式）

站点与相关系数汇总

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

导入中国地图

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

添加经纬度线

白化

展示

全部代码

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

最新文章

热门文章

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

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

效果图

代码实现

导入包

定义计算相关系数的函数

设置中文字体

读取文本数据（txt格式）

站点与相关系数汇总

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

导入中国地图

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

添加经纬度线

白化

展示

全部代码

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

最新文章

热门文章

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

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

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