GMT绘时间序列及相分析的方法

小二

GMT绘时间序列及相分析

三目标：时间序列，[url=]加趋势[/url][h1] ，频谱分析。本节为其一：时间序列的绘制。

   使用matlab也可以很方便的绘制时间序列，但是对于时间轴的设置却比较困难。例如minor ticks的设置和标签位置剧中要求的设置。因此使用GMT绘制时间序列。主要是为了方便次级月标签的绘制。

第一节：时间序列

1绘制时间序列，手动调整坐标刻度，一年为一个主刻度（main ticks），一月为一个次刻度（minor ticks），并且绘制Y坐标的grid。

参考样式：如下2图（选自中国知网，以及内部资料）所示.

我们依照上面的时间序列样式，使用GMT同样绘制类似图形。

最终的结果将下图：

具体步骤如下：

1首先是准备好规则数据。如果数据格式是matlab处理和存储的矩阵格式，使用MATLAB则需要提取出来时间序列，并保存为ASCII格式数据。注意此时的数据不带有时间信息，只是简单的一列数据而已，我们将在在后面添加时间信息。

2为sigma0数据添加时间列[h2]

这个方法不是固定的，可以自己编程实现，也可以使用已有的简单程序。注意时间数据的格式是固定的几种。

如果数据不多的话，可以使用Windows自带的excel添加。

本例子的最终数据格式如下：第一列是时间，第二列到第四列是sigma0的属性值，（time，sigma0Ku，sigma0C，sigma0C-Ku）。

3 GMT绘图

将三列sigma0的数据放到一幅图上，即同时绘制三条时间序列，且使用同一个坐标系。横坐标为年月，具体设置为：年为大刻度，月为小刻度，界限设置为1993.1.1-2004.12.1。竖坐标为sigma0的值，范围可当已知，设置为0-15。

本例子参考GMT Cookbook中的程序21：

GMT中Time格式需要按照数据的时间格式设置。在GMT的帮助文件中可以找到，有格里高利，ISO之类。还要区分输入的时间格式和输出的时间格式。

我们选择例21中使用的时间标准dd-o-yy。（GMT有细致的时间设置参数，注意在这一步别出错误。）

注意：数据文件中时间序列的日期坐标可以和标注的日期坐标格式不一样，但是时间上要一致。

如本例数据格式：

rem 主要代码，

rem []内容参考文后注释。

gmtset[url=]INPUT_DATE_FORMAT dd-o-yy[/url][刘晓辉3] [url=]PLOT_DATE_FORMAT[/url][刘晓辉4] o ANNOT_FONT_SIZE_PRIMARY +10p

gmtset[url=]TIME_FORMAT_PRIMARY[/url][刘晓辉5] abbreviated CHAR_ENCODING ISOLatin1+

坐标轴设置，前面的是时间起点和终点，后面的0/300是竖坐标的界限。

[url=]set[/url][刘晓辉6] R=-R1999-08-11T00:00:00/2006-11-01T00:00:00/0/300

gmtconvert为数据提取程序，可以提取多列数据文件中的某一列或者几列数据。

[url=]gmtconvert[/url][刘晓辉7] -F0,2 -f0T -Hi RHAT_price.csv > RHAT.env

psxy为绘制曲线的程序，可以设置点和线的样式。

[url=]psxy[/url][刘晓辉8] -R -J -Gyellow -O -K RHAT.env >> example_21.ps

4最终的GMT bat文件

代码如下：

gmtset INPUT_DATE_FORMAT dd-o-yy PLOT_DATE_FORMAT o ANNOT_FONT_SIZE_PRIMARY +10p

gmtset TIME_FORMAT_PRIMARY abbreviated CHAR_ENCODING ISOLatin1+

REM Create a suitable -R string

set R=-R1993-01-01T00:00:00/2004-12-31T00:00:00/0/15

REM Lay down the basemap:

psbasemap %R% -JX9i/4i -Gwhite -K -Bs1Y:"t/a":/:,db:WSen -Bpf1o/2WSen:."Time series of sigma0 over Xinjiang China": > example_21.ps

gmtset OUTPUT_DATE_FORMAT dd-o-yy

gmtconvert -F0,2 -f0T XJsigma0_TS.txt > RHAT.env

gmtconvert -F0,3 -f0T XJsigma0_TS.txt > RHAT2.env

gmtconvert -F0,1 -f0T XJsigma0_TS.txt > RHAT1.env

rem red C blue Ku black C-Ku

psxy -R -J -O -K -Wthin -Sc0.1c -Gblack RHAT2.env >> example_21.ps

psxy -R -J -O -K -Wthinner RHAT2.env >> example_21.ps

psxy -R -J -O -K -Wthin,blue RHAT.env >> example_21.ps

psxy -R -J -O -K -St0.2c -Gblue RHAT.env >> example_21.ps

psxy -R -J RHAT1.env -Wthin,red -O -K >> example_21.ps

psxy -R -J -O -K -Si0.2c -Gred RHAT1.env >> example_21.ps

rem gmtset INPUT_DATE_FORMAT dd-o-yy

REM Clean up after ourselves:

del RHAT.*

del .gmt*

5时间序列的结果如下：

至此，时间序列完成了。

PS：网格虚线，也可以自动化实现。修改default中的grid primary size，在粗细后面加上“，-”。然后在psbasemap中的-B属性中加g控制网格密度。

下一步是使用GMT分析其频率域以及趋势分析。(建议自己尝试)

2添加趋势。如下图所示。
趋势不必在所有的时间序列中添加，目前只在明显的添加。但是要给出趋势线的主要参数。并且在时间序列图中显示。

使用MATLAB可以做到，但是不能与第一步的图结合。

因此先全部用GMT作图。

trend 1d

3 FFT绘图。

116个数据FFT的MATLAB程序还没调好。这个也尝试用GMT绘制。

spectrum 1d

grdfft?

Matlab的字符串数组不知道怎么设置。

用excel却可以很好的完成时间的添加。法为：先写出第一个时间起点，然后设置单元格样式，设置为

然后托动鼠标下拉即可。

然后再修改时间格式为dd-mmm-yy.

使用matlab编程难以实现，却使用excel实现。

输入格式

输出格式In addition, you may use a singleoinstead of mm (to plot month name)

You may specify either Gregorian calendar format or ISO week calendar format.

As examples, try yy/mm/dd, yyyy=jjj, dd-o-yyyy, dd-mm-yy,yy-mm, etc.

Controls how primary month-, week-, and weekday-names are formatted. Choose amongfull, abbreviated, and character.If the leading f, a, or c are replaced with F, A, and C the entire annotation will be in upper case.

设置标注使用的时间界限。

后面的/0/300则是Y轴的界限。

提取第0，2列，-f0T表示时间