python绘制回归分析图
本帖最后由 第八星系人气爱 于 2024-2-29 22:18 编辑作者：第八星系-Gugu
邮箱：gu2001yi@163.com
引入库并处理
import numpy as np
import xarray as xr
#import cmaps
import pandas as pd
#from eofs.standard import Eof
import matplotlib.pyplot as plt
from scipy import signal
from sklearn import preprocessing
#from scipy.stats import pearsonr
from scipy.stats import linregress
import datetime as dt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
from cartopy.util import add_cyclic_point
import matplotlib.patches as patches
g = xr.open_dataset('D:\\SSTdata\\sst.mnmean.nc')
print(g)
sst = g['sst'].loc)].loc['1954-01-01':'2003-09-01'].loc[:, 88:-88, 0:358]#选取时间及范围
sst_lat = g['lat'].loc
sst_lon = g['lon'].loc
sstx = np.array(sst).reshape(50,3, 89, 180)
sstx=0 #空值替换成0
sstx=np.array(sstx).mean(1)
ssty = signal.detrend(np.array(sstx).reshape(50,89, 180),axis=0,type='linear')+np.array(sstx).mean(0)#去趋势
ssti= np.zeros(50)
for i in range(0,50):
ssti=np.mean(ssty)#计算海温指数（此处计算的是北太平洋区域）
path = "D:\\NCEP\\NCEPdata\\hgt.mon.mean.nc" #读取位势高度场
f_z = xr.open_dataset(path)
z =np.array(f_z.hgt.loc)].loc['1954-05-01':'2003-09-01',500]).reshape(50,3,73,144).mean(1)
lat = f_z.lat
lon = f_z.lon
print(z.shape)
s,r,p = np.zeros((73,144)),np.zeros((73,144)),np.zeros((73,144))
for i in range(len(lat)):
for j in range(len(lon)):
s,_,r, p,_ = linregress(ssti,z[:,i,j])
绘图
fig = plt.figure(figsize=(12,10))
proj = ccrs.PlateCarree(central_longitude=180)
leftlon, rightlon, lowerlat, upperlat = (120,240,15,60)
img_extent =
ax = fig.add_axes(,projection = proj)
ax.set_extent(img_extent, crs=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE)
ax.set_xticks(np.arange(leftlon,rightlon+30,30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(lowerlat,upperlat+15,15), crs=ccrs.PlateCarree())
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_title('hgt(500hPa)-SSTI',loc='left',fontsize =14)
c1 = ax.contourf(lon,lat, s, zorder=0,levels =np.arange(-18,18,1) ,
extend = 'both', transform=ccrs.PlateCarree(), cmap=plt.cm.bwr)
c1b = ax.contourf(lon,lat, p,levels=, zorder=1,hatches=['...', None],colors="None", transform=ccrs.PlateCarree())
position=fig.add_axes()
fig.colorbar(c1,cax=position,orientation='horizontal',format='%.1f',)
plt.show()
图象：hgt海温回归指数
完整代码
#引入库
import numpy as np
import xarray as xr
#import cmaps
import pandas as pd
#from eofs.standard import Eof
import matplotlib.pyplot as plt
from scipy import signal
from sklearn import preprocessing
#from scipy.stats import pearsonr
from scipy.stats import linregress
import datetime as dt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
from cartopy.util import add_cyclic_point
import matplotlib.patches as patches
g = xr.open_dataset('D:\\SSTdata\\sst.mnmean.nc')
print(g)
sst = g['sst'].loc)].loc['1954-01-01':'2003-09-01'].loc[:, 88:-88, 0:358]#选取时间及范围
sst_lat = g['lat'].loc
sst_lon = g['lon'].loc
sstx = np.array(sst).reshape(50,3, 89, 180)
sstx=0 #空值替换成0
sstx=np.array(sstx).mean(1)
ssty = signal.detrend(np.array(sstx).reshape(50,89, 180),axis=0,type='linear')+np.array(sstx).mean(0)#去趋势
ssti= np.zeros(50)
for i in range(0,50):
ssti=np.mean(ssty)#计算海温指数（此处计算的是北太平洋区域）
path = "D:\\NCEP\\NCEPdata\\hgt.mon.mean.nc" #读取位势高度场
f_z = xr.open_dataset(path)
z =np.array(f_z.hgt.loc)].loc['1954-05-01':'2003-09-01',500]).reshape(50,3,73,144).mean(1)
lat = f_z.lat
lon = f_z.lon
print(z.shape)
s,r,p = np.zeros((73,144)),np.zeros((73,144)),np.zeros((73,144))
for i in range(len(lat)):
for j in range(len(lon)):
s,_,r, p,_ = linregress(ssti,z[:,i,j])
#绘图
fig = plt.figure(figsize=(12,10))
proj = ccrs.PlateCarree(central_longitude=180)
leftlon, rightlon, lowerlat, upperlat = (120,240,15,60)
img_extent =
ax = fig.add_axes(,projection = proj)
ax.set_extent(img_extent, crs=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE)
ax.set_xticks(np.arange(leftlon,rightlon+30,30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(lowerlat,upperlat+15,15), crs=ccrs.PlateCarree())
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_title('hgt(500hPa)-SSTI',loc='left',fontsize =14)
c1 = ax.contourf(lon,lat, s, zorder=0,levels =np.arange(-18,18,1) ,
extend = 'both', transform=ccrs.PlateCarree(), cmap=plt.cm.bwr)
c1b = ax.contourf(lon,lat, p,levels=, zorder=1,hatches=['...', None],colors="None", transform=ccrs.PlateCarree())
position=fig.add_axes()
fig.colorbar(c1,cax=position,orientation='horizontal',format='%.1f',)
plt.show()
微信搜索“第八星系人造大气理论爱好者”公众号，关注获取文章数据
页:
[1]