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中箭头的秘密（关于设置箭头属性）

写在最后

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