溫馨提示×
您好，登錄后才能下訂單哦！
點擊 登錄注冊 即表示同意《億速云用戶服務條款》
您好，登錄后才能下訂單哦！
這篇文章主要介紹“如何使用Python對NetCDF數據做空間相關分析”，在日常操作中，相信很多人在如何使用Python對NetCDF數據做空間相關分析問題上存在疑惑，小編查閱了各式資料，整理出簡單好用的操作方法，希望對大家解答”如何使用Python對NetCDF數據做空間相關分析”的疑惑有所幫助！接下來，請跟著小編一起來學習吧！
引言：我一直想理解空間相關分析的計算思維，于是今天又拿起Python腳本和數據來做練習。首先需要說明的是，這次實驗的數據和Python腳本均來自于[好久不見]大佬，在跟大佬說明之后，允許我寫到公眾號來與大家共享，在此對大佬的指點表示感謝，這次實驗的腳本可在氣象家園或簡書app(如果沒記錯的話)搜索到這次實驗的相關內容，也可以微信或者后臺發消息給我獲取。在此之前我覺得自己還沒理解這個方法的計算思維，檢驗的標準就是我能否迅速運用到其他方面。于是今天又重新回來溫習一遍，我把自己的理解與大伙共同交流。
首先，數據的格式是NetCDF(.nc)數據，兩個數據分別是[哈德來中心海溫sst數據，pc數據是對東太平洋SSTA做的EOF獲取]。知道數據信息之后我們就準備開始去運行程序。原始腳本包括了回歸分析和相關分析兩部分，但是今天我做了空間相關分析這一部分，有興趣的可以到[好久不見]大佬的氣象家園閱讀喔！如果還沒有安裝Cartopy包的話請在后臺聯系我喔
為了方便理解每一步，我選擇去Jupyter運行，因為可以一段一段程序的運行，這是比較方便的。繪圖部分并不是很難，關鍵還是在于數據預處理部分。
空間相關分析的腳本如下：
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代替NaNarea = 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'] = 16ax2 = 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())# 設置yax2.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()# 顯示出圖形
那么就運行看看效果吧
如果覺得這個color不喜歡的話，就換一下ncl的來吧，ncl的顏色多而漂亮，喜歡啥就換啥
想要理解這個方法的計算思維，有必要觀察原始數據和數據處理之后的樣式，理解了數據樣式之后可能更有助于我們理解整個程序
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 cmapswith 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
看看格網化后的經緯度多規范啊。畫張圖來看看可能也會直觀一些。
到此，關于“如何使用Python對NetCDF數據做空間相關分析”的學習就結束了，希望能夠解決大家的疑惑。理論與實踐的搭配能更好的幫助大家學習，快去試試吧！若想繼續學習更多相關知識，請繼續關注億速云網站，小編會繼續努力為大家帶來更多實用的文章！
免責聲明：本站發布的內容（圖片、視頻和文字）以原創、轉載和分享為主，文章觀點不代表本網站立場，如果涉及侵權請聯系站長郵箱：is@yisu.com進行舉報，并提供相關證據，一經查實，將立刻刪除涉嫌侵權內容。
