如何使用Python對(duì)NetCDF數(shù)據(jù)做空間相關(guān)分析

import numpy as np #數(shù)值計(jì)算用，如相關(guān)系數(shù)import xarray as xr #讀取.nc文件用from sklearn.feature_selection import f_regression #做顯著性檢驗(yàn)import matplotlib.pyplot as plt #繪制和展示圖形用import cartopy.crs as ccrs #繪制地圖用，如果沒有安裝好的話，請(qǐng)?jiān)诤笈_(tái)聯(lián)系我import cartopy.feature as cfeature #添加一些矢量用，這里沒用到，因?yàn)槲覜]數(shù)據(jù)from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter #經(jīng)緯度格式設(shè)置import cmaps #ncl的color，如果沒有的話，請(qǐng)聯(lián)系我，也可以在氣象家園找到
#使用上下文管理器讀取.nc數(shù)據(jù)，并提取數(shù)據(jù)中的變量，可以提前用NASA的panoply這個(gè)軟件查看.nc信息with xr.open_dataset(r'D:\inuyasha\codeX\codeLEARN\sst.DJF.mean.anom.nc') as f1:      pre = f1['sst_anom'][:-1, :, :]  # 三維數(shù)據(jù)全取，時(shí)間，緯度+經(jīng)度      lat, lon = f1['lat'], f1['lon'] #提取經(jīng)緯度，后面格網(wǎng)化需要用到pre2d = np.array(pre).reshape(pre.shape[0], pre.shape[1]*pre.shape[2])#0表示行個(gè)數(shù)，1列代表的個(gè)數(shù)，2經(jīng)度代表個(gè)數(shù)with xr.open_dataset(r'D:\inuyasha\codeX\codeLEARN\pc.DJF.sst.nc') as f2:      pc = f2['pc'][0, :]
# 相關(guān)系數(shù)計(jì)算pre_cor = np.corrcoef(pre2d.T, pc)[:-1, -1].reshape(len(lat), len(lon))
# 做顯著性檢驗(yàn)pre_cor_sig = f_regression(np.nan_to_num(pre2d), pc)[1].reshape(len(lat), len(lon))#用0代替NaNarea = np.where(pre_cor_sig < 0.05)# numpy的作用又來了 nx, ny = np.meshgrid(lon, lat)  # 格網(wǎng)化經(jīng)緯度，打印出來看看就知道為什么要這么做了plt.figure(figsize=(16, 8)) #創(chuàng)建一個(gè)空畫布#讓colorbar字體設(shè)置為新羅馬字符plt.rcParams['font.family'] = 'Times New Roman'plt.rcParams['font.size'] = 16
ax2 = plt.subplot(projection=ccrs.PlateCarree(central_longitude=180))# 在畫布上繪圖，這個(gè)叫axes，這不是坐標(biāo)軸喔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關(guān)鍵字設(shè)置colorbar的形狀，both為兩端尖的，pad是距離主圖的距離，其他參數(shù)web搜索
# 顯著性打點(diǎn)sig2 = ax2.scatter(nx[area], ny[area], marker='+', s=1, c='k', alpha=0.6, transform=ccrs.PlateCarree())# 凸顯顯著性區(qū)域plt.title('Correlation Analysis', fontdict={'family' : 'Times New Roman', 'size'   : 16})#標(biāo)題字體也修改為新羅馬字符，數(shù)字和因?yàn)榻ㄗh都用新羅馬字符ax2.set_xticks(np.arange(0, 361, 30),crs=ccrs.PlateCarree())# 經(jīng)度范圍設(shè)置，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())# 設(shè)置yax2.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label = False))#經(jīng)度0度不加?xùn)|西ax2.yaxis.set_major_formatter(LatitudeFormatter())# 設(shè)置經(jīng)緯度格式，就是多少度顯示那樣的，而不是一些數(shù)字ax2.set_extent([-178, 178, -70, 70], crs=ccrs.PlateCarree())# 設(shè)置空間范圍plt.grid(color='k')# 畫一個(gè)網(wǎng)格吧plt.show()# 顯示出圖形

import numpy as npimport xarray as xrfrom sklearn.feature_selection import f_regressionimport matplotlib.pyplot as pltimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatterimport cmaps
with xr.open_dataset(r'D:\inuyasha\codeX\codeLEARN\sst.DJF.mean.anom.nc') as f1:      pre = f1['sst_anom'][:-1, :, :]  # 三維數(shù)據(jù)全取，時(shí)間，緯度+經(jīng)度      lat, lon = f1['lat'], f1['lon']pre2d = np.array(pre).reshape(pre.shape[0], pre.shape[1]*pre.shape[2])#0行代表的個(gè)數(shù)，1緯度，2經(jīng)度#pre2d.shape是一個(gè)39行，16020列的矩陣，T之后就變?yōu)榱?6020行，39列
with xr.open_dataset(r'D:\inuyasha\codeX\codeLEARN\pc.DJF.sst.nc') as f2:      pc = f2['pc'][0, :]#pc是一個(gè)39行的數(shù)組
# # 相關(guān)系數(shù)pre_cor = np.corrcoef(pre2d.T, pc)[:-1, -1].reshape(len(lat), len(lon))#pre_cor.shape,(16020,)->reshape(89,180)# # 顯著性檢驗(yàn)
# 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)  # 格網(wǎng)化nx,ny