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()
謝謝