使用python对海洋气象数据做显著性检验,并绘制空间pattern
选择数据集:
1 SST (Daily Sea Surface Temperature)
NOAA High-resolution Blended Analysis
- daily
- 分辨率:2.5
- 时间:2010
- http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v…
2 OLR (Outgoing Longwave Radiation)- NOAA Interpolated OLR - daily
- 分辨率:2.5
- 时间覆盖范围:1974-2013
使用编程工具:
- python
主要使用函数:
- stats.linregress() 函数说明点这里
内容:
- 对2010年的SST和OLR数据分布进行显著性检验,并绘制空间分布
- 先各自对数据计算一年的趋势以及相关,再计算两个数据之间的相关
- 绘制空间pattern
具体过程:
- 1 导入库
- 2 读取数据
- 3 计算相关和趋势
- 4 绘图
- 5 保存数据
过程还是比较清晰的,直接附上结果,
ps(标题时间打错了、横轴的label也搞错了,懒得重画了,代码中应该没问题了)
:
SST和OLR的pattern还是比较容易理解的,这两者的相关的空间分布,属实有点没看懂。
- 图中异常空白处是由于数据存在nan值导致的
- 包括打的点有些存在缺失,也和数据缺测有关联
最好,还是附上完整的代码吧,也没啥好保留的:
# -*- coding: utf-8 -*- """ Created on Sat Oct 22 20:50:09 2022 @author: Administrator """ import cartopy.feature as cfeature import numpy as np import xarray as xr from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter import cmaps import matplotlib.pyplot as plt import cartopy.crs as ccrs import matplotlib.ticker as mticker from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER from scipy import stats ############################################################################### olr_path = r'J:/olr.day.mean.nc' sst_path = r'J:/sst.intep.nc' da1 = xr.open_dataset(sst_path).sortby('lat') da2 = xr.open_dataset(olr_path).sortby('lat') sst = da1.sst.sel(lat=slice(-30,30),lon=slice(100,200)) olr = da2.olr.sel(lat=slice(-30,30),lon=slice(100,200), time=slice('2010','2010')) ############### calculate ################################################ trend = np.zeros((sst.lat.shape[0],sst.lon.shape[0])) p_value = np.zeros((sst.lat.shape[0],sst.lon.shape[0])) for i in range (0,sst.lat.shape[0]): for j in range (0,sst.lon.shape[0]): trend[i,j], intercept, r_value, p_value[i,j], std_err=stats.linregress(np.arange(1,366),sst[:,i,j]) ############################################################################## ################ plot ####################################################### lon = sst.lon.data lat = sst.lat.data ############################################################################## box = [100,200,-20,20] xstep,ystep = 20,10 proj = ccrs.PlateCarree(central_longitude=180) plt.rcParams['font.family'] = 'Times New Roman', ############################################################################## fig = plt.figure(figsize=(8,7),dpi=200) fig.tight_layout() ax = fig.add_axes([0.1,0.2,0.8,0.7],projection = proj) ax.set_extent(box,crs=ccrs.PlateCarree()) ax.coastlines('50m') ax.set_xticks(np.arange(box[0],box[1]+xstep, xstep),crs=ccrs.PlateCarree()) ax.set_yticks(np.arange(box[2], box[3]+1, ystep),crs=ccrs.PlateCarree()) lon_formatter = LongitudeFormatter(zero_direction_label=False)#True/False lat_formatter = LatitudeFormatter() ax.xaxis.set_major_formatter(lon_formatter) ax.yaxis.set_major_formatter(lat_formatter) ax.spines[['right','left','top','bottom']].set_linewidth(1.1) ax.spines[['right','left','top','bottom']].set_visible(True) ax.set_xlabel('Lon',fontsize=14) ax.set_title('Significance Test',fontsize=16,pad=8,loc='right') ax.set_title('this is title',fontsize=16,pad=8,loc='left') ax.tick_params( which='both',direction='in', width=0.7, pad=5, labelsize=14, bottom=True, left=True, right=True, top=True) c = ax.contourf(lon,lat,trend, # levels=np.linspace(-0.008,0.008,17), # levels=np.linspace(-0.32,0.32,17), #levels=np.linspace(-0.012,0.012,17), extend = 'both', transform=ccrs.PlateCarree(), cmap=cmaps.BlueWhiteOrangeRed) c1b = ax.contourf(lon,lat, p_value, [np.nanmin(p_value),0.05,np.nanmax(p_value)], hatches=['.', None], colors="none", transform=ccrs.PlateCarree()) cb=plt.colorbar(c, shrink=0.85, pad=0.15, orientation='horizontal', aspect=25, ) cb.ax.tick_params(labelsize=10,which='both',direction='in',) plt.show() # save picture # fig.savefig(r'D:\SST_OLR.png',format='png',dpi=500) ############################################################################## # xtick = np.arange(100, 200, 20) # ytick = np.arange(-30,31, 10) # gl=ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, # xlocs=[120,140,160,180,], # ylocs=[-30,-20,-10,0,10,20,30,], # x_inline=False, # y_inline=False, # linewidth=1.5, color='gray', alpha=0.5, linestyle='--', # ) # gl.xlabels_top = False # gl.ylabels_right = False # gl.xlines = True # plt.xlabel('Lon',fontsize=15) # plt.ylabel('Lat',fontsize=15) # gl.xformatter = LONGITUDE_FORMATTER # gl.yformatter = LATITUDE_FORMATTER