写在前面
和之前类似的内容,只不过换了一下数据为GPCP月平均资料。主要还是简单的趋势分析以及气候态空间分布,将绘制方法这里调整为basemap,也算是温习一下basemap的绘图函数
GPCP(月平均资料)
对 1979 年至 2023 年 月平均 GCPCP 降水数据进行了初步分析:
- 趋势
- 气候态年平均降水量
- 降水率的面积加权平均值
- 纬向平均降水量
数据
月平均 GPCP(全球降水气候学项目)数据融合了测量站、卫星和探空观测的降雨量,以估计 1979 年至今 2.5 度全球网格上的月降雨量。
数据地址为
- https://climatedataguide.ucar.edu/climate-data/gpcp-monthly-global-precipitation-climatology-project 上找到。
空间趋势
数据读取
import numpy as np
import scipy.stats as stats
import datetime
from netCDF4 import Dataset as netcdf # netcdf4-python module
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
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
import xarray as xr
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6
import warnings
warnings.simplefilter('ignore')
path = r'I:/precip.mon.mean.nc'
# 读取数据
ds = xr.open_dataset(path).sel(time=slice('1979-01-01', '2023-12-01'))
lons = ds ['lon'][:]
lats = ds ['lat'][:]
time = ds ['time'][:]
pr = ds ['precip'][:]ds
趋势计算
计算趋势并分析
- 将降水数据 pr.data 重塑为 (nt, ngrd) 的形状,其中 nt 是时间步数,ngrd 是网格点数。order=‘F’ 表示按列优先重塑
- 初始化两个空数组 pr_rate 和 pr_val,用于存储每个网格点的降水变化率和p值,并将它们的值设为 NaN
- 遍历每个网格点 i
- 提取该网格点的时间序列数据 y,并去除 NaN 值,得到 y0
- 创建一个与 y0 长度相同的线性空间 x
- 使用 stats.linregress 进行线性回归,计算斜率 slope、截距 intercept、相关系数 r_value、p值 p_value 和标准误差 std_err
- 将斜率乘以120(假设数据是月度数据,转换为每十年的变化率)并存储在 pr_rate 中
- 将p值存储在 pr_val 中
undef = -99999.0
pr.data[pr.data==undef] = np.nan
nt,nlat,nlon = pr.shape
ngrd = nlat*nlon
nt,nlat,nlon
pr_grd = pr.data.reshape((nt, ngrd), order='F')
pr_rate = np.empty((ngrd,1))
pr_rate[:,:] = np.nan
pr_val = np.empty((ngrd,1))
pr_val[:,:] = np.nanfor i in range(ngrd): print(i)y = pr_grd[:,i] y0 = y[~np.isnan(y)] x = np.linspace(1,len(y0), len(y0))slope, intercept, r_value, p_value, std_err = stats.linregress(x, y0)pr_rate[i,0] = slope*120.0 pr_val[i,0] = p_value pr_rate = pr_rate.reshape((nlat,nlon), order='F')
pr_val = pr_val.reshape((nlat,nlon), order='F')
使用basemap可视化空间趋势分布
plt.figure(figsize=(12,6),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.6, 0.6, 25)
cs = m.contourf(x, y, pr_rate.squeeze(), clevs, cmap='RdYlBu_r')
m.drawcoastlines()
m.drawparallels(np.arange(-90.,91.,30.), labels=[1,0,0,0], fontsize=10)
m.drawmeridians(np.arange(-180.,181.,60.), labels=[0,0,0,1], fontsize=10)
cb = m.colorbar(cs)
plt.title('GPCP Changing Rate (mm/day/decade)', fontsize=16)sig_lats, sig_lons = np.where(pr_val < 0.05)
m.scatter(lons[sig_lons], lats[sig_lats], marker='o', color='k', s=1)
plt.show()
气候态平均
pr_ann_clm = np.nanmean(pr, axis=0)
plt.figure(figsize=(12,6),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.0, 12.0, 25)
cs = m.contourf(x, y, pr_ann_clm.squeeze(), clevs, cmap=plt.cm.RdBu_r)
m.drawcoastlines()
m.drawparallels(np.arange(-90.,91.,30.), labels=[1,0,0,0], fontsize=10)
m.drawmeridians(np.arange(-180.,181.,60.), labels=[0,0,0,1], fontsize=10)
cb = m.colorbar(cs)
plt.title('GPCP v2.2 Annual Mean (mm/day)', fontsize=16)
全球平均 | 纬向平均 | 经向平均
- 计算每个网格点的权重,权重是纬度的余弦值,用于加权平均。纬度越高,权重越小。
- 遍历每个时间步 it。
- 提取该时间步的降水数据 pone。
- 将数据中的 NaN 值掩盖,创建一个掩码数组 pone。
- 使用权重计算该时间步的加权平均降水量,并存储在 pr_glb_avg 中。
lonx, latx = np.meshgrid(lons, lats)
weights = np.cos(latx * np.pi / 180.)pr_glb_avg = np.zeros(nt)for it in np.arange(nt):pone = pr[it,:, :]pone = np.ma.masked_array(pone, mask=np.isnan(pone))pr_glb_avg[it] = np.average(pone, weights=weights) glb_avg = np.mean(pr_glb_avg)
print(glb_avg)
print(pr_glb_avg.shape)fig, ax = plt.subplots(1, 1, figsize=(15,6))ax.plot(time, pr_glb_avg, color='b', linewidth=2)ax.axhline(glb_avg, linewidth=2, color='r', label="Global areal mean: " + str(np.round(glb_avg*100)/100))
ax.legend()
ax.set_title('GPCP areal mean (1979-2023)', fontsize=16)
ax.set_xlabel('Month/Year #', fontsize=12)
ax.set_ylabel('Preicipitation [mm/day]', fontsize=12)
ax.set_ylim(2.4, 2.9)
ax.minorticks_on()
fig.autofmt_xdate()
ax.fmt_xdata = mdates.DateFormatter('%Y')
纬向平均
pr_ann_clm = np.nanmean(pr, axis=0)
pr_ann_zonal = np.nanmean(pr_ann_clm, axis=1)
pr_ann_clm.shape
fig, ax = plt.subplots(1, 1, figsize=(15,6),dpi=200)ax.plot(lats, pr_ann_zonal, color='b', linewidth=2)ax.axhline(glb_avg, linewidth=2, color='r', label="Global areal mean: " + str(np.round(glb_avg*100)/100))
ax.legend()
ax.set_title('GPCP zonal mean (1979-2023)', fontsize=16)
ax.set_xlabel('Month/Year #', fontsize=12)
ax.set_ylabel('Preicipitation [mm/day]', fontsize=12)
ax.set_ylim(0.0, 6.0)
ax.set_xlim(-90, 90)
ax.xaxis.set_major_formatter(LatitudeFormatter())
ax.minorticks_on()
经向平均
pr_ann_mer = np.nanmean(pr_ann_clm, axis=0)
fig, ax = plt.subplots(1, 1, figsize=(15,6),dpi=200)
ax.plot(lons, pr_ann_mer, color='b', linewidth=2)
ax.axhline(glb_avg, linewidth=2, color='r', label="Global areal mean: " + str(np.round(glb_avg*100)/100))
ax.legend()
ax.set_title('GPCP Meridional mean (1979-2023)', fontsize=16)
ax.set_xlabel('Month/Year #', fontsize=12)
ax.set_ylabel('Preicipitation [mm/day]', fontsize=12)
ax.set_ylim(1, 3.5)
ax.set_xlim(0, 361)
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.minorticks_on()
总结
巩固了趋势分析、显著性打点、计算了纬向加权平均的全球分布、纬向分布、经向分布,并使用basemap绘制了相关空间分布图以及时间序列。
相关数据和line by line的代码ipynb文件放到了GitHub上:
- https://github.com/Blissful-Jasper/jianpu_record