https://imgur.com/a/FADxycQ
如上图所示
档案内座标系统的经度是0-360
cartopy使用的座标系统经度应该是-180到180(我没看错的话)
我原先没有转换座标系统就出问题了
但转换座标系统后后还是有问题
请各位大神帮我看看这个应该怎么处理
程式如下
import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
from scipy.signal import detrend
import cartopy.crs as ccrs
# 读取.nc档案
dataset = xr.open_dataset('D:\PCCU\python\ERA5 monthly averaged data on single levels from 1940 to present.nc')
# 选择1979到2020年间的数据
dataset = dataset.sel(time=slice('1979', '2020'))
# 选择目标范围
precip_dataset = dataset.sel(longitude=slice(105, 125), latitude=slice(28, 20))
sst_dataset = dataset.sel(longitude=slice(80, 300), latitude=slice(60, -60))
# 将经度资料转换到 -180 到 180 的范围
longitude = sst_dataset.longitude.values
longitude[longitude > 180] = longitude[longitude > 180] - 360
sst_dataset = sst_dataset.assign_coords(longitude=longitude)
# 提取并转换降水资料单位
precip = precip_dataset['tp'].values * 1000 # 单位转换
# 提取SST资料
sst = sst_dataset['sst'].values
# 进行空间平均
precip_mean = precip.mean(axis=(1, 2))
# 降水资料去趋势
precip_detrended = detrend(precip_mean)
# 为绘制散布图和回归直线准备资料
time = np.arange(len(precip_detrended)) # 创建一个与 precip_detrended 相同长度的时间数组
# 执行回归并获取斜率和截距
slope, intercept, _, _, _ = stats.linregress(time, precip_detrended)
# 绘制散布图和回归直线
plt.figure(figsize=(10, 6))
plt.scatter(time, precip_detrended, label='Detrended Precipitation') # 绘制散布图
plt.plot(time, slope*time + intercept, color='red', label='Regression Line') # 绘制回归直线
plt.xlabel('Time')
plt.ylabel('Detrended Precipitation')
plt.legend()
plt.title('Scatter plot of Detrended Precipitation and Regression Line')
plt.show()
# 初始化一个空的阵列来存储斜率 (与 sst 没有时间维度的形状相同)
slopes = np.empty(sst.shape[1:])
# 对每个网格点进行回归
for i in range(sst.shape[1]):
for j in range(sst.shape[2]):
finite_mask = np.isfinite(sst[:, i, j]) & np.isfinite(precip_detrended)
if finite_mask.sum() > 0:
slope, _, _, _, _ = stats.linregress(sst[finite_mask, i, j], precip_detrended[finite_mask])
slopes[i, j] = slope
else:
slopes[i, j] = np.nan
# 绘图
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=180)) #设定图中央位在太平洋上
p = ax.pcolormesh(sst_dataset.longitude, sst_dataset.latitude, slopes, cmap='RdBu_r')
ax.coastlines()
ax.set_extent([80, -60, -60, 60]) # 将绘图范围设定为东经80度到西经60度，北纬60度到南纬60度
plt.colorbar(p, ax=ax, orientation='vertical', label='Slope')
plt.title('Regression slopes of SST and Detrended Precipitation')
plt.show()
谢谢