雖然近些年來氣象數(shù)據(jù)基本都以netcdf (.nc)格式為主隧熙，但是不同機(jī)構(gòu)婉弹，不同模式生成的數(shù)據(jù)仍存在差異睬魂，比如說有的經(jīng)度坐標(biāo)是0°-360°，有的是-180°-180°镀赌，有的緯度坐標(biāo)是從南到北氯哮，有的緯度坐標(biāo)是從北到南，然而這些處理起來都還算很簡單的商佛。唯獨(dú)時(shí)間模塊喉钢，真的是讓人麻爪。尤其是在python中良姆，時(shí)間格式的數(shù)據(jù)有可能是string肠虽，有可能是datetime64，有可能是object等等玛追，每種都需要用不同的方法去調(diào)整税课，就算是同一種格式，也可能在處理的過程中遇到各種各樣的問題痊剖。那么我就把自己親身經(jīng)歷過的各種崩潰瞬間分享一下韩玩，提供的解決辦法不一定是最優(yōu)解，僅供參考陆馁，大家有更好的辦法也可以評論留言找颓。

一、多模式結(jié)果處理中的時(shí)間格式不統(tǒng)一

這是之前在處理CMIP6模式中遇到的問題叮贩，我一共是處理了24個(gè)CMIP6模式(SSP245击狮，以及historical兩個(gè)情景)，里邊就有那么幾個(gè)模式獨(dú)立特行妇汗，跟別人格格不入帘不。大部分的模式很乖，很好處理杨箭，時(shí)間格式默認(rèn)就是datetime64格式的寞焙，比如以ACCESS-CM2模式的tas變量為例(僅利用CDO通過mergetime和線性插值預(yù)處理了原始文件)：

import xarray as xr
f = xr.open_dataset('./ACCESS-CM2/tas/t2m.nc')
print(f)

文件信息如下：

<pre style="box-sizing: border-box; overflow: auto; font-family: monospace; font-size: 14px; display: block; padding: 1px 0px; margin: 0px; line-height: inherit; color: rgb(0, 0, 0); word-break: break-all; overflow-wrap: break-word; background-color: rgb(255, 255, 255); border: 0px; border-radius: 0px; white-space: pre-wrap; vertical-align: baseline; font-style: normal; font-variant-ligatures: normal; font-variant-caps: normal; font-weight: 400; letter-spacing: normal; orphans: 2; text-align: left; text-indent: 0px; text-transform: none; widows: 2; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration-style: initial; text-decoration-color: initial;"><xarray.Dataset>
Dimensions:    (lat: 73, lon: 144, nb2: 2, time: 31390)
Coordinates:
  * lon        (lon) float64 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * lat        (lat) float64 -90.0 -87.5 -85.0 -82.5 ... 82.5 85.0 87.5 90.0
  * time       (time) datetime64[ns] 2015-01-01T12:00:00 ... 2100-12-31T12:00:00
Dimensions without coordinates: nb2
Data variables:
    time_bnds  (time, nb2) datetime64[ns] ...
    tas        (time, lat, lon) float32 ...

那我們可以看到，其時(shí)間格式為datetime64[ns]互婿，這是最方便的格式捣郊，可以直接索引。比如說我需要選取2015年到2099年冬季的數(shù)據(jù)慈参，那么直接通過以下指令就可以實(shí)現(xiàn)：

tas = f['tas'].loc[f.time.dt.month.isin([12,1,2])].loc['2015-12-01':'2100-02-28']
print(tas)

輸出信息如下：

<xarray.DataArray 'tas' (time: 7650, lat: 73, lon: 144)>
[80416800 values with dtype=float32]
Coordinates:
  * lon      (lon) float64 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * lat      (lat) float64 -90.0 -87.5 -85.0 -82.5 -80.0 ... 82.5 85.0 87.5 90.0
  * time     (time) datetime64[ns] 2015-12-01T12:00:00 ... 2100-02-28T12:00:00
Attributes:
    standard_name:  air_temperature
    long_name:      Near-Surface Air Temperature
    units:          K
    comment:        near-surface (usually, 2 meter) air temperature
    cell_methods:   area: time: mean
    history:        2019-11-08T10:42:10Z altered by CMOR: Treated scalar dime...

那么呛牲，當(dāng)遇到另一種時(shí)間格式時(shí)，就不能直接這樣索引了驮配，比如說接下來的這個(gè)：

f = xr.open_dataset('/data/home/zenggang/yxy/CMIP6/SSP245/CanESM5/tas/t2m.nc')    
print(f)

其文件信息如下：

<xarray.Dataset>
Dimensions:    (lat: 73, lon: 144, nb2: 2, time: 31390)
Coordinates:
  * lon        (lon) float64 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * lat        (lat) float64 -90.0 -87.5 -85.0 -82.5 ... 82.5 85.0 87.5 90.0
  * time       (time) object 2015-01-01 12:00:00 ... 2100-12-31 12:00:00
Dimensions without coordinates: nb2
Data variables:
    time_bnds  (time, nb2) object ...
    tas        (time, lat, lon) float32 ...

注意咯娘扩，這里是object着茸，那么我們還以之前的方法索引會發(fā)生什么呢？

tas = f['tas'].loc[f.time.dt.month.isin([12,1,2])].loc['2015-12-01':'2100-02-28']

image.png

哎琐旁？格式不支持涮阔，好氣哦！
不過Xarray庫提供了一個(gè)函數(shù)來進(jìn)行轉(zhuǎn)換灰殴，因此這個(gè)格式也還算友好敬特。

import xarray as xr
f = xr.open_dataset('/data/home/zenggang/yxy/CMIP6/SSP245/CanESM5/tas/t2m.nc')    
f= f.assign_coords(time = f.indexes['time'].to_datetimeindex())
tas = f['tas'].loc[f.time.dt.month.isin([12,1,2])].loc['2015-12-01':'2100-02-28']

后來，我又遇到了另一種情況牺陶，就是有的模式的時(shí)間是0時(shí)(比如說 1979-01-01 00:00:00)伟阔，而有的模式的時(shí)間是12時(shí)(比如說 1979-01-01 12:00:00)，如果需要統(tǒng)一的話掰伸，只需要搭配datetime庫進(jìn)行調(diào)整即可皱炉，比如說將所有的12時(shí)變?yōu)?時(shí)：

from datetime import datetime, timedelta
f = f.assign_coords(time = (f.indexes['time'].to_datetimeindex()- timedelta(hours=12)))        

利用assign_coords重新聲明坐標(biāo)時(shí)，將所有的時(shí)間減去了12小時(shí)碱工。(利用assign_coords也可以重新聲明經(jīng)緯度坐標(biāo)等等)

還有一類模式有點(diǎn)反人類(不點(diǎn)名了)娃承，我目前也沒想到太好的辦法來處理，它所使用的不是標(biāo)準(zhǔn)日歷怕篷，他每個(gè)月都有30天，2月份也有30天酗昼，很離譜廊谓，這在xarray庫中無法被識別(轉(zhuǎn)換)，簡單粗暴點(diǎn)就直接用CDO刪掉了每年的這一天麻削，或者只能將全年的數(shù)據(jù)索引出來蒸痹，通過reshape成(年，月呛哟，日叠荠，經(jīng)度，緯度)的方法通過下標(biāo)序號索引扫责。

二榛鼎、非標(biāo)準(zhǔn)日歷數(shù)據(jù)

這個(gè)是我處理CAM模式數(shù)據(jù)時(shí)遇到的問題(run了30年的敏感性試驗(yàn)，生成的數(shù)據(jù)的時(shí)間是從0001-01-01到0032-01-01)鳖孤，后面處理時(shí)讓我瞬間崩潰了者娱，生成的數(shù)據(jù)是cftime庫中的時(shí)間格式，但是表面還是寫著object類型苏揣，通過前邊介紹的assign_coords的方法根本行不通黄鳍，報(bào)錯(cuò)理由是

image.png

深層的原因就不在這里解釋了，只提供一個(gè)解決的方法平匈，思路還是將時(shí)間轉(zhuǎn)為datetime64格式：

import xarray as xr
import pandas as pd
f = xr.open_dataset('./TS.nc')
time_tmp1 = f.indexes['time']
time_tmp2 = []
j = 0
for i in range(time_tmp1.shape[0]):
    if ((i+1)%365==0):
        j = j+1
    else:
        j = j 
    tmp = time_tmp1[i].replace(year = 2000+j)
    a = tmp.year
    b = tmp.month
    c = tmp.day
    time_tmp2.append(pd.to_datetime('{}/{}/{}'.format(a,b,c),format='%Y/%m/%d')) 
f = f.assign_coords(time = time_tmp2)

中間那一大段的循環(huán)的目的框沟，是將每個(gè)時(shí)間坐標(biāo)統(tǒng)統(tǒng)加上2000年藏古，使其從0001-01-01到0030-12-31變?yōu)?000-01-01到2032-01-01，這樣就又可以通過正常的時(shí)間索引方法進(jìn)行索引了忍燥。
該文件默認(rèn)格式的時(shí)間是cftime拧晕，其存在一個(gè)方法，.replace可以替換年灾前，月防症，日，然后通過.year, .month,.day獲取新的年月日哎甲，再通過pd.to_datetime函數(shù)變?yōu)閐atetime64格式的日期蔫敲，最后重新賦值給f即可。
原始f信息：

<xarray.Dataset>
Dimensions:    (lat: 96, lon: 144, nb2: 2, time: 11681)
Coordinates:
  * lon        (lon) float64 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * lat        (lat) float64 -90.0 -88.11 -86.21 -84.32 ... 86.21 88.11 90.0
  * time       (time) object 0001-01-01 00:00:00 ... 0033-01-01 00:00:00
Dimensions without coordinates: nb2
Data variables:
    time_bnds  (time, nb2) object ...
    TREFHT     (time, lat, lon) float32 ...
    TS         (time, lat, lon) float32 ...
    TSMN       (time, lat, lon) float32 ...
    TSMX       (time, lat, lon) float32 ...

修改后f信息：

<xarray.Dataset>
Dimensions:    (lat: 96, lon: 144, nb2: 2, time: 11681)
Coordinates:
  * lon        (lon) float64 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * lat        (lat) float64 -90.0 -88.11 -86.21 -84.32 ... 86.21 88.11 90.0
  * time       (time) datetime64[ns] 2000-01-01 2000-01-02 ... 2032-01-01
Dimensions without coordinates: nb2
Data variables:
    time_bnds  (time, nb2) object ...
    TREFHT     (time, lat, lon) float32 ...
    TS         (time, lat, lon) float32 ...
    TSMN       (time, lat, lon) float32 ...
    TSMX       (time, lat, lon) float32 ...

三炭玫、批量索引特定時(shí)間的數(shù)據(jù)

我主要的方向是極端天氣氣候事件奈嘿，因此嘗嘗需要批量索引特定時(shí)間的數(shù)據(jù)，比如說先統(tǒng)計(jì)40年里發(fā)生的極端(高溫吞加，低溫)事件的日期裙犹，然后從40年的逐日位勢高度數(shù)據(jù)里提取出發(fā)生事件當(dāng)天的位勢高度異常，有沒有不寫一大長串的判斷循環(huán)就能實(shí)現(xiàn)的方法呢衔憨？

方法1：

在統(tǒng)計(jì)事件時(shí)叶圃，順便統(tǒng)計(jì)好每次事件發(fā)生時(shí)的時(shí)間坐標(biāo)索引序號(是這40年里的第幾天發(fā)生的)
然后通過：

z = np.array([z_all[i]  for i in number])

實(shí)現(xiàn)索引，其中践图，z為取出的發(fā)生極端事件當(dāng)天的位勢高度掺冠，z_all為這40年的逐日位勢高度，number是每次事件發(fā)生時(shí)的時(shí)間坐標(biāo)索引序號码党，比如說事件是這40年里的第3,6,8...天發(fā)生的德崭，那么number就是[2,5,7,...]，這里z和z_all都是np.array
實(shí)際上這里的本質(zhì)還是通過循環(huán)解決的揖盘。

方法2：通過xarray直接索引特定時(shí)間

首先在統(tǒng)計(jì)事件時(shí)眉厨，統(tǒng)計(jì)好發(fā)生每次事件的時(shí)間，datetime64格式的兽狭！憾股，索引時(shí)通過

z = z_all.loc[date]

直接完成。
其中z和z_all是xr.DataArray格式椭符，date是統(tǒng)計(jì)好的時(shí)間列表荔燎。
并且，對于分析事件個(gè)例來說销钝，通常還會提取爆發(fā)前后的一段時(shí)間進(jìn)行逐日的分析有咨，那么通過如下方法實(shí)現(xiàn)：

z_1day_before = z_all.loc[date - np.timedelta64(1,'D')]

這里以發(fā)生前一天的數(shù)據(jù)提取為例。D表示day

小結(jié)

看到這里大家應(yīng)該可以發(fā)現(xiàn)蒸健，能否用python處理好數(shù)據(jù)的關(guān)鍵在于編程者能否將不同的庫的方法和函數(shù)靈活的組合起來座享，單純依靠基礎(chǔ)語言或者單一的庫婉商，很難實(shí)現(xiàn)一些靈活的操作，也無法體現(xiàn)出py的方便和強(qiáng)大之處渣叛。有時(shí)候我看到曾經(jīng)寫過的一些代碼丈秩，都想不起來當(dāng)時(shí)是怎么能想到這樣處理的，因?yàn)殪`活就意味著多種解決方法淳衙。實(shí)在不行的話蘑秽，遇事不決就簡單粗暴的使用循環(huán)判斷吧！