使用python对海洋气象数据做显著性检验,并绘制空间pattern

选择数据集: 1 SST (Daily Sea Surface Temperature) NOAA High-resolution Blended Analysis

daily分辨率:2.5时间:2010http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v… 2 OLR (Outgoing Longwave Radiation)- NOAA Interpolated OLRdaily分辨率:2.5时间覆盖范围:1974-2013http://www.esrl.noaa.gov/psd/data/gridded/data.interp_OLR.html

使用编程工具:

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

文章链接

评论可见,请评论后查看内容,谢谢!!!
 您阅读本篇文章共花了: