--------------------?數(shù)據(jù)處理專題一 --------------------?

;說(shuō)在前面：在處理數(shù)據(jù)之前，你首先需要明確你的目的-你的數(shù)據(jù)是最終是用來(lái)干嘛的。第二需要了解你的數(shù)據(jù)是啥樣的，針對(duì)不同類型的數(shù)據(jù)有不同的處理方式-因材施教。

;;; 使用目的：氣候態(tài)？趨勢(shì)？趨勢(shì)是季節(jié)、年際或是年代際。相關(guān)？等等

;;; 原始數(shù)據(jù)：一般的氣象數(shù)據(jù)格式以.nc居多，但也有g(shù)rib,tiff,csv,txt等等。

;; 如何對(duì)原始數(shù)據(jù)有一個(gè)初步的了解？3個(gè)方法：

; 1、一般情況下，數(shù)據(jù)網(wǎng)站會(huì)有描述，請(qǐng)仔細(xì)閱讀。

; 2、利用可視化軟件panoply或者linux端安裝cdo或者nc_filedump(only for .NC)。

; 3、利用ncl-printVarsummary

;;; 第一步讀取數(shù)據(jù),下面主要以nc數(shù)據(jù)為例 橙色字體為NCL自帶函數(shù)，可在官網(wǎng)搜索

f = addfile("precip.nc","r")

precip = f->precip ;傻瓜式讀取

; 與之相反的是有意識(shí)地讀取。那么問(wèn)題來(lái)了什么是有意識(shí)地讀取呢？這里就體現(xiàn)提前理解原始數(shù)據(jù)的重要性了

; 各數(shù)據(jù)網(wǎng)站，數(shù)據(jù)存儲(chǔ)方式五花八門(mén)

; 但NCL最喜歡的方式是(time,level,lat,lon),所有維度按順序排列，這里對(duì)于各個(gè)維度代表什么不作過(guò)多解釋

; 因此，非常有必要將個(gè)性張揚(yáng)的數(shù)據(jù)進(jìn)行變動(dòng)一番

; eg. GLEAM(time,lat,lon)外表看著挺正常一數(shù)據(jù)，printVarsummary一下，你會(huì)發(fā)現(xiàn)

printVarsummary(precip)

;; 它是這樣的

; Variable: precip

; Type: float

; Total Size: 378432000 values

; ? ? ? ? ? ? 1513728000 bytes

; Number of Dimensions: 3

; Dimensions and sizes: ? ? ? [ 365 <time | unlimited> x 720 <lat> x 1440 <lon> ]

; Chunking Info: ? ? ?[ 1 <time> x 45 <lat> x 1440 <lon> ]

; Coordinates:

; ? ? ? ? ? ? time: [ 0..364]

; ? ? ? ? ? ? lat: [89.875..-89.875]

; ? ? ? ? ? ? lon: [-179.875..179.875]

; ? ? Number of Attributes: ? ? ? ?5

; ? ? ? ? standard_name ? ? ? : ? ? ? Precipitation

; ? ? ? ? long_name ? : ? ? ? Precipitation from GLEAM

; ? ? ? ? units ? ? ? : ? ? ? mm/day

; ? ? ? ? _FillValue ?: ? ? ? -999

; ? ? ? ? missing_value ? ? ? : ? ? ? -999

;; 第二維 lat 逆序了，此時(shí)需要轉(zhuǎn)個(gè)序。其實(shí)如果你不處理也沒(méi)關(guān)系，但不保證后面進(jìn)行計(jì)算和plot時(shí)不出問(wèn)題。

precip = f->precip(:,::-1,:) ;有意識(shí)地讀取

printVarsummary(precip)

; Variable: precip

; Type: float

; Total Size: 378432000 values

; ? ? ? ? ? ? 1513728000 bytes

; Number of Dimensions: 3

; Dimensions and sizes: ? ? ? [ 365 <time | unlimited> x 720 <lat> x 1440 <lon> ]

; Chunking Info: ? ? ?[ 1 <time> x 45 <lat> x 1440 <lon> ]

; Coordinates:

; ? ? ? ? ? ? time: [ 0..364]

; ? ? ? ? ? ? lat: [-89.875..89.875]

; ? ? ? ? ? ? lon: [-179.875..179.875]

; ? ? Number of Attributes: ? ? ? ?5

; ? ? ? ? standard_name ? ? ? : ? ? ? Precipitation

; ? ? ? ? long_name ? : ? ? ? Precipitation from GLEAM

; ? ? ? ? units ? ? ? : ? ? ? mm/day

; ? ? ? ? _FillValue ?: ? ? ? -999

; ? ? ? ? missing_value ? ? ? : ? ? ? -999

;; 以上是其一

; 其二 假如該數(shù)據(jù) 時(shí)間尺度 19820101-19821231 而你只需要其中一段數(shù)據(jù)，eg. 19820403-19820422，此時(shí)你需要借助別的函數(shù)進(jìn)行處理

time = f->time

yyyymmddd = cd_calendar(time, -2);gregorian轉(zhuǎn)換成標(biāo)準(zhǔn)時(shí) 年月日 格式, 第二維參數(shù)可以修改，根據(jù)數(shù)據(jù)或者需求進(jìn)行修改，詳情見(jiàn)NCL官網(wǎng)

ymd_ind = ind(yyyymmddd.ge.19820402.and.yyyymmddd.le.19820422) ;ind為索引函數(shù)，這一行的目的是通過(guò)ind函數(shù)截取上述時(shí)間范圍內(nèi)的數(shù)據(jù)。

; 重新讀取

precip = f->precip(ymd_ind,::-1,:) ;有意識(shí)地讀取

printVarsummary(precip)

; 除此之外，還可以按研究區(qū)范圍進(jìn)行讀取 eg.(10-20°N，-20-20°E)

precip = f->precip(ymd_ind,{10:20:-1},{-20:20}) ;有意識(shí)地讀取

printVarsummary(precip)

temp = f->t2m(ymd_ind,{10:20:-1},{-20:20})

;; 氣候態(tài) 有兩種途徑：第一通過(guò)循環(huán)；第二利用NCL自帶函數(shù)

; 這里只介紹用函數(shù) 多一事不如少一事，畢竟NCL跑起循環(huán)來(lái)效率特別低

; 年氣候態(tài)

;1、保證數(shù)據(jù)時(shí)間尺度為年尺度 下面以日尺度數(shù)據(jù)為例

;a 用cdo將日換成年 cdo -yearmean -cat 'files*nc' yearlyMeanOutput.nc

temp_clmyear = dim_avg_n_Wrap(temp_year, 0) ;

; 月氣候態(tài)

;b 用cdo將日換成月 cdo monavg in.nc out.nc

temp_clmmon = clmMonTLL(temp_mon, yyyyddd) ;

;季節(jié)氣候態(tài),夏季

temp_jja = month_to_season(temp_mon, "JJA")

; 日氣候態(tài)

; a 創(chuàng)建所需的 yyyyddd

TIME = cd_calendar(time, 0) ; 類型為 float

year = toint(TIME(:,0)) ; toint 去除元數(shù)據(jù)

month = toint(TIME(:,1))

day = toint(TIME(:,2))

; b 檢查日歷屬性

if (isatt(TIME,"calendar")) then ; 默認(rèn)為公歷

? ? year@calendar = TIME@calendar

end if

ddd = day_of_year(year, month, day)

if (isatt(year,"calendar")) then ; 默認(rèn)為公歷

? ? ddd@calendar = year@calendar

end if

yyyyddd = year*1000 + ddd ; 輸入所需

if (isatt(ddd,"calendar")) then ; 默認(rèn)為公歷

? ? yyyyddd@calendar = ddd@calendar

end if

; 計(jì)算日氣候態(tài)：原始日均值

pre_clmday = clmDayTLL(pre, yyyyddd) ;每個(gè)網(wǎng)格點(diǎn)的日氣候?qū)W

; ?。?！ 利用日氣候態(tài)計(jì)算月氣候態(tài)

monthly_climatology = new((/12,pre_clmday&lat,pre_clmday&lon/), typeof(pre_clmday))

do month = 0, 11

? ? days_in_month = days_in_month(month + 1)

? ? start_day = sum(days_in_month(0:month-1))

? ? end_day = start_day + days_in_month(month) - 1

? ? monthly_climatology(month,:,:) = dim_avg_n_Wrap(pre_clmday(start_day:end_day,:,:), 0)

end do

; ?。。±萌諝夂驊B(tài)計(jì)算每個(gè)季節(jié)的平均季節(jié)氣候態(tài)

seasonal_climatology = new((/4,pre_clmday&lat,pre_clmday&lon/), typeof(pre_clmday))

do season = 0, 3

? ? start_month = season * 3

? ? end_month = start_month + 2

? ? start_day = sum(days_in_month(0:start_month-1))

? ? end_day = start_day + sum(days_in_month(start_month:end_month)) - 1

? ? seasonal_climatology(season,:,:) =?dim_avg_n_Wrap(pre_clmday(start_day:end_day,:,:), 0)

end do

; ?。。?！利用日氣候態(tài)計(jì)算年氣候態(tài)

annual_climatology = dim_sum_n_Wrap(pre_clmday, 0)

;;;; 方法多變，具體使用什么，如何使用請(qǐng)結(jié)合需求進(jìn)行選擇。以上代碼僅供參考，歡迎糾錯(cuò)。

;;------------------------------------ Ending ---------------------------------------

寫(xiě)在后面：一般情況下，私信我看到都會(huì)回復(fù)，但可能不及時(shí)，如果遇到?jīng)]回復(fù)的，請(qǐng)先想一下你是不是一上來(lái)就問(wèn)代碼在哪？脾氣不好，對(duì)于這種情況我一般會(huì)視而不見(jiàn)。所以愿意看您就看，不愿看請(qǐng)您退出。

鏈接：https://pan.baidu.com/s/1YwGu9IBkttFNDkD1lFAuvw?pwd=NCL0?

提取碼：NCL0?