利用Python+xarray实现遥感数据 |
您所在的位置：网站首页 › EOF气象 › 利用Python+xarray实现遥感数据 |
利用Python+xarray实现遥感数据
|
利用python实现气象遥感数据的EOF分解，研究其空间变化特征。EOF的原理不再细讲。 1. 安装eofs conda install eofs or pip install eofs 2. 进行EOF分解 2.1 import 模块 import xarray as xr import numpy as np from eofs.standard import Eof import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER,mticker import cartopy.mpl.ticker as cticker 2.2 EOF的实现 SST= xr.open_dataset(r"*****mean_16.nc")['sst'].fillna(0.0) cc = np.array(SST) print(cc.shape) lat = SST['lat'] lon = SST['lon'] coslat = np.cos(np.deg2rad(lat)) wgts = np.sqrt(coslat)[..., np.newaxis] solver = Eof(cc, weights=wgts) eof = solver.eofsAsCorrelation(neofs=3) # spatial pc = solver.pcs(npcs=3, pcscaling=1) # temporal var = solver.varianceFraction() # 解释方差 2.3 制图 fig_ax1 = fig.add_axes([0.1, 0.8, 0.5, 0.3],projection = proj) fig_ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree()) # 设置地图属性:加载国界、海岸线、 fig_ax1.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.8, zorder=1) fig_ax1.add_feature(cfeature.BORDERS.with_scale('50m'), linewidth=0.8, zorder=1) fig_ax1.set_xticks(np.arange(leftlon,rightlon+5,10), crs=ccrs.PlateCarree()) fig_ax1.set_yticks(np.arange(lowerlat,upperlat+5,10), crs=ccrs.PlateCarree()) fig_ax1.xaxis.set_major_formatter(lon_formatter) fig_ax1.yaxis.set_major_formatter(lat_formatter) fig_ax1.set_title('(a) EOF1',loc='left',fontsize =15,weight = 'semibold') fig_ax1.set_title( '%.2f%%' % (var[0]*100),loc='right',fontsize =15,weight = 'semibold') c1 = data0[0].plot(ax = fig_ax1, transform =ccrs.PlateCarree(), cmap = plt.cm.get_cmap('RdBu_r'), add_colorbar = False, ) fig_ax2 = fig.add_axes([0.1, 0.45, 0.5, 0.3],projection = proj) fig_ax2.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree()) #fig_ax2.add_feature(cfeature.OCEAN, zorder=0,color = 'white',edgecolor='black') fig_ax2.add_feature(cfeature.LAND, color = 'lightgray',zorder=1,edgecolor='black') fig_ax2.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.8, zorder=1) fig_ax2.add_feature(cfeature.BORDERS.with_scale('50m'), linewidth=0.8, zorder=1) fig_ax2.set_xlabel('x',fontweight ='bold') fig_ax2.set_ylabel('y',fontweight ='bold') fig_ax2.set_xticks(np.arange(leftlon,rightlon+5,10), crs=ccrs.PlateCarree()) fig_ax2.set_yticks(np.arange(lowerlat,upperlat+5,10), crs=ccrs.PlateCarree()) fig_ax2.xaxis.set_major_formatter(lon_formatter) fig_ax2.yaxis.set_major_formatter(lat_formatter) fig_ax2.set_title('(c) EOF2',loc='left',fontsize =15,weight = 'semibold') fig_ax2.set_title( '%.2f%%' % (var[1]*100),loc='right',fontsize =15,weight = 'semibold') c2 = data0[1].plot(ax = fig_ax2, transform =ccrs.PlateCarree(), cmap = plt.cm.get_cmap('RdBu_r'), add_colorbar = False, ) fig_ax3 = fig.add_axes([0.1, 0.1, 0.5, 0.3],projection = proj) fig_ax3.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree()) fig_ax3.add_feature(cfeature.LAND, color = 'lightgray',zorder=1,edgecolor='black') fig_ax3.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.8, zorder=1) fig_ax3.add_feature(cfeature.BORDERS.with_scale('50m'), linewidth=0.8, zorder=1) fig_ax3.set_xlabel('x ',fontweight ='bold') fig_ax3.set_ylabel('y ',fontweight ='bold') fig_ax3.set_xticks(np.arange(leftlon,rightlon+5,10), crs=ccrs.PlateCarree()) fig_ax3.set_yticks(np.arange(lowerlat,upperlat+5,10), crs=ccrs.PlateCarree()) fig_ax3.xaxis.set_major_formatter(lon_formatter) fig_ax3.yaxis.set_major_formatter(lat_formatter) fig_ax3.set_title('(e) EOF3',loc='left',fontsize =15,weight = 'semibold') fig_ax3.set_title( '%.2f%%' % (var[2]*100),loc='right',fontsize =15,weight = 'semibold') c3 = data0[2].plot(ax = fig_ax3, transform =ccrs.PlateCarree(), cmap = plt.cm.get_cmap('RdBu_r'), add_colorbar = False, ) cbposition=fig.add_axes([0.13, 0.04, 0.4, 0.015]) fig.colorbar(c1,cax=cbposition,orientation='horizontal',format='%.1f',) fig_ax4 = fig.add_axes([0.55, 0.808, 0.40, 0.285]) fig_ax4.set_title('(b) PC1',loc='left',fontsize = 15,weight = 'semibold') fig_ax4.set_ylim(-3.0,3.0) fig_ax4.axhline(0,linestyle="--") fig_ax4.plot(np.arange(1,len(cc)+1,1),pc[:,0]) fig_ax5 = fig.add_axes([0.55, 0.458, 0.40, 0.285]) fig_ax5.set_title('(d) PC2',loc='left',fontsize = 15,weight = 'semibold') fig_ax5.set_ylim(-3.0,3.0) fig_ax5.axhline(0,linestyle="--") fig_ax5.plot(np.arange(1,len(cc)+1,1),pc[:,1]) fig_ax6 = fig.add_axes([0.55, 0.108, 0.40, 0.285]) fig_ax6.set_title('(f) PC3',loc='left',fontsize = 15,weight = 'semibold') fig_ax6.set_ylim(-3.0,3) fig_ax6.axhline(0,linestyle="--") fig_ax6.plot(np.arange(1,len(cc)+1,1),pc[:,2]) #plt.tight_layout() plt.show() 2.4 结果展示
前面的EOF分析还有待完善，结果图的详细分析后续更新… |
【本文地址】