Grads打点信度检验的绘制方法

游子

不想科研，就发张技术帖子吧。

做合成分析的同学，是不是很想在统计显著区域打点呢？

使用grads的传统方法是用灰度覆盖显著区域,但是很不美观。

如果我们希望使用shaded来显示数值，那么绘制显著区域只能使用等值线了。

如果我们希望同时使用contour和shaded来显示数值，而通过信度检验的区域打点，grads该如何

处理呢？

亲，这已经能够在2.0+版本中实现了。那么尝试一下呗。

先来看一下shapefile是个什么东西吧？

简言之，shapefile是一种数据集合，包含点、线、多边形等信息。通常由海岸线、国界、

州（省）界、气候区域，交通路线，河流等地理数据组成。每种空间特征都会被赋予几何学特征，以矢量的形式存储。一个shapefile文件通常包括四个部分：（1）主文件*.shp;(2)索引文件*.shx；（3）数据库文件*.dbf；(4)投影文件*.prj

我们可以在grads中创建shapefile并且将创建的文件用于打点绘图：

程序是合成分析，所以信度检验用到了t检验，挑选了10个异常年，所以样本可以看做10。

t检验的公式是t=(ave-u0)/sqrt(s)*sqrt(n)

其中n是样本量,ave是合成的均值，u0为总体均值，总体异常的均值是0，所以u0=0,公式退化为

t=ave/sqrt(s)*sqrt(n)

当n=10，我们查到自由度为n-1=9，达到95%信度的t值为2.26。所以欲达到信度t数值需要满足

abs(t)>2,26。程序中该部分的代码如下：

'define t1 = ave(tta,time=dec1853,time=feb1854)'

'define t2 = ave(tta,time=dec1855,time=feb1856)'

'define t3 = ave(tta,time=dec1862,time=feb1863)'

'define t4 = ave(tta,time=dec1888,time=feb1889)'

'define t5 = ave(tta,time=dec1933,time=feb1934)'

'define t6 = ave(tta,time=dec1954,time=feb1955)'

'define t7 = ave(tta,time=dec1960,time=feb1961)'

'define t8 = ave(tta,time=dec1970,time=feb1971)'

'define t9 = ave(tta,time=dec1979,time=feb1980)'

'define t10 = ave(tta,time=dec1988,time=feb1989)'

'define tave = (t1+t2+t3+t4+t5+t6+t7+t8+t9+t10)/10'

'define ss = pow(t1-tave,2)+pow(t2-tave,2)+pow(t3-tave,2)+pow(t4-tave,2)'

'define ss = pow(t5-tave,2)+pow(t6-tave,2)+pow(t7-tave,2)+pow(t8-tave,2)+ss'

'define ss = pow(t9-tave,2)+pow(t10-tave,2)+ss'

'ss = ss/10'

'define test = tave/sqrt(ss)*sqrt(10)'

挑选了10个个例，且计算了t的数值test

在grads绘图要用到maskout

'myp=maskout(test,test*test-2.26*2.26)'

*** or

'myp=maskout(test,abs(test)-2.26)'

'set gxout shp'

'set shp -pt shppt'

'd myp'

将通过信度的t值标记出来，并且存储成shapefile.

其中-pt可以理解为输出的数据将以点（PoinT）的形式予以存储。很关键喔。

在下一步将显著区域绘制出来即可。

grads越来越强大了。期待grads的2.1版本(支持多种字体了。。）

绘制出来的图的效果如下：

temperature T

由于这未完成的工作，数据就不给大家了。

实在想要数据的可以QQ联系

670143899

hujinggao

师弟，关于Shapefile还有很多用途，希望你统统研究一下和大家分享啊，顺便也教教我，呵呵。另外grads2.1不仅对字体，而且对出图的形式也作了很多修改，出图的质量更高，eps，ps等各种矢量图获取更方便。

平流层的萝卜

沙发，很给力！！

LOSER

楼主，您好，看了您算冬季的10个时次的平均，我如果我想计算30年平均，不会也要写t1，t2，…，t30……
因此想请教您，计算平均，还有其他方法吗

游子

LOSER 发表于 2013-11-6 08:32

                               
登录/注册后可看大图

楼主，您好，看了您算冬季的10个时次的平均，我如果我想计算30年平均，不会也要写t1，t2，…，t30……
因此 ...

如果是纯粹计算30年气候态，程序很简单。这个程序我是有的，循环即可以解决。

topmad

在R里面应该很容易就实现了  不过grads还是气象界的大头  

mofangbao

shapeFile用的很灵活呀，我也很期待新版的，不过开发人员老是跳票。。。

游子

hujinggao 发表于 2013-11-6 09:06

                               
登录/注册后可看大图

师弟，关于Shapefile还有很多用途，希望你统统研究一下和大家分享啊，顺便也教教我，呵呵。另外grads2.1不仅 ...

师兄，见笑了。

river

楼主很给力，支持一个！之前研究过一段时间，不过后来有事儿就耽搁了。再赞一下！

游子

river 发表于 2013-11-6 14:09

                               
登录/注册后可看大图

楼主很给力，支持一个！之前研究过一段时间，不过后来有事儿就耽搁了。再赞一下！

额，是么。我也是突然在文献中看到有人绘制grads打点图