关于python处理气象数据过程中常遇到的关于时间处理的问题

虽然近些年来气象数据基本都以netcdf (.nc)格式为主,但是不同机构,不同模式生成的数据仍存在差异,比如说有的经度坐标是0°-360°,有的是-180°-180°,有的纬度坐标是从南到北,有的纬度坐标是从北到南,然而这些处理起来都还算很简单的。唯独时间模块,真的是让人麻爪。尤其是在python中,时间格式的数据有可能是string,有可能是datetime64,有可能是object等等,每种都需要用不同的方法去调整,就算是同一种格式,也可能在处理的过程中遇到各种各样的问题。那么我就把自己亲身经历过的各种崩溃瞬间分享一下,提供的解决办法不一定是最优解,仅供参考,大家有更好的办法也可以评论留言。

一、多模式结果处理中的时间格式不统一

这是之前在处理CMIP6模式中遇到的问题,我一共是处理了24个CMIP6模式(SSP245,以及historical两个情景),里边就有那么几个模式独立特行,跟别人格格不入。大部分的模式很乖,很好处理,时间格式默认就是datetime64格式的,比如以ACCESS-CM2模式的tas变量为例(仅利用CDO通过mergetime和线性插值预处理了原始文件):

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年冬季的数据,那么直接通过以下指令就可以实现:

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...

那么,当遇到另一种时间格式时,就不能直接这样索引了,比如说接下来的这个:

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,那么我们还以之前的方法索引会发生什么呢?

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

哎?格式不支持,好气哦!
不过Xarray库提供了一个函数来进行转换,因此这个格式也还算友好。

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),如果需要统一的话,只需要搭配datetime库进行调整即可,比如说将所有的12时变为0时:

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

利用assign_coords重新声明坐标时,将所有的时间减去了12小时。(利用assign_coords也可以重新声明经纬度坐标等等)

还有一类模式有点反人类(不点名了),我目前也没想到太好的办法来处理,它所使用的不是标准日历,他每个月都有30天,2月份也有30天,很离谱,这在xarray库中无法被识别(转换),简单粗暴点就直接用CDO删掉了每年的这一天,或者只能将全年的数据索引出来,通过reshape成(年,月,日,经度,纬度)的方法通过下标序号索引。

二、非标准日历数据

这个是我处理CAM模式数据时遇到的问题(run了30年的敏感性试验,生成的数据的时间是从0001-01-01到0032-01-01),后面处理时让我瞬间崩溃了,生成的数据是cftime库中的时间格式,但是表面还是写着object类型,通过前边介绍的assign_coords的方法根本行不通,报错理由是


image.png

深层的原因就不在这里解释了,只提供一个解决的方法,思路还是将时间转为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)

中间那一大段的循环的目的,是将每个时间坐标统统加上2000年,使其从0001-01-01到0030-12-31变为2000-01-01到2032-01-01,这样就又可以通过正常的时间索引方法进行索引了。
该文件默认格式的时间是cftime,其存在一个方法,.replace可以替换年,月,日,然后通过.year, .month,.day获取新的年月日,再通过pd.to_datetime函数变为datetime64格式的日期,最后重新赋值给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 ...

三、批量索引特定时间的数据

我主要的方向是极端天气气候事件,因此尝尝需要批量索引特定时间的数据,比如说先统计40年里发生的极端(高温,低温)事件的日期,然后从40年的逐日位势高度数据里提取出发生事件当天的位势高度异常,有没有不写一大长串的判断循环就能实现的方法呢?

方法1:

在统计事件时,顺便统计好每次事件发生时的时间坐标索引序号(是这40年里的第几天发生的)
然后通过:

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

实现索引,其中,z为取出的发生极端事件当天的位势高度,z_all为这40年的逐日位势高度,number是每次事件发生时的时间坐标索引序号,比如说事件是这40年里的第3,6,8...天发生的,那么number就是[2,5,7,...],这里z和z_all都是np.array
实际上这里的本质还是通过循环解决的。

方法2:通过xarray直接索引特定时间

首先在统计事件时,统计好发生每次事件的时间,datetime64格式的!,索引时通过

z = z_all.loc[date]

直接完成。
其中z和z_all是xr.DataArray格式,date是统计好的时间列表。
并且,对于分析事件个例来说,通常还会提取爆发前后的一段时间进行逐日的分析,那么通过如下方法实现:

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

这里以发生前一天的数据提取为例。D表示day

小结

看到这里大家应该可以发现,能否用python处理好数据的关键在于编程者能否将不同的库的方法和函数灵活的组合起来,单纯依靠基础语言或者单一的库,很难实现一些灵活的操作,也无法体现出py的方便和强大之处。有时候我看到曾经写过的一些代码,都想不起来当时是怎么能想到这样处理的,因为灵活就意味着多种解决方法。实在不行的话,遇事不决就简单粗暴的使用循环判断吧!

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,294评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,493评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,790评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,595评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,718评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,906评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,053评论 3 410
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,797评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,250评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,570评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,711评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,388评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,018评论 3 316
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,796评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,023评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,461评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,595评论 2 350

推荐阅读更多精彩内容