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

简介: 最近,师弟在学习使用python复现毕设论文,正好之前没有处理过站点数据,也没咋用过EOF,特此记录下使用python处理站的数据的过程。

前言



最近,师弟在学习使用python复现毕设论文,正好之前没有处理过站点数据,也没咋用过EOF,特此记录下使用python处理站的数据的过程。


  • 读取站点资料数据
  • 对站点数据进行插值,插值到规则网格上
  • 绘制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的绘图结果:


9fc7b76ff91e46f08d5cd45fd79c9943.png


使用metpt插值函数的结果:


image.png


附上全部的绘图代码:


# -*- 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的ocean er
              水平有限,欢迎指正!!!
              欢迎评论、收藏、点赞、转发、关注。
            关注我不后悔,记录学习进步的过程~~


相关文章
|
10天前
|
缓存 API 网络架构
淘宝item_search_similar - 搜索相似的商品API接口,用python返回数据
淘宝联盟开放平台中,可通过“物料优选接口”(taobao.tbk.dg.optimus.material)实现“搜索相似商品”功能。该接口支持根据商品 ID 获取相似推荐商品,并返回商品信息、价格、优惠等数据,适用于商品推荐、比价等场景。本文提供基于 Python 的实现示例,包含接口调用、数据解析及结果展示。使用时需配置淘宝联盟的 appkey、appsecret 和 adzone_id,并注意接口调用频率限制和使用规范。
|
2月前
|
存储 Web App开发 前端开发
Python + Requests库爬取动态Ajax分页数据
Python + Requests库爬取动态Ajax分页数据
|
12天前
|
JSON 安全 API
Python处理JSON数据的最佳实践:从基础到进阶的实用指南
JSON作为数据交换通用格式,广泛应用于Web开发与API交互。本文详解Python处理JSON的10个关键实践,涵盖序列化、复杂结构处理、性能优化与安全编程,助开发者高效应对各类JSON数据挑战。
70 1
|
2月前
|
数据采集 监控 调度
干货分享“用 多线程 爬取数据”:单线程 + 协程的效率反超 3 倍,这才是 Python 异步的正确打开方式
在 Python 爬虫中,多线程因 GIL 和切换开销效率低下,而协程通过用户态调度实现高并发,大幅提升爬取效率。本文详解协程原理、实战对比多线程性能,并提供最佳实践,助你掌握异步爬虫核心技术。
|
2月前
|
JSON 数据挖掘 API
闲鱼商品列表API响应数据python解析
闲鱼商品列表API(Goodfish.item_list)提供标准化数据接口,支持GET请求,返回商品标题、价格、图片、卖家信息等。适用于电商比价、数据分析,支持多语言调用,附Python示例代码,便于开发者快速集成。
|
2月前
|
JSON 自然语言处理 API
闲鱼商品详情API响应数据python解析
闲鱼商品详情API(goodfish.item_get)通过商品ID获取标题、价格、描述、图片等信息,支持Python等多语言调用。本文提供Python请求示例,包含请求构造与数据处理方法。
|
2月前
|
JSON API 数据格式
微店商品列表API响应数据python解析
微店商品列表API为开发者提供稳定高效获取商品信息的途径,支持HTTP GET/POST请求,返回JSON格式数据,含商品ID、名称、价格、库存等字段,适用于电商数据分析与展示平台搭建等场景。本文提供Python调用示例,助您快速上手。
|
6月前
|
机器学习/深度学习 存储 设计模式
Python 高级编程与实战:深入理解性能优化与调试技巧
本文深入探讨了Python的性能优化与调试技巧,涵盖profiling、caching、Cython等优化工具,以及pdb、logging、assert等调试方法。通过实战项目,如优化斐波那契数列计算和调试Web应用,帮助读者掌握这些技术,提升编程效率。附有进一步学习资源,助力读者深入学习。
|
3月前
|
Python
Python编程基石:整型、浮点、字符串与布尔值完全解读
本文介绍了Python中的四种基本数据类型:整型(int)、浮点型(float)、字符串(str)和布尔型(bool)。整型表示无大小限制的整数,支持各类运算;浮点型遵循IEEE 754标准,需注意精度问题;字符串是不可变序列,支持多种操作与方法;布尔型仅有True和False两个值,可与其他类型转换。掌握这些类型及其转换规则是Python编程的基础。
196 33
|
2月前
|
数据采集 分布式计算 大数据
不会Python,还敢说搞大数据?一文带你入门大数据编程的“硬核”真相
不会Python,还敢说搞大数据?一文带你入门大数据编程的“硬核”真相
71 1

推荐镜像

更多