Python | 海表面温度(SST) | 长期趋势和异常分析

06-11 1130阅读

趋势和异常分析（Trend and anomaly)在大气和海洋学研究中被广泛用于探测长期变化。

趋势分析（Trend Analysis）：

趋势分析是一种用于检测数据随时间的变化趋势的方法。在海洋学和大气学中，常见的趋势分析包括海表面温度（SST）、海平面上升、气温变化等。趋势分析通常包括以下步骤：

数据预处理：首先需要对数据进行预处理，包括去除季节循环、填补缺失值等。
计算趋势：采用统计方法（如线性回归、非线性回归）来计算数据随时间的变化趋势。常见的趋势计算方法包括最小二乘法、曲线拟合等。
趋势显著性检验：对计算得到的趋势进行显著性检验，以确定趋势是否具有统计显著性。常见的显著性检验方法包括t检验、F检验等。

趋势可视化：将计算得到的趋势以图形方式呈现，通常使用折线图或柱状图来展示数据随时间的变化趋势。

趋势分析的结果可以帮助科学家们了解气候系统的长期演变趋势，从而预测未来可能的变化情况。

异常分析（Anomaly Analysis）：

异常分析是一种用于检测数据中非正常事件或突发事件的方法。在海洋学和大气学中，异常分析通常用于检测气候系统中的异常事件，如El Niño事件、极端气候事件等。异常分析通常包括以下步骤：

基准确定：选择一个合适的基准期，通常是一段相对稳定的时间段，用于计算异常值。
计算异常：将观测数据与基准期的平均值进行比较，计算出每个时间点的异常值。异常值表示该时间点的数据与基准期相比的偏离程度。
异常检测：对计算得到的异常值进行检测，识别出突发事件或非正常事件。

异常可视化：将计算得到的异常值以图形方式呈现，通常使用折线图或柱状图来展示异常事件的发生情况。

异常分析的结果可以帮助科学家们理解气候系统中的非正常事件，从而采取相应的应对措施或预测未来可能发生的异常情况。

本案例分析以海表温度为例，计算了1982年至2016年全球每十年的温度变化率。此外，还给出了其面积加权的全球月海温异常时间序列。

数据来源：
NOAA Optimum Interpolation (OI) Sea Surface Temperature (SST) V2

下载地址：

https://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.html

NOAA Optimum Interpolation (OI) Sea Surface Temperature (SST) V2

空间分辨率:

1.0纬度 x 1.0经度全球网格(180x360)。

覆盖范围

89.5 N-89.5 S 
0.5 E-359.5 E.

因为oisst是一个插值数据，所以它覆盖了海洋和陆地。

因此，必须同时使用陆地-海洋掩膜数据，可以从如下网站获得：

https://psl.noaa.gov/repository/entry/show?entryid=b5492d1c-7d9c-47f7-b058-e84030622bbd

1. 加载基础库

import numpy as np
import datetime 
import cftime
from netCDF4 import Dataset as netcdf # netcdf4-python module
import netCDF4 as nc
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import matplotlib.dates as mdates
from matplotlib.dates import MonthLocator, WeekdayLocator, DateFormatter
import matplotlib.ticker as ticker
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6
import warnings
warnings.simplefilter('ignore')
# Your code continues here
# ================================================================================================
# Author: %(Jianpu)s | Affiliation: Hohai
# Email : %(email)s
# Last modified: 2024-04-04 12:28:06
# Filename: Trend and Anomaly Analyses.py
# Description:
# =================================================================================================

2. 读取数据、提取变量

2.1 Read SST

ncset= netcdf(r'G:/code_daily/sst.mnmean.nc')
lons = ncset['lon'][:]  
lats = ncset['lat'][:]           
sst  = ncset['sst'][1:421,:,:]    # 1982-2016 to make it divisible by 12
nctime = ncset['time'][1:421]
t_unit = ncset['time'].units
try :
    t_cal =ncset['time'].calendar
except AttributeError : # Attribute doesn't exist
    t_cal = u"gregorian" # or standard
nt, nlat, nlon = sst.shape
ngrd = nlon*nlat

2.2 解析时间

datevar = nc.num2date(nctime,units = "days since 1800-1-1 00:00:00")
print(datevar.shape)
datevar[0:5]

2.3读取 mask (1=ocen, 0=land)

lmfile = 'G:/code_daily/lsmask.nc'
lmset  = netcdf(lmfile)
lsmask = lmset['mask'][0,:,:]
lsmask = lsmask-1
num_repeats = nt
lsm = np.stack([lsmask]*num_repeats,axis=-1).transpose((2,0,1))
lsm.shape

2.3 将海温的陆地区域进行掩膜

sst = np.ma.masked_array(sst, mask=lsm)
3. Trend Analysis
3.1 计算线性趋势
sst_grd  = sst.reshape((nt, ngrd), order='F') 
x        = np.linspace(1,nt,nt)#.reshape((nt,1))
sst_rate = np.empty((ngrd,1))
sst_rate[:,:] = np.nan
for i in range(ngrd): 
    y = sst_grd[:,i]   
    if(not np.ma.is_masked(y)):         
        z = np.polyfit(x, y, 1)
        sst_rate[i,0] = z[0]*120.0
        #slope, intercept, r_value, p_value, std_err = stats.linregress(x, sst_grd[:,i])
        #sst_rate[i,0] = slope*120.0     
    
sst_rate = sst_rate.reshape((nlat,nlon), order='F')

3 绘制趋势空间分布

plt.figure(dpi=200)
m = Basemap(projection='cyl', llcrnrlon=min(lons), llcrnrlat=min(lats),
        urcrnrlon=max(lons), urcrnrlat=max(lats))
x, y = m(*np.meshgrid(lons, lats))
clevs = np.linspace(-0.5, 0.5, 21)
cs = m.contourf(x, y, sst_rate.squeeze(), clevs, cmap=plt.cm.RdBu_r)
m.drawcoastlines()
#m.fillcontinents(color='#000000',lake_color='#99ffff')
cb = m.colorbar(cs)
cb.set_label('SST Changing Rate ($^oC$/decade)', fontsize=12)
plt.title('SST Changing Rate ($^oC$/decade)', fontsize=16)

4. Anomaly analysis

4.1 转换数据大小为： (nyear) x (12) x (lat x lon)

sst_grd_ym  = sst.reshape((12,round(nt/12), ngrd), order='F').transpose((1,0,2))
sst_grd_ym.shape

4.2 计算季节趋势

sst_grd_clm = np.mean(sst_grd_ym, axis=0)
sst_grd_clm.shape

4.3 去除季节趋势

sst_grd_anom = (sst_grd_ym - sst_grd_clm).transpose((1,0,2)).reshape((nt, nlat, nlon), order='F')
sst_grd_anom.shape

4.4 计算区域权重

4.4.1 确认经纬度的方向

print(lats[0:12])
print(lons[0:12])

4.4.2 计算随纬度变化的区域权重

lonx, latx = np.meshgrid(lons, lats)
weights = np.cos(latx * np.pi / 180.)
print(weights.shape)

4.4.3 计算全球、北半球、南半球的有效网格总面积

sst_glb_avg = np.zeros(nt)
sst_nh_avg  = np.zeros(nt)
sst_sh_avg  = np.zeros(nt)
for it in np.arange(nt):
    sst_glb_avg[it] = np.ma.average(sst_grd_anom[it, :], weights=weights)
    sst_nh_avg[it]  = np.ma.average(sst_grd_anom[it,0:round(nlat/2),:],    weights=weights[0:round(nlat/2),:])
    sst_sh_avg[it]  = np.ma.average(sst_grd_anom[it,round(nlat/2):nlat,:], weights=weights[round(nlat/2):nlat,:])

4.5 转换时间为字符串格式

datestr = [date.strftime('%Y-%m-%d') for date in datevar]

5. 绘制海温异常时间序列

fig, ax = plt.subplots(1, 1 , figsize=(15,5),dpi=200)
ax.plot(datestr[::12], sst_glb_avg[::12], color='b', linewidth=2, label='GLB')
ax.plot(datestr[::12], sst_nh_avg[::12],  color='r', linewidth=2, label='NH')
ax.plot(datestr[::12], sst_sh_avg[::12],  color='g', linewidth=2, label='SH')
ax.set_xticklabels(datestr[::12], rotation=45)
ax.axhline(0, linewidth=1, color='k')
ax.legend()
ax.set_title('Monthly SST Anomaly Time Series (1982 - 2016)', fontsize=16)
ax.set_xlabel('Month/Year ', fontsize=12)
ax.set_ylabel('$^oC$', fontsize=12)
ax.set_ylim(-0.6, 0.6)
fig.set_figheight(9)
# rotate and align the tick labels so they look better
fig.autofmt_xdate()
# use a more precise date string for the x axis locations in the toolbar
ax.fmt_xdata = mdates.DateFormatter('%Y')

http://unidata.github.io/netcdf4-python/

http://www.scipy.org/

本文由mdnice多平台发布

VPS购买请点击我