關(guān)于python處理氣象數(shù)據(jù)過程中常遇到的關(guān)于時間處理的問題

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

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

這是之前在處理CMIP6模式中遇到的問題，我一共是處理了24個CMIP6模式(SSP245，以及historical兩個情景)，里邊就有那么幾個模式獨立特行，跟別人格格不入。大部分的模式很乖，很好處理，時間格式默認就是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 ...

那我們可以看到，其時間格式為datetime64[ns]，這是最方便的格式，可以直接索引。比如說我需要選取2015年到2099年冬季的數(shù)據(jù)，那么直接通過以下指令就可以實現(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)遇到另一種時間格式時，就不能直接這樣索引了，比如說接下來的這個：

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庫提供了一個函數(shù)來進行轉(zhuǎn)換，因此這個格式也還算友好。

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']

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

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

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

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

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

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


image.png

深層的原因就不在這里解釋了，只提供一個解決的方法，思路還是將時間轉(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)的目的，是將每個時間坐標(biāo)統(tǒng)統(tǒng)加上2000年，使其從0001-01-01到0030-12-31變?yōu)?000-01-01到2032-01-01，這樣就又可以通過正常的時間索引方法進行索引了。
該文件默認格式的時間是cftime，其存在一個方法，.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ù)據(jù)

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

方法1：

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

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

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

方法2：通過xarray直接索引特定時間

首先在統(tǒng)計事件時，統(tǒng)計好發(fā)生每次事件的時間，datetime64格式的！，索引時通過

z = z_all.loc[date]

直接完成。
其中z和z_all是xr.DataArray格式，date是統(tǒng)計好的時間列表。
并且，對于分析事件個例來說，通常還會提取爆發(fā)前后的一段時間進行逐日的分析，那么通過如下方法實現(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ǔ)語言或者單一的庫，很難實現(xiàn)一些靈活的操作，也無法體現(xiàn)出py的方便和強大之處。有時候我看到曾經(jīng)寫過的一些代碼，都想不起來當(dāng)時是怎么能想到這樣處理的，因為靈活就意味著多種解決方法。實在不行的話，遇事不決就簡單粗暴的使用循環(huán)判斷吧！

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容