python对站点数据做EOF且做插值绘制填 *
作者:oceanography-Rookie 发布时间:2023-03-05 03:30:56
前言
读取站点资料数据对站点数据进行插值,插值到规则网格上绘制EOF第一模态和第二模态的空间分布图绘制PC序列
关于插值,这里主要提供了两个插值函数,一个是一般常用的规则网格插值:
griddata
另一个是metpy中的:
inverse_distance_to_grid()
本来只是测验一下不同插值方法,但是发现两种插值方法的结果差别很大,由于对于站点数据处理较少,所以不太清楚具体原因。如果有知道的朋友可以告知一下,不甚感谢!
本次数据存储的文件格式为.txt
,读取的方法是通过pandas.read_csv()
同时,之前一直尝试使用proplot进行绘图,较长时间不用matplotlib.pyplot绘图了,也当做是复习一下绘图过程。
绘图中的代码主要通过封装函数,这样做的好处是大大减少了代码量。
导入必要的库:
from eofs.standard import Eof
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from metpy.interpolate import inverse_distance_to_grid
出现找不到库的报错,这里使用conda install packagename
安装一下就好
读取存储的数据
##################### read station data ##########################################
path = r'D:/data.txt'
file = pd.read_csv(path,sep= "\t",
header=None,
names=['station','lat','lon','year','data'],
na_values=-99.90)
data = file['data'].to_numpy()
lon = file['lon'].to_numpy()
lat = file['lat'].to_numpy()
year = file['year'].to_numpy()
station = file['station'].to_numpy()
year_max = np.max(year)
year_min = np.min(year)
year_range = np.arange(year_min,year_max+1,1)
data_all = data.reshape(70,53)
lon_all = lon.reshape(70,53)/100
lat_all = lat.reshape(70,53)/100
lon_real = lon_all[:,0]
lat_real = lat_all[:,0]
这里将读取的数据全部转为array
格式,方便查看以及后续处理。本来存储的文件中是没有相关的经度、纬度、站点、时间的名称的,这里我是自己加在上面方面数据读取的。
本次处理的数据包含70个站点,53年
插值
##################### interp data ##########################################
### interp data to target grid
### set target grid
lon_target = np.arange(115,135.5,0.5)
lat_target = np.arange(38,55.5,0.5)
x_t, y_t = np.meshgrid(lon_target, lat_target)
z = np.zeros((len(year_range),lat_target.shape[0],lon_target.shape[0]))
for i in range(len(year_range)):
print(i)
# z[i] = inverse_distance_to_grid(lon_real,lat_real,
# data_all[:,i],
# x_t,y_t, r=15, min_neighbors=3)
z[i] = griddata((lon_real,lat_real),
data_all[:,i],
(x_t,y_t),method='cubic')
这里显示了使用griddata()
的插值过程,metpy的过程注释掉了,需要测试的同学之间取消注释即可。
注意点:插值过程需要先设置目标的插值网格。
EOF处理
#计算纬度权重
lat_new = np.array(lat_target)
coslat=np.cos(np.deg2rad(lat_new))
wgts = np.sqrt(coslat)[..., np.newaxis]
#创建EOF分解器
solver=Eof(z,weights=wgts)
eof=solver.eofsAsCorrelation(neofs=2)
#此处的neofs的值是我们需要的空间模态数
pc=solver.pcs(npcs=2,pcscaling=1)#方差
var=solver.varianceFraction(neigs=2)
这里没啥好说的,需要几个模态自由选择即可
定义绘图函数并绘图
##################### def plot function ##########################################
def contourf(ax,i,level,cmap):
extents = [115,135,35,55]
ax.set_extent(extents, crs=proj)
ax.add_feature(cfeature.LAND, edgecolor='black',facecolor='silver',
)
ax.add_feature(cfeature.LAKES, edgecolor='black',facecolor='w',
)
ax.add_feature(cfeature.BORDERS)
xtick = np.arange(extents[0], extents[1], 5)
ytick = np.arange(extents[2], extents[3], 5)
ax.coastlines()
tick_proj = ccrs.PlateCarree()
c = ax.contourf(lon_target,lat_target,eof[i],
transform=ccrs.PlateCarree(),
levels=level,
extend='both',
cmap=cmap)
ax.set_xticks(xtick, crs=tick_proj)
ax.set_xticks(xtick, crs=tick_proj)
ax.set_yticks(ytick, crs=tick_proj)
ax.set_yticks(ytick, crs=tick_proj)
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())
plt.yticks(fontproperties='Times New Roman',size=10)
plt.xticks(fontproperties='Times New Roman',size=10)
ax.tick_params(which='major', direction='out',
length=4, width=0.5,
pad=5, bottom=True, left=True, right=True, top=True)
ax.tick_params(which='minor', direction='out',
bottom=True, left=True, right=True, top=True)
ax.set_title( 'EOF'+str(i),loc='left',fontsize =15)
return c
def p_line(ax,i):
ax.set_title('pc'+str(i)+'',loc='left',fontsize = 15)
# ax.set_ylim(-3.5,3.5)
ax.axhline(0,linestyle="--")
ax.plot(year_range,pc[:,i],color='blue')
ax.set_ylim(-3,3)
##################### plot ##########################################
fig = plt.figure(figsize=(8, 6), dpi=200)
proj = ccrs.PlateCarree()
contourf_1 = fig.add_axes([0.02,0.63,0.5,0.3],projection=proj)
c1=contourf(contourf_1,0,np.arange(0.7,1,0.05),plt.cm.bwr)
plot_1 = fig.add_axes([0.45,0.63,0.5,0.3])
p_line(plot_1,0)
contourf_2 = fig.add_axes([0.02,0.15,0.5,0.3],projection=proj)
c2= contourf(contourf_2,1,np.arange(-0.5,0.6,0.1),plt.cm.bwr)
plot_2 = fig.add_axes([0.45,0.15,0.5,0.3],)
p_line(plot_2,1)
cbposition1 = fig.add_axes([0.16, 0.55, 0.22, 0.02])
cb = fig.colorbar(c1,cax=cbposition1,
orientation='horizontal',format='%.1f')
cb.ax.tick_params(which='both',direction='in')
cbposition2=fig.add_axes([0.16, 0.08, 0.22, 0.02])
cb2 = fig.colorbar(c2,cax=cbposition2,
orientation='horizontal',format='%.1f')
cb2.ax.tick_params(which='both',direction='in')
plt.show()
这里将大部分重复的绘图代码,进行了封装,通过封装好的函数进行调用,大大简洁了代码量。相关的封装过程之前也有讲过,可以翻看之前的记录。
展示结果
使用griddata
的绘图结果:
使用metpt插值函数
的结果:
附上全部的绘图代码:
# -*- coding: utf-8 -*-
"""
Created on Fri Sep 23 17:46:42 2022
@author: Administrator
"""
from eofs.standard import Eof
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from metpy.interpolate import inverse_distance_to_grid
##################### read station data ##########################################
path = r'D:/data.txt'
file = pd.read_csv(path,sep= "\t",
header=None,
names=['station','lat','lon','year','data'],
na_values=-99.90)
data = file['data'].to_numpy()
lon = file['lon'].to_numpy()
lat = file['lat'].to_numpy()
year = file['year'].to_numpy()
station = file['station'].to_numpy()
year_max = np.max(year)
year_min = np.min(year)
year_range = np.arange(year_min,year_max+1,1)
data_all = data.reshape(70,53)
lon_all = lon.reshape(70,53)/100
lat_all = lat.reshape(70,53)/100
lon_real = lon_all[:,0]
lat_real = lat_all[:,0]
##################### interp data ##########################################
### interp data to target grid
### set target grid
lon_target = np.arange(115,135.5,0.5)
lat_target = np.arange(38,55.5,0.5)
x_t, y_t = np.meshgrid(lon_target, lat_target)
z = np.zeros((len(year_range),lat_target.shape[0],lon_target.shape[0]))
for i in range(len(year_range)):
print(i)
# z[i] = inverse_distance_to_grid(lon_real,lat_real,
# data_all[:,i],
# x_t,y_t, r=15, min_neighbors=3)
z[i] = griddata((lon_real,lat_real),
data_all[:,i],
(x_t,y_t),method='cubic')
#计算纬度权重
lat_new = np.array(lat_target)
coslat=np.cos(np.deg2rad(lat_new))
wgts = np.sqrt(coslat)[..., np.newaxis]
#创建EOF分解器
solver=Eof(z,weights=wgts)
eof=solver.eofsAsCorrelation(neofs=2)
#此处的neofs的值是我们需要的空间模态数
pc=solver.pcs(npcs=2,pcscaling=1)#方差
var=solver.varianceFraction(neigs=2)
##################### def plot function ##########################################
def contourf(ax,i,level,cmap):
extents = [115,135,35,55]
ax.set_extent(extents, crs=proj)
ax.add_feature(cfeature.LAND, edgecolor='black',facecolor='silver',
)
ax.add_feature(cfeature.LAKES, edgecolor='black',facecolor='w',
)
ax.add_feature(cfeature.BORDERS)
xtick = np.arange(extents[0], extents[1], 5)
ytick = np.arange(extents[2], extents[3], 5)
ax.coastlines()
tick_proj = ccrs.PlateCarree()
c = ax.contourf(lon_target,lat_target,eof[i],
transform=ccrs.PlateCarree(),
levels=level,
extend='both',
cmap=cmap)
ax.set_xticks(xtick, crs=tick_proj)
ax.set_xticks(xtick, crs=tick_proj)
ax.set_yticks(ytick, crs=tick_proj)
ax.set_yticks(ytick, crs=tick_proj)
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())
plt.yticks(fontproperties='Times New Roman',size=10)
plt.xticks(fontproperties='Times New Roman',size=10)
ax.tick_params(which='major', direction='out',
length=4, width=0.5,
pad=5, bottom=True, left=True, right=True, top=True)
ax.tick_params(which='minor', direction='out',
bottom=True, left=True, right=True, top=True)
ax.set_title( 'EOF'+str(i),loc='left',fontsize =15)
return c
def p_line(ax,i):
ax.set_title('pc'+str(i)+'',loc='left',fontsize = 15)
# ax.set_ylim(-3.5,3.5)
ax.axhline(0,linestyle="--")
ax.plot(year_range,pc[:,i],color='blue')
ax.set_ylim(-3,3)
##################### plot ##########################################
fig = plt.figure(figsize=(8, 6), dpi=200)
proj = ccrs.PlateCarree()
contourf_1 = fig.add_axes([0.02,0.63,0.5,0.3],projection=proj)
c1=contourf(contourf_1,0,np.arange(0.7,1,0.05),plt.cm.bwr)
plot_1 = fig.add_axes([0.45,0.63,0.5,0.3])
p_line(plot_1,0)
contourf_2 = fig.add_axes([0.02,0.15,0.5,0.3],projection=proj)
c2= contourf(contourf_2,1,np.arange(-0.5,0.6,0.1),plt.cm.bwr)
plot_2 = fig.add_axes([0.45,0.15,0.5,0.3],)
p_line(plot_2,1)
cbposition1 = fig.add_axes([0.16, 0.55, 0.22, 0.02])
cb = fig.colorbar(c1,cax=cbposition1,
orientation='horizontal',format='%.1f')
cb.ax.tick_params(which='both',direction='in')
cbposition2=fig.add_axes([0.16, 0.08, 0.22, 0.02])
cb2 = fig.colorbar(c2,cax=cbposition2,
orientation='horizontal',format='%.1f')
cb2.ax.tick_params(which='both',direction='in')
plt.show()
来源:https://blog.csdn.net/weixin_44237337/article/details/127018335
![](https://www.aspxhome.com/images/zang.png)
![](https://www.aspxhome.com/images/jiucuo.png)
猜你喜欢
- Linux终端里面可谓是奇妙无限,很多优秀的软件都诞生在终端里面。相较之下,Windows本身的理念和Linux就不一致,所以,你懂得。 下
- 前言随着我们不断地在一个文件中添加新的功能, 就会使得文件变得很长。 即便使用了继承,也抑制不住类的成长。为了解决这一问题,我们可以将类存储
- 获取一个类的所有子类def itersubclasses(cls, _seen=None): "
- 大量的多行段落本身就会降低可读性,同时空行分段也比空格分段有更高的可适应性...前文讨论的热烈程度远超我预期,正好还有之前查阅资料拍的几张实
- 今天有一位同学给了我一个excel文件,要求读取某些行,某些列,然后我试着做了一个demo,这里分享出来,希望能帮到大家:首先安装xlrd:
- 如下所示:import dateutildef before_month_lastday(ti): today=dateutil
- 如何制作一个分页程序?确实,翻页程序可以相互借鉴,但具体到每一需求,还是有较大差别的。代码入下,供参考:<%language=&quo
- 本文实例讲述了Python实现阿拉伯数字和罗马数字的互相转换功能。分享给大家供大家参考,具体如下:前面一篇介绍了《Java实现的求解经典罗马
- 测试代码1:def test(self): data = {"add": {"doc":
- 代码如下:DECLARE @c INT DECLARE @c2 INT SELECT @c = COUNT(1) FROM dbo.Spli
- 本文实例讲述了Python有序查找算法之二分法。分享给大家供大家参考,具体如下:二分法是一种快速查找的方法,时间复杂度低,逻辑简单易懂,总的
- 最近 W3C 一口气推出 7 个 HTML 工作草案,涵盖了 HTML5,HTML RDF,HTML Microdata,HTM
- docutils 的官方工具地址为:https://docutils.sourceforge.io/目前的更新主要是在版本和使用手册的更新上
- 可以查看: 代码如下:OPEN SYMMETRIC KEY 命令关于 对称密钥加密使用证书解密 CREATE MASTER KEY ENC
- 从Access数据库中选取记录有件最令人丧气的事情,它们是以怎样的顺序输入到数据库内就按照怎样的顺序出来。就算你在Access环境内采用So
- 这个问题用了我整整一晚上的时间才解决,希望有人遇到和我一样的时能少走些弯路。启动Django,服务器拒绝访问,可以尝试以下方法解决:1. 没
- 使用python爬虫库requests,urllib爬取今日头条街拍美图代码均有注释import re,json,requests,osfr
- zipfile是python里用来做zip格式编码的压缩和解压缩的,由于是很常见的zip格式,所以这个模块使用频率也是比较高的zipfile
- threading.Timer一次timer只生效一次,不会反复循环,如果实现循环触发,代码如下:import timeimport thr
- python提高图像质量概述调研了一些提高图像质量的方式深度学习方法,如微软的Bringing-Old-Photos-Back-to-Lif