如何使用Python對NetCDF數據做空間相關分析

import numpy as np #數值計算用，如相關系數
import xarray as xr #讀取.nc文件用
from sklearn.feature_selection import f_regression #做顯著性檢驗
import matplotlib.pyplot as plt #繪制和展示圖形用
import cartopy.crs as ccrs #繪制地圖用，如果沒有安裝好的話，請在後臺聯系我
import cartopy.feature as cfeature #添加一些矢量用，這裡沒用到，因為我沒數據
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter #經緯度格式設置
import cmaps #ncl的color，如果沒有的話，請聯系我，也可以在氣象傢園找到

#使用上下文管理器讀取.nc數據，並提取數據中的變量，可以提前用NASA的panoply這個軟件查看.nc信息
with xr.open_dataset(r'D:\inuyasha\codeX\codeLEARN\sst.DJF.mean.anom.nc') as f1:
      pre = f1['sst_anom'][:-1, :, :]  # 三維數據全取，時間，緯度+經度
      lat, lon = f1['lat'], f1['lon'] #提取經緯度，後面格網化需要用到
pre2d = np.array(pre).reshape(pre.shape[0], pre.shape[1]*pre.shape[2])
#0表示行個數，1列代表的個數，2經度代表個數
with xr.open_dataset(r'D:\inuyasha\codeX\codeLEARN\pc.DJF.sst.nc') as f2:
      pc = f2['pc'][0, :]

# 相關系數計算
pre_cor = np.corrcoef(pre2d.T, pc)[:-1, -1].reshape(len(lat), len(lon))

# 做顯著性檢驗
pre_cor_sig = f_regression(np.nan_to_num(pre2d), pc)[1].reshape(len(lat), len(lon))#用0代替NaN
area = np.where(pre_cor_sig < 0.05)
# numpy的作用又來瞭 
nx, ny = np.meshgrid(lon, lat)  
# 格網化經緯度，打印出來看看就知道為什麼要這麼做瞭
plt.figure(figsize=(16, 8)) #創建一個空畫佈
#讓colorbar字體設置為新羅馬字符
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['font.size'] = 16

ax2 = plt.subplot(projection=ccrs.PlateCarree(central_longitude=180))
# 在畫佈上繪圖，這個叫axes，這不是坐標軸喔
ax2.coastlines(lw=0.4)
ax2.set_global()
c2 = ax2.contourf(nx, ny, pre_cor, extend='both', cmap=cmaps.nrl_sirkes, transform=ccrs.PlateCarree())
plt.colorbar(c2,fraction=0.05,orientation='horizontal', shrink=0.4, pad=0.06)
# extend關鍵字設置colorbar的形狀，both為兩端尖的，pad是距離主圖的距離，其他參數web搜索

# 顯著性打點
sig2 = ax2.scatter(nx[area], ny[area], marker='+', s=1, c='k', alpha=0.6, transform=ccrs.PlateCarree())
# 凸顯顯著性區域
plt.title('Correlation Analysis', fontdict={'family' : 'Times New Roman', 'size'   : 16})
#標題字體也修改為新羅馬字符，數字和因為建議都用新羅馬字符
ax2.set_xticks(np.arange(0, 361, 30),crs=ccrs.PlateCarree())
# 經度范圍設置，nunpy的作用這不就又來瞭嘛
plt.xticks(fontproperties = 'Times New Roman',size=16) #修改xy刻度字體為新羅馬字符
plt.yticks(fontproperties = 'Times New Roman',size=16)
ax2.set_yticks(np.arange(-90, 90, 15),crs=ccrs.PlateCarree())
# 設置y
ax2.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label = False))#經度0度不加東西
ax2.yaxis.set_major_formatter(LatitudeFormatter())
# 設置經緯度格式，就是多少度顯示那樣的，而不是一些數字
ax2.set_extent([-178, 178, -70, 70], crs=ccrs.PlateCarree())
# 設置空間范圍
plt.grid(color='k')
# 畫一個網格吧
plt.show()
# 顯示出圖形

import numpy as np
import xarray as xr
from sklearn.feature_selection import f_regression
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cmaps

with xr.open_dataset(r'D:\inuyasha\codeX\codeLEARN\sst.DJF.mean.anom.nc') as f1:
      pre = f1['sst_anom'][:-1, :, :]  # 三維數據全取，時間，緯度+經度
      lat, lon = f1['lat'], f1['lon']
pre2d = np.array(pre).reshape(pre.shape[0], pre.shape[1]*pre.shape[2])#0行代表的個數，1緯度，2經度
#pre2d.shape是一個39行，16020列的矩陣，T之後就變為瞭16020行，39列

with xr.open_dataset(r'D:\inuyasha\codeX\codeLEARN\pc.DJF.sst.nc') as f2:
      pc = f2['pc'][0, :]
#pc是一個39行的數組

# # 相關系數
pre_cor = np.corrcoef(pre2d.T, pc)[:-1, -1].reshape(len(lat), len(lon))
#pre_cor.shape,(16020,)->reshape(89,180)
# # 顯著性檢驗

# pre_cor_sig = f_regression(np.nan_to_num(pre2d), pc)[1].reshape(len(lat), len(lon))#用0代替NaN
# area = np.where(pre_cor_sig < 0.05)

nx, ny = np.meshgrid(lon, lat)  # 格網化
nx,ny