您好,登錄后才能下訂單哦!
這篇文章主要介紹了怎么使用Python對NetCDF數(shù)據(jù)做空間相關(guān)分析,具有一定借鑒價值,感興趣的朋友可以參考下,希望大家閱讀完這篇文章之后大有收獲,下面讓小編帶著大家一起了解一下。
python常用的庫:1.requesuts;2.scrapy;3.pillow;4.twisted;5.numpy;6.matplotlib;7.pygama;8.ipyhton等。
引言:我一直想理解空間相關(guān)分析的計算思維,于是今天又拿起Python腳本和數(shù)據(jù)來做練習(xí)。首先需要說明的是,這次實驗的數(shù)據(jù)和Python腳本均來自于[好久不見]大佬,在跟大佬說明之后,允許我寫到公眾號來與大家共享,在此對大佬的指點表示感謝,這次實驗的腳本可在氣象家園或簡書app(如果沒記錯的話)搜索到這次實驗的相關(guān)內(nèi)容,也可以微信或者后臺發(fā)消息給我獲取。在此之前我覺得自己還沒理解這個方法的計算思維,檢驗的標(biāo)準(zhǔn)就是我能否迅速運用到其他方面。于是今天又重新回來溫習(xí)一遍,我把自己的理解與大伙共同交流。
首先,數(shù)據(jù)的格式是NetCDF(.nc)數(shù)據(jù),兩個數(shù)據(jù)分別是[哈德來中心海溫sst數(shù)據(jù),pc數(shù)據(jù)是對東太平洋SSTA做的EOF獲取]。知道數(shù)據(jù)信息之后我們就準(zhǔn)備開始去運行程序。原始腳本包括了回歸分析和相關(guān)分析兩部分,但是今天我做了空間相關(guān)分析這一部分,有興趣的可以到[好久不見]大佬的氣象家園閱讀喔!如果還沒有安裝Cartopy包的話請在后臺聯(lián)系我喔
為了方便理解每一步,我選擇去Jupyter運行,因為可以一段一段程序的運行,這是比較方便的。繪圖部分并不是很難,關(guān)鍵還是在于數(shù)據(jù)預(yù)處理部分。
空間相關(guān)分析的腳本如下:
import numpy as np #數(shù)值計算用,如相關(guān)系數(shù) import xarray as xr #讀取.nc文件用 from sklearn.feature_selection import f_regression #做顯著性檢驗 import matplotlib.pyplot as plt #繪制和展示圖形用 import cartopy.crs as ccrs #繪制地圖用,如果沒有安裝好的話,請在后臺聯(lián)系我 import cartopy.feature as cfeature #添加一些矢量用,這里沒用到,因為我沒數(shù)據(jù) from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter #經(jīng)緯度格式設(shè)置 import cmaps #ncl的color,如果沒有的話,請聯(lián)系我,也可以在氣象家園找到 #使用上下文管理器讀取.nc數(shù)據(jù),并提取數(shù)據(jù)中的變量,可以提前用NASA的panoply這個軟件查看.nc信息 with xr.open_dataset(r'D:\inuyasha\codeX\codeLEARN\sst.DJF.mean.anom.nc') as f1: pre = f1['sst_anom'][:-1, :, :] # 三維數(shù)據(jù)全取,時間,緯度+經(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表示行個數(shù),1列代表的個數(shù),2經(jīng)度代表個數(shù) with xr.open_dataset(r'D:\inuyasha\codeX\codeLEARN\pc.DJF.sst.nc') as f2: pc = f2['pc'][0, :] # 相關(guān)系數(shù)計算 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) # 格網(wǎng)化經(jīng)緯度,打印出來看看就知道為什么要這么做了 plt.figure(figsize=(16, 8)) #創(chuàng)建一個空畫布 #讓colorbar字體設(shè)置為新羅馬字符 plt.rcParams['font.family'] = 'Times New Roman' plt.rcParams['font.size'] = 16 ax2 = plt.subplot(projection=ccrs.PlateCarree(central_longitude=180)) # 在畫布上繪圖,這個叫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搜索 # 顯著性打點 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ù)字和因為建議都用新羅馬字符 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è)置y ax2.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') # 畫一個網(wǎng)格吧 plt.show() # 顯示出圖形
那么就運行看看效果吧
如果覺得這個color不喜歡的話,就換一下ncl的來吧,ncl的顏色多而漂亮,喜歡啥就換啥
想要理解這個方法的計算思維,有必要觀察原始數(shù)據(jù)和數(shù)據(jù)處理之后的樣式,理解了數(shù)據(jù)樣式之后可能更有助于我們理解整個程序
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, :, :] # 三維數(shù)據(jù)全取,時間,緯度+經(jīng)度 lat, lon = f1['lat'], f1['lon'] pre2d = np.array(pre).reshape(pre.shape[0], pre.shape[1]*pre.shape[2])#0行代表的個數(shù),1緯度,2經(jīng)度 #pre2d.shape是一個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是一個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) # # 顯著性檢驗 # 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
看看格網(wǎng)化后的經(jīng)緯度多規(guī)范啊。畫張圖來看看可能也會直觀一些。
感謝你能夠認(rèn)真閱讀完這篇文章,希望小編分享的“怎么使用Python對NetCDF數(shù)據(jù)做空間相關(guān)分析”這篇文章對大家有幫助,同時也希望大家多多支持億速云,關(guān)注億速云行業(yè)資訊頻道,更多相關(guān)知識等著你來學(xué)習(xí)!
免責(zé)聲明:本站發(fā)布的內(nèi)容(圖片、視頻和文字)以原創(chuàng)、轉(zhuǎn)載和分享為主,文章觀點不代表本網(wǎng)站立場,如果涉及侵權(quán)請聯(lián)系站長郵箱:is@yisu.com進(jìn)行舉報,并提供相關(guān)證據(jù),一經(jīng)查實,將立刻刪除涉嫌侵權(quán)內(nèi)容。