TimeSeriesAnalysis

时间序列分析

这个例子的作者是Nikolay Koldunov，节取自他的系列文章Python for Geosciences notes.原文请参见 https://github.com/koldunovn/python_for_geosciences/blob/master/06%20-%20Time%20series%20analysis%20(Pandas).ipynb

本文介绍了利用pandas包进行时间序列分析的例子，作者指出，事实上一维数据都可以利用pandas高效地解决

注：Docklet上pandas的版本号和原作者不同，部分实例输出结果有细节上的差别，这里我们以Docklet上的运行结果为准，但不影响数据处理的结果

引用需要的模块

首先我们引用必要的模块

In [1]:
import pandas as pd
import numpy as np
pd.set_option('max_rows',15) # 这里限制了行数的最大值

打开笔记本内嵌的显示图形功能：

In [2]:
%matplotlib inline

pandas发展快速，即便我们只使用基本的功能，一些实现上的细节在新版本中仍可能有变化

In [4]:
pd.__version__
Out[4]:
'0.15.0'

加载数据

现在我们做些准备工作，获取数据。可以直接使用wget命令：

In [5]:
!wget http://www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ao_index/monthly.ao.index.b50.current.ascii

--2016-04-24 01:14:51--  http://www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ao_index/monthly.ao.index.b50.current.ascii
Resolving www.cpc.ncep.noaa.gov (www.cpc.ncep.noaa.gov)... 140.90.101.63
Connecting to www.cpc.ncep.noaa.gov (www.cpc.ncep.noaa.gov)|140.90.101.63|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 19875 (19K) [text/plain]
Saving to: ‘monthly.ao.index.b50.current.ascii’

monthly.ao.index.b5 100%[===================>]  19.41K  44.2KB/s    in 0.4s    

2016-04-24 01:14:52 (44.2 KB/s) - ‘monthly.ao.index.b50.current.ascii’ saved [19875/19875]

实际上pandas拥有强大的I/O能力，但这个例子中我们并不涉及。这里简单地用numpy loadtxt打开文件：

In [6]:
ao = np.loadtxt('monthly.ao.index.b50.current.ascii')

文件中的每一行包含三项：年、月、值：

In [7]:
ao[0:2]
Out[7]:
array([[  1.95000000e+03,   1.00000000e+00,  -6.03100000e-02],
       [  1.95000000e+03,   2.00000000e+00,   6.26810000e-01]])

可以查看数组的形式(shape):

In [8]:
ao.shape
Out[8]:
(795, 3)

时间序列

我们希望自然而简单地将这些数据转换成时间序列。作为第一步，我们选定我们的时间序列的时间范围。从文件可以很清晰地看出时间起点是2013年1月（此时作者开始写作本系列文章），你可以根据自己的文件调整截止时间。数据的频率（两次数据的时间间隔）是月。

In [9]:
dates = pd.date_range('1950-01', '2014-01', freq='M')

如你所见，Pandas语法简单，这正是作者喜欢它的原因之一。另外请注意，截止时间参数设成了2014-01而不是2013-12，因为区间右侧是开的

In [10]:
dates
Out[10]:
<class 'pandas.tseries.index.DatetimeIndex'>
[1950-01-31, ..., 2013-12-31]
Length: 768, Freq: M, Timezone: None
In [11]:
dates.shape
Out[11]:
(768,)

现在我们开始创建我们的第一个事件序列。变量dates所标明的日期将作为索引，AO值则作为我们的值。下面我们就用截止到2013年底的数据进行分析。

In [12]:
AO = pd.Series(ao[:768,2], index=dates)
AO
Out[12]:
1950-01-31   -0.060310
1950-02-28    0.626810
1950-03-31   -0.008127
1950-04-30    0.555100
1950-05-31    0.071577
...
2013-07-31   -0.011112
2013-08-31    0.154250
2013-09-30   -0.460880
2013-10-31    0.262760
2013-11-30    2.029000
2013-12-31    1.474900
Freq: M, Length: 768

我们可以绘制出完整的时间序列：

In [13]:
AO.plot()
Out[13]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f80d25a2e48>

也可以绘制出其中的一部分：

In [14]:
AO['1980':'1990'].plot()
Out[14]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f80d1ebfd30>

甚至更小的区间：

In [15]:
AO['1980-05':'1981-03'].plot()
Out[15]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f80d1e9bcc0>

参照时间以一种很自然的方式完成了。你也可以找某一个时间点，可以通过序号：

In [16]:
AO[120]
Out[16]:
-2.4842

或者通过索引（我们选定了date所标明的日期）：

In [17]:
AO['1960-01']
Out[17]:
1960-01-31   -2.4842
Freq: M, dtype: float64

也可以尝试查看一年的数据：

In [18]:
AO['1960']
Out[18]:
1960-01-31   -2.484200
1960-02-29   -2.212400
1960-03-31   -1.624600
1960-04-30   -0.297310
1960-05-31   -0.857430
1960-06-30    0.054978
1960-07-31   -0.619060
1960-08-31   -1.007900
1960-09-30   -0.381640
1960-10-31   -1.187000
1960-11-30   -0.553230
1960-12-31   -0.342950
Freq: M, dtype: float64

Data Frame

下面我们做一些更有趣的尝试。首先下载更多数据：

In [19]:
!wget http://www.cpc.ncep.noaa.gov/products/precip/CWlink/pna/norm.nao.monthly.b5001.current.ascii

--2016-04-24 01:22:15--  http://www.cpc.ncep.noaa.gov/products/precip/CWlink/pna/norm.nao.monthly.b5001.current.ascii
Resolving www.cpc.ncep.noaa.gov (www.cpc.ncep.noaa.gov)... 140.90.101.63
Connecting to www.cpc.ncep.noaa.gov (www.cpc.ncep.noaa.gov)|140.90.101.63|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 19080 (19K) [text/plain]
Saving to: ‘norm.nao.monthly.b5001.current.ascii’

norm.nao.monthly.b5 100%[===================>]  18.63K  42.4KB/s    in 0.4s    

2016-04-24 01:22:16 (42.4 KB/s) - ‘norm.nao.monthly.b5001.current.ascii’ saved [19080/19080]

创建nao序列，就像之前创建AO序列一样，时间区间也一样：

In [21]:
nao = np.loadtxt('norm.nao.monthly.b5001.current.ascii')
dates_nao = pd.date_range('1950-01', '2014-01', freq='M')
NAO = pd.Series(nao[:768,2], index=dates_nao)
NAO.index
Out[21]:
<class 'pandas.tseries.index.DatetimeIndex'>
[1950-01-31, ..., 2013-12-31]
Length: 768, Freq: M, Timezone: None

现在我们创建Data Frame，它将包含AO和NAO数据。可以直接绘制出数据：

In [22]:
aonao = pd.DataFrame({'AO' : AO, 'NAO' : NAO})
aonao.plot()
Out[22]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f80d1c899e8>

我们可以查看最初的几行：

In [23]:
aonao.head()
Out[23]:
AO NAO
1950-01-31 -0.060310 0.92
1950-02-28 0.626810 0.40
1950-03-31 -0.008127 -0.36
1950-04-30 0.555100 0.73
1950-05-31 0.071577 -0.59

我们可以通过列名查看每一列：

In [24]:
aonao['NAO']
Out[24]:
1950-01-31    0.92
1950-02-28    0.40
1950-03-31   -0.36
1950-04-30    0.73
1950-05-31   -0.59
...
2013-07-31    0.67216
2013-08-31    0.97019
2013-09-30    0.24060
2013-10-31   -1.28010
2013-11-30    0.90082
2013-12-31    0.94566
Freq: M, Name: NAO, Length: 768

或者通过DATA FRAME变量的方法查看：

In [25]:
aonao.NAO
Out[25]:
1950-01-31    0.92
1950-02-28    0.40
1950-03-31   -0.36
1950-04-30    0.73
1950-05-31   -0.59
...
2013-07-31    0.67216
2013-08-31    0.97019
2013-09-30    0.24060
2013-10-31   -1.28010
2013-11-30    0.90082
2013-12-31    0.94566
Freq: M, Name: NAO, Length: 768

可以很容易地添加列：

In [26]:
aonao['Diff'] = aonao['AO'] - aonao['NAO']
aonao.head()
Out[26]:
AO NAO Diff
1950-01-31 -0.060310 0.92 -0.980310
1950-02-28 0.626810 0.40 0.226810
1950-03-31 -0.008127 -0.36 0.351872
1950-04-30 0.555100 0.73 -0.174900
1950-05-31 0.071577 -0.59 0.661577

也可以很容易地删除列：

In [27]:
del aonao['Diff']
aonao.tail()
Out[27]:
AO NAO
2013-08-31 0.15425 0.97019
2013-09-30 -0.46088 0.24060
2013-10-31 0.26276 -1.28010
2013-11-30 2.02900 0.90082
2013-12-31 1.47490 0.94566

依然可以查看子区间：

In [28]:
aonao['1981-01':'1981-03']
Out[28]:
AO NAO
1981-01-31 -0.11634 0.37
1981-02-28 -0.33158 0.92
1981-03-31 -1.64470 -1.19

甚至以一组组合条件约束：

In [29]:
import datetime
aonao.ix[(aonao.AO > 0) & (aonao.NAO < 0) 
        & (aonao.index > datetime.datetime(1980,1,1)) 
        & (aonao.index < datetime.datetime(1989,1,1)),
        'NAO'].plot(kind='barh')
Out[29]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f80d259d668>

重采样

这里列出作者的最后一个综合性例子。Pandas提供了简单的重采样方法，两个主要参数分别是采样频率（“A”是一年，类似地，三年则是“3A”），采样方式，后者默认值是mean.

可以将多个重采样一并展示出来：

In [30]:
AO_mm = AO.resample("A", how=['mean', np.min, np.max])
AO_mm['1900':'2020'].plot()
Out[30]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f80d1c2aa20>
In [ ]:

results matching ""

    No results matching ""