引言：我一直想理解空間相關(guān)分析的計(jì)算思維，于是今天又拿起Python腳本和數(shù)據(jù)來做練習(xí)。首先需要說明的是，這次實(shí)驗(yàn)的數(shù)據(jù)和Python腳本均來自于[好久不見]大佬，在跟大佬說明之后，允許我寫到公眾號(hào)來與大家共享，在此對(duì)大佬的指點(diǎn)表示感謝，這次實(shí)驗(yàn)的腳本可在氣象家園或簡(jiǎn)書app(如果沒記錯(cuò)的話)搜索到這次實(shí)驗(yàn)的相關(guān)內(nèi)容，也可以微信或者后臺(tái)發(fā)消息給我獲取。在此之前我覺得自己還沒理解這個(gè)方法的計(jì)算思維，檢驗(yàn)的標(biāo)準(zhǔn)就是我能否迅速運(yùn)用到其他方面。于是今天又重新回來溫習(xí)一遍，我把自己的理解與大伙共同交流。
首先，數(shù)據(jù)的格式是NetCDF(.nc)數(shù)據(jù)，兩個(gè)數(shù)據(jù)分別是[哈德來中心海溫sst數(shù)據(jù)，pc數(shù)據(jù)是對(duì)東太平洋SSTA做的EOF獲取]。知道數(shù)據(jù)信息之后我們就準(zhǔn)備開始去運(yùn)行程序。原始腳本包括了回歸分析和相關(guān)分析兩部分，但是今天我做了空間相關(guān)分析這一部分，有興趣的可以到[好久不見]大佬的氣象家園閱讀喔！如果還沒有安裝Cartopy包的話請(qǐng)?jiān)诤笈_(tái)聯(lián)系我喔
為了方便理解每一步，我選擇去Jupyter運(yùn)行，因?yàn)榭梢砸欢我欢纬绦虻倪\(yùn)行，這是比較方便的。繪圖部分并不是很難，關(guān)鍵還是在于數(shù)據(jù)預(yù)處理部分。
空間相關(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代替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)建一個(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è)置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')
# 畫一個(gè)網(wǎng)格吧
plt.show()
# 顯示出圖形
那么就運(yùn)行看看效果吧
如果覺得這個(gè)color不喜歡的話，就換一下ncl的來吧，ncl的顏色多而漂亮，喜歡啥就換啥
想要理解這個(gè)方法的計(jì)算思維，有必要觀察原始數(shù)據(jù)和數(shù)據(jù)處理之后的樣式，理解了數(shù)據(jù)樣式之后可能更有助于我們理解整個(gè)程序
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ù)全取，時(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
看看格網(wǎng)化后的經(jīng)緯度多規(guī)范啊。畫張圖來看看可能也會(huì)直觀一些。
好吧，今天的分享就到這里了，理解了這個(gè)計(jì)算思維，能更好地遷移運(yùn)用到其他研究方面，如果還沒有安裝Cartopy包的話請(qǐng)?jiān)诤笈_(tái)聯(lián)系我喔，如果需要測(cè)試數(shù)據(jù)和腳本請(qǐng)?jiān)诤笈_(tái)聯(lián)系我，當(dāng)然也可以去[好久不見]大佬的主頁(yè)。如果覺得這次分享不錯(cuò)的話，還請(qǐng)老鐵們點(diǎn)個(gè)贊，多多分享，歡迎交流學(xué)習(xí)，感謝各位！
原始資料：
http://bbs.06climate.com/forum.php?mod=viewthreadtid=92816highlight=%CF%D4%D6%F8%D0%D4%BC%EC%D1%E9%2B%CF%E0%B9%D8%B7%D6%CE%F6
以上就是如何使用Python對(duì)NetCDF數(shù)據(jù)做空間相關(guān)分析的詳細(xì)內(nèi)容，更多關(guān)于Python對(duì)NetCDF數(shù)據(jù)做空間分析的資料請(qǐng)關(guān)注腳本之家其它相關(guān)文章！
您可能感興趣的文章:- python 微信好友特征數(shù)據(jù)分析及可視化
- python3對(duì)拉勾數(shù)據(jù)進(jìn)行可視化分析的方法詳解
- Python數(shù)據(jù)分析:手把手教你用Pandas生成可視化圖表的教程
- Python數(shù)據(jù)可視化正態(tài)分布簡(jiǎn)單分析及實(shí)現(xiàn)代碼
- python數(shù)據(jù)分析之公交IC卡刷卡分析
- python數(shù)據(jù)分析之用sklearn預(yù)測(cè)糖尿病
- Python數(shù)據(jù)分析之pandas函數(shù)詳解
- python爬蟲之你好,李煥英電影票房數(shù)據(jù)分析
- python基于scrapy爬取京東筆記本電腦數(shù)據(jù)并進(jìn)行簡(jiǎn)單處理和分析
- python數(shù)據(jù)分析之員工個(gè)人信息可視化
