python对站点数据做EOF且做插值绘制填色图

更新时间：2022年09月26日 11:46:22 作者：oceanography-Rookie

这篇文章主要介绍了python对站点数据做EOF且做插值绘制填色图，文章围绕主题展开详细的内容介绍，具有一定的参考价值，，需要的小伙伴可以参考一下

前言

读取站点资料数据对站点数据进行插值，插值到规则网格上绘制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()

总结

metpy的插值函数好处在于可以自由填充整个绘图区域，但是感觉griddata函数的插值结果更加符合预期，虽然也有点怪怪的。

到此这篇关于python对站点数据做EOF且做插值绘制填色图的文章就介绍到这了,更多相关python EOF插值绘制内容请搜索脚本之家以前的文章或继续浏览下面的相关文章希望大家以后多多支持脚本之家！

您可能感兴趣的文章:

python3实现全角和半角字符转换的方法示例
在自然语言处理过程中，全角、半角的的不一致会导致信息抽取不一致，因此需要统一，下面这篇文章主要给大家介绍了关于python3中全角和半角字符转换的方法，需要的朋友可以参考借鉴，下面来一起看看吧。
2017-09-09
详解Python中import模块导入的实现原理
这篇文章主要给大家介绍了Python中import模块导入的实现原理,主要从什么是模块,import搜索路径以及导入原理这三个方面给大家介绍,感兴趣的小伙伴跟着小编一起来看看吧
2023-08-08
Python Requests 基础入门
首先，Python 标准库中的 urllib2 模块提供了你所需要的大多数 HTTP 功能，但是它的 API 不友好。通过本文给大家介绍Python Requests 入门基础，感兴趣的朋友一起学习吧
2016-04-04
解决Python报错:SyntaxError: invalid character ‘，‘ (U+FF0C)
Python中的 SyntaxError错误是Python语言中常见的异常错误类型之一,表示语法错误,下面这篇文章主要给大家介绍了关于解决Python报错:SyntaxError: invalid character ‘，‘ (U+FF0C)的相关资料,需要的朋友可以参考下
2022-12-12
python对象与json相互转换的方法
这篇文章主要介绍了python对象与json相互转换，文中通过示例代码介绍的非常详细，对大家的学习或者工作具有一定的参考学习价值，需要的朋友们下面随着小编来一起学习学习吧
2019-05-05
python统计指定目录内文件的代码行数
这篇文章主要介绍了python统计指定目录内文件的代码行数
2019-09-09
python如何利用traceback获取详细的异常信息
异常信息的获取对于程序的调试非常重要，可以有助于快速定位有错误程序语句的位置。这篇文章主要给大家介绍了关于python如何利用traceback获取详细的异常信息的相关资料,需要的朋友可以参考下
2021-06-06
Python如何执行系统命令
这篇文章主要介绍了Python如何执行系统命令，帮助大家更好的利用python操作系统，感兴趣的朋友可以了解下
2020-09-09
Python中读取和加解密PDF文件的详细教程
在Python中读取和加密PDF文件是一项常见且实用的任务,尤其对于需要处理大量文档自动化处理的场景,本文将详细介绍如何使用Python读取PDF文件内容以及如何使用不同的库来给PDF文件加密,需要的朋友可以参考下
2024-08-08
python中以函数作为参数(回调函数)的实现方法
这篇文章主要介绍了python中以函数作为参数(回调函数)的实现方法，具有很好的参考价值，希望对大家有所帮助。如有错误或未考虑完全的地方，望不吝赐教
2023-01-01