气象数据处理与绘图_EOF
本帖最后由 月咏弥織 于 2020-8-13 09:08 编辑对于气象数据处理人而言，由于NCAR推出的NCL目前停更，取而代之的是更为普及的Python，因此将搜集来的各种j基于Python气象数据处理方式搬运并与大家分享。第一篇是数理中最为常见的EOF(经验正交分解)
本文将直接介绍该库的安装及使用，关于EOF的原理不做介绍。一、安装1
conda install -c conda-forge eofs
二、使用介绍
首先import
from eofs.standard import Eof
该库有几个基本函数是必须掌握的，我们一一介绍。1
solver = Eof(x, weights)
eof = solver.eofsAsCorrelation(neofs=3)
pc = solver.pcs(npcs=3, pcscaling=1)
var = solver.varianceFraction()
solver = Eof()建立一个EOF分解器，x为要进行分解的变量，weights为权重，通常指纬度权重solver.eofsAsCorrelation，solver.pcs，solver.varianceFraction分别取出空间模态，PC和方差。三、示例我们以中国夏季降水三类雨型的分解为例，展示EOF分析完整的Python实现。首先上图
数据下载：1961-2016夏季降水 以及 边界shp文件
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.mpl.ticker as cticker
import cartopy.io.shapereader as shpreader
import xarray as xr
from eofs.standard import Eof
#读取数据
f = xr.open_dataset('./pre.nc')
pre = np.array(f['pre'])
lat = f['lat']
lon = f['lon']
#计算纬度权重
lat = np.array(lat)
coslat = np.cos(np.deg2rad(lat))
wgts = np.sqrt(coslat)[..., np.newaxis]
#创建EOF分解器
solver = Eof(summer_mean_tmp, weights=wgts)
#获取前三个模态，获取对应的PC序列和解释方差
eof = solver.eofsAsCorrelation(neofs=3)
pc = solver.pcs(npcs=3, pcscaling=1)
var = solver.varianceFraction()以下是绘图部分，我们先给出其中不重复的部分，文章最末会给出完整代码：
#设置色标颜色，用于绘制PC柱状图
color1=[]
color2=[]
color3=[]
for i in range(1961,2017):
if pc >=0:
color1.append('red')
elif pc <0:
color1.append('blue')
if pc >=0:
color2.append('red')
elif pc <0:
color2.append('blue')
if pc >=0:
color3.append('red')
elif pc <0:
color3.append('blue')
#建立图形以及基本设置
fig = plt.figure(figsize=(15,15))
proj = ccrs.PlateCarree(central_longitude=115)
leftlon, rightlon, lowerlat, upperlat = (70,140,15,55)
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
#绘制第一模态
fig_ax1 = fig.add_axes(,projection = proj)
fig_ax1.set_extent(, crs=ccrs.PlateCarree())
fig_ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
fig_ax1.add_feature(cfeature.LAKES, alpha=0.5)
fig_ax1.set_xticks(np.arange(leftlon,rightlon+10,10), crs=ccrs.PlateCarree())
fig_ax1.set_yticks(np.arange(lowerlat,upperlat+10,10), crs=ccrs.PlateCarree())
fig_ax1.xaxis.set_major_formatter(lon_formatter)
fig_ax1.yaxis.set_major_formatter(lat_formatter)
china = shpreader.Reader('./bou2_4l.dbf').geometries()
fig_ax1.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
fig_ax1.set_title('(a) EOF1',loc='left',fontsize =15)
fig_ax1.set_title( '%.2f%%' % (var*100),loc='right',fontsize =15)
c1=fig_ax1.contourf(pre_lon,pre_lat, eof, levels=np.arange(-0.9,1.0,0.1), zorder=0, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
#添加南海
fig_ax11 = fig.add_axes(,projection = proj)
fig_ax11.set_extent(, crs=ccrs.PlateCarree())
fig_ax11.add_feature(cfeature.COASTLINE.with_scale('50m'))
china = shpreader.Reader('./bou2_4l.dbf').geometries()
fig_ax11.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
#添加色标
cbposition=fig.add_axes()
fig.colorbar(c1,cax=cbposition,orientation='horizontal',format='%.1f',)
#绘制PC
fig_ax4 = fig.add_axes()
fig_ax4.set_title('(b) PC1',loc='left',fontsize = 15)
fig_ax4.set_ylim(-2.5,2.5)
fig_ax4.axhline(0,linestyle="--")
fig_ax4.bar(np.arange(1961,2017,1),pc[:,0],color=color1)
plt.show()下面给出完整脚本的代码。看上去很长，但实际上大部分都是重复的：
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.mpl.ticker as cticker
import cartopy.io.shapereader as shpreader
import xarray as xr
from eofs.standard import Eof
f = xr.open_dataset('./pre.nc')
pre = np.array(f['pre'])
lat = f['lat']
lon = f['lon']
lat = np.array(lat)
coslat = np.cos(np.deg2rad(lat))
wgts = np.sqrt(coslat)[..., np.newaxis]
solver = Eof(pre, weights=wgts)
eof = solver.eofsAsCorrelation(neofs=3)
pc = solver.pcs(npcs=3, pcscaling=1)
var = solver.varianceFraction()
color1=[]
color2=[]
color3=[]
for i in range(1961,2017):
if pc >=0:
color1.append('red')
elif pc <0:
color1.append('blue')
if pc >=0:
color2.append('red')
elif pc <0:
color2.append('blue')
if pc >=0:
color3.append('red')
elif pc <0:
color3.append('blue')
fig = plt.figure(figsize=(15,15))
proj = ccrs.PlateCarree(central_longitude=115)
leftlon, rightlon, lowerlat, upperlat = (70,140,15,55)
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
fig_ax1 = fig.add_axes(,projection = proj)
fig_ax1.set_extent(, crs=ccrs.PlateCarree())
fig_ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
fig_ax1.add_feature(cfeature.LAKES, alpha=0.5)
fig_ax1.set_xticks(np.arange(leftlon,rightlon+10,10), crs=ccrs.PlateCarree())
fig_ax1.set_yticks(np.arange(lowerlat,upperlat+10,10), crs=ccrs.PlateCarree())
fig_ax1.xaxis.set_major_formatter(lon_formatter)
fig_ax1.yaxis.set_major_formatter(lat_formatter)
china = shpreader.Reader('./bou2_4l.dbf').geometries()
fig_ax1.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
fig_ax1.set_title('(a) EOF1',loc='left',fontsize =15)
fig_ax1.set_title( '%.2f%%' % (var*100),loc='right',fontsize =15)
c1=fig_ax1.contourf(pre_lon,pre_lat, eof, levels=np.arange(-0.9,1.0,0.1), zorder=0, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
fig_ax2 = fig.add_axes(,projection = proj)
fig_ax2.set_extent(, crs=ccrs.PlateCarree())
fig_ax2.add_feature(cfeature.COASTLINE.with_scale('50m'))
fig_ax2.add_feature(cfeature.LAKES, alpha=0.5)
fig_ax2.set_xticks(np.arange(leftlon,rightlon+10,10), crs=ccrs.PlateCarree())
fig_ax2.set_yticks(np.arange(lowerlat,upperlat+10,10), crs=ccrs.PlateCarree())
fig_ax2.xaxis.set_major_formatter(lon_formatter)
fig_ax2.yaxis.set_major_formatter(lat_formatter)
china = shpreader.Reader('./bou2_4l.dbf').geometries()
fig_ax2.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
fig_ax2.set_title('(c) EOF2',loc='left',fontsize =15)
fig_ax2.set_title( '%.2f%%' % (var*100),loc='right',fontsize =15)
c2=fig_ax2.contourf(pre_lon,pre_lat, eof, levels=np.arange(-0.9,1.0,0.1), zorder=0, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
fig_ax3 = fig.add_axes(,projection = proj)
fig_ax3.set_extent(, crs=ccrs.PlateCarree())
fig_ax3.add_feature(cfeature.COASTLINE.with_scale('50m'))
fig_ax3.add_feature(cfeature.LAKES, alpha=0.5)
fig_ax3.set_xticks(np.arange(leftlon,rightlon+10,10), crs=ccrs.PlateCarree())
fig_ax3.set_yticks(np.arange(lowerlat,upperlat+10,10), crs=ccrs.PlateCarree())
fig_ax3.xaxis.set_major_formatter(lon_formatter)
fig_ax3.yaxis.set_major_formatter(lat_formatter)
china = shpreader.Reader('./bou2_4l.dbf').geometries()
fig_ax3.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
fig_ax3.set_title('(e) EOF3',loc='left',fontsize =15)
fig_ax3.set_title( '%.2f%%' % (var*100),loc='right',fontsize =15)
c3=fig_ax3.contourf(pre_lon,pre_lat, eof, levels=np.arange(-0.9,1.0,0.1), zorder=0, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
fig_ax11 = fig.add_axes(,projection = proj)
fig_ax11.set_extent(, crs=ccrs.PlateCarree())
fig_ax11.add_feature(cfeature.COASTLINE.with_scale('50m'))
china = shpreader.Reader('./bou2_4l.dbf').geometries()
fig_ax11.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
fig_ax22 = fig.add_axes(,projection = proj)
fig_ax22.set_extent(, crs=ccrs.PlateCarree())
fig_ax22.add_feature(cfeature.COASTLINE.with_scale('50m'))
china = shpreader.Reader('./bou2_4l.dbf').geometries()
fig_ax22.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
fig_ax33 = fig.add_axes(,projection = proj)
fig_ax33.set_extent(, crs=ccrs.PlateCarree())
fig_ax33.add_feature(cfeature.COASTLINE.with_scale('50m'))
china = shpreader.Reader('./bou2_4l.dbf').geometries()
fig_ax33.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
cbposition=fig.add_axes()
fig.colorbar(c1,cax=cbposition,orientation='horizontal',format='%.1f',)
fig_ax4 = fig.add_axes()
fig_ax4.set_title('(b) PC1',loc='left',fontsize = 15)
fig_ax4.set_ylim(-2.5,2.5)
fig_ax4.axhline(0,linestyle="--")
fig_ax4.bar(np.arange(1961,2017,1),pc[:,0],color=color1)
fig_ax5 = fig.add_axes()
fig_ax5.set_title('(d) PC2',loc='left',fontsize = 15)
fig_ax5.set_ylim(-2.5,2.5)
fig_ax5.axhline(0,linestyle="--")
fig_ax5.bar(np.arange(1961,2017,1),pc[:,1],color=color2)
fig_ax6 = fig.add_axes()
fig_ax6.set_title('(f) PC3',loc='left',fontsize = 15)
fig_ax6.set_ylim(-2.5,2.5)
fig_ax6.axhline(0,linestyle="--")
fig_ax6.bar(np.arange(1961,2017,1),pc[:,2],color=color3)
plt.show() 学气象的，支持下，虽然看不懂 谢谢分享。有点难，小白慢慢学习。 正好工作需要，好好学习一下 楼主是搞气象的？请教个问题呢，针对grib文件你是用的哪种方式读取？ 谢谢分享！
页:
[1]