IDL遥感应用入门(11):利用QuikScat数据绘制风速矢量图(色标表示大小箭头表示方向)

在之前的文章中我介绍了几种方法来绘制风场的矢量图,包括直接图形法和对象图形法,并且还做出了一系列的改进,具体的内容可以查看教程10教程9,这是两个经过改进修改的版本,算是比较完善。

现在我们换一种思路,表示二维矢量场的时候,我们可以这么来做,风速的标量值我们用色标来表示,而方向图我们用箭头来表示,这里的箭头将只表示方向,而不表示大小,也就是说我们把所有的箭头长短设为一致。

经过跟@杨铭伦 大神的讨论和我的一些尝试,这里给出一种方法可以实现上面所述的思路。下面是整个过程。

QuikScat数据说明

QuikScat提供3个方面数据信息:表面风速,表面风向和降雨信息。

其中数据在0到250属于有效信息,大于250的值均属于无效信息。

计算风速的时候需要乘以0.2,风向需要乘以1.5

更多详细信息可以查看这里,这个说明很重要,因为这决定了我们数据的读取方式。

效果预览

先来一张预览的效果图

所有代码

;Author:Sailor
;2013-11-15

PRO Course_10
;定义文件路径
MyRootDir='D:\3\'

;遍历文件夹
filearr = file_search(MyRootDir,'*.gz',count=num);
FOR fileindex=0,num-1,1 DO BEGIN

  ;打开文件
  OPENR , LUN , filearr[fileindex] , /COMPRESS,/GET_LUN

  ;构造显示标题
  filename=FILE_BASENAME(filearr[fileindex]);
  filename=STRMID(filename,6,10)
  year=STRMID(filename,0,4)+' '
  month=STRMID(filename,4,2)+' '
  day=STRMID(filename,6,2)+' '
  diplayname=year+month+day

  ;读入数组&抽样
  FleArr = BYTARR(1440,720,3)
  READU , LUN , FleArr

  DisSpd=FleArr(*,*,0)*0.2
  WdSpd = congrid(FleArr(*,*,0),90,45)
  WdDir = congrid(FleArr(*,*,1),90,45)

  ;过滤数据
  GT250IdxDir = where(WdDir GT 250);
  WdSpd[GT250IdxDir] = 0;

  GT250IdxDir = where(DisSpd GT 250*0.2);
  DisSpd[GT250IdxDir] = 0;
  DisSpd[GT250IdxDir] = max(DisSpd)+0.1;

  ;计算U和V,X和y
  WdSpd = WdSpd*0.2
  WdDir = WdDir*1.5

  WdSpd[0,*]=0

  U = WdSpd*sin(WdDir/180*!pi)
  V = WdSpd*cos(WdDir/180*!pi)
  x = DINDGEN(90)*16
  y = DINDGEN(45)*16

  ;绘图
  SpdImg = IMAGE(DisSpd,RGB_TABLE=39,TITLE=diplayname,grid_units=1,POSITION=[0.1,0.2,0.9,0.9])
  ;加载vector
  vec = VECTOR(COLOR='BLACK',U,V,x,y,/OVERPLOT,LENGTH_SCALE=0.0000001,HEAD_ANGLE=20,HEAD_SIZE=0.4,THICK=0.3)
  grdx= axis(TICKDIR=1,'X',LOCATION=[0,0],AXIS_RANGE=[0,360],COORD_TRANSFORM=[0,0.25],MINOR=0,MAJOR=13)
  grdy= axis(TICKDIR=1,'Y',LOCATION=[0,0],AXIS_RANGE=[-90,90],COORD_TRANSFORM=[-90,0.25],MINOR=0,MAJOR=7)
  ;COLORBAR
  c = COLORBAR(target=SpdImg,RGB_TABLE=39,TITLE='Wind Speed (m/s)',POSITION=[0.15,0.1,0.85,0.15])

  ;关闭文件释放设备号
  CLOSE , LUN
  FREE_LUN , LUN

ENDFOR
END

代码解释

MyRootDir='D:\3\'

这里定义了数据的路径,用的几井QuikScat的数据,你可以去官方下载,也可以从这里下载几幅测试数据,当然路径要改成你自己的。

filearr = file_search(MyRootDir,'*.gz',count=num);
FOR fileindex=0,num-1,1 DO BEGIN

这个代码段是对上面的路径做一个遍历,找出所有的.gz数据(也就是QuikScat的数据格式),之所以这么写是为了方便批处理,因为不可能只处理一个文件。不过如果仅仅是做测试,这里可以改为num-num只处理一个文件测试效果,效果满意再做批处理。关于这个函数的详细使用,你可以看这里

OPENR , LUN , filearr[fileindex] , /COMPRESS,/GET_LUN

这段代码是指读取某一个文件,关于LUN和GET_LUN,如果你不熟悉,可以看这里

;构造显示标题
filename=FILE_BASENAME(filearr[fileindex]);
filename=STRMID(filename,6,10)
year=STRMID(filename,0,4)+' '
month=STRMID(filename,4,2)+' '
day=STRMID(filename,6,2)+' '
diplayname=year+month+day 

这个代码快说的很清楚,就是在从文件名中构造显示标题,说到底就是一个字符串处理函数,关于IDL字符串处理的内容,你可以看这里

;读入数组&抽样
FleArr = BYTARR(1440,720,3)
READU , LUN , FleArr

DisSpd=FleArr(*,*,0)*0.2
WdSpd = congrid(FleArr(*,*,0),90,45)
WdDir = congrid(FleArr(*,*,1),90,45)

这里的代码是读数据了,这里我们用到了两个数据,一个是风速,一个是风向,同时因为不可能显示每一个点的风向(其实也可以,但是一个箭头是有大小的,那样的所有的箭头挤在一起,会有密集恐惧症的),所以这里做了一个抽样,把1440720抽样成了9045(不必惊慌,觉得够用,甚至这都会很密集)。关于抽样的函数,你可以看这里

这里可能你会感觉奇怪,为什么我创建了一个DisSpd的数组?其实应该可以猜想出来,这个数组没有进行抽样,是用来显示风速的,也就是那个色标。

;过滤数据
GT250IdxDir = where(WdDir GT 250);
WdSpd[GT250IdxDir] = 0;

GT250IdxDir = where(DisSpd GT 250*0.2);
DisSpd[GT250IdxDir] = 0;
DisSpd[GT250IdxDir] = max(DisSpd)+0.1;

这里的代码是过滤了一次数据,对于第一块代码,大多数人能看懂,因为因为大于250的是无效数据,所以赋值为0。但是第二段代码最后为什么要先赋值0,再赋值最大值+0.1呢?其实这与后面用到的色表有关系,因为我习惯用的是39号色表,也就是rainbow+white,也就是说最大值为白色,如果我不加那句的话,所有的陆地(无效数据中大部分是陆地)就变成黑色,但是会和经纬度的标注颜色重叠,所以我赋给一个最大值多一点,这样陆地就变成白色的了。

;计算U和V,X和y
WdSpd = WdSpd*0.2
WdDir = WdDir*1.5

WdSpd[0,*]=0

U = WdSpd*sin(WdDir/180*!pi)
V = WdSpd*cos(WdDir/180*!pi)
x = DINDGEN(90)*16
y = DINDGEN(45)*16

这个地方就是计算矢量四个重要参数了,U,V,X,Y。但是这里多了一个WdSpd[0,*]=0 这是为什么呢?如果去掉,你可以看到图的效果是这样的。

有的箭头就跑出去了,这里仅仅是为了出图的美观。关于U,V,X,Y的设置,如果你有不明白的地方,可以看这篇文章:关于IDL对象图形法的VECTOR函数

;绘图
SpdImg = IMAGE(DisSpd,RGB_TABLE=39,TITLE=diplayname,grid_units=1,POSITION=[0.1,0.2,0.9,0.9])
;加载vector
vec = VECTOR(COLOR='BLACK',U,V,x,y,/OVERPLOT,LENGTH_SCALE=0.0000001,HEAD_ANGLE=20,HEAD_SIZE=0.4,THICK=0.3)
grdx= axis(TICKDIR=1,'X',LOCATION=[0,0],AXIS_RANGE=[0,360],COORD_TRANSFORM=[0,0.25],MINOR=0,MAJOR=13)
grdy= axis(TICKDIR=1,'Y',LOCATION=[0,0],AXIS_RANGE=[-90,90],COORD_TRANSFORM=[-90,0.25],MINOR=0,MAJOR=7)
;COLORBAR
c = COLORBAR(target=SpdImg,RGB_TABLE=39,TITLE='Wind Speed (m/s)',POSITION=[0.15,0.1,0.85,0.15])

最后就是出图了,这里的代码跟之前的没有什么区别,主要的一个技巧就是在LENGTH_SCALE=0.0000001,这句的作用就是把箭头做一个缩放,缩放的比例设到很小,就难以看出差别,所以就像图中那样,显示的为没有长度的箭头。关于箭头的设置,推荐你看这篇文章:IDL中箭头的秘密(关于设置箭头属性)

写在最后

这就是全部了,如果还有什么意见或者建议,欢迎留言。

IDL中箭头的秘密(关于设置箭头属性) 关于HDF文件的一点概述(HDF4,HDF5)

作者:,GIS爱好者。
分享本文,请您带上本文链接
分享到:

已有 8 条评论

  1. 小良
    1#
    小良  · 2015-01-15 16:46

    测试数据哪里可以下载?

  2. dataiyang
    2#
    dataiyang  · 2015-02-24 09:29

    也想问问数据从哪里下载的啊??

  3. Carollina
    3#
    Carollina  · 2015-12-10 23:18

    有的部分代码看不清,请问截图可以再发一次么~

      1. Sailor
        Sailor  · 2015-12-11 23:09

        代码是文本,直接复制即可。

  4. rice317
    4#
    rice317  · 2016-05-23 11:51

    非常感谢分享,很有用~
    请教一下,计算U和V,X和y那部分下面两句代码有什么意义?
    WdSpd = WdSpd*0.2
    WdDir = WdDir*1.5

      1. Sailor
        Sailor  · 2016-05-23 13:50

        是QuikScat的数据说明中的参数。

  5. 王先生
    5#
    王先生  · 2016-08-12 12:53

    你好,我现在非常急需您这个画矢量图的方法,但是您文章中数据的链接不存在了,能否给我一个测试的数据呢?多谢,如果可以的话,能发我邮箱吗?感谢感谢

      1. Sailor
        Sailor  · 2016-08-12 13:31

        随便去下载一景QuickSat的数据就可以了。下载地址:http://malagis.com/surface-wind-field-data-download.html

发表评论