任意两点间斜剖面做法总结【利用GRADS自带函数】

river

   有的时候正常的经纬度剖面图不能很好地反映某个物理量的变化，需要将其沿某一不平行于经纬度的直线做剖面才能反映出其变化。这也是很多人都问的一个问题。前些日子看动力论坛上有相关的帖子，学习了一下。今天就在这儿总结一下，希望对大家有一些用处。
传送门 【利用GRADS自带函数】制作沿任意“折线”的剖面图，并叠加地形的做法

   其实主要就是利用grads自带的几个函数实现。用到的主要function函数如下：

collect 、coll2gr 、gr2stn

其中gr2stn清风大神@mofangbao已经讲的比较清楚了，简单点儿说就是把格点的值插值到选定的经纬度上，具体见清风的帖子http://bbs.06climate.com/forum.php?mod=viewthread&tid=1780。

      那另外两个函数干什么用的呢，简单说一下。collect 命令可以把站点数据和不同的时间序列的数据暂时存放在内存中作为可以调用的一个数据集，用法很简单：collect cnum expr ，cnum就是个编号，为0-31之间任意一个值，expr就是grads中任意一个表达或者表达式。

   coll2gr就是将collect在内存中的站点数据集转化成格点数据。

    上面这两个函数的解释和用法是我自己的理解。还有很多注意事项什么的，这里不细说，有兴趣的就看这个附件吧 用到的几个grads命令的说明.docx (19.13 KB, 下载次数: 1046)

2012-11-10 19:41 上传

点击文件名下载附件 。

注意事项：

1.coll2gr does not yet support time slices; currently, it will only work when the collection of stations is a collection of vertical profiles.

coll2gr函数现在还不能用来做时间的剖面（好多人都问，特此注明）

2.The "-u" option will only work for data in pressure vertical coordinates. If your vertical coordinate is height, num should be the actual number of vertical levels in the collection of profiles. Note: This restriction has been removed as of version 2.1.0.  （一定要注意这个参数，很多人都因为这个参数没搞明白而栽跟头了）

好了，进入正题，任意两点间斜剖面的做法。具体思路是这样的：

一，利用gr2stn将需要的数据插值到做斜剖面的直线上

二，利用collect将插值好的站点数据保存在内存中，供下一步调用

三，利用coll2gr将内存中的数据插值到格点，作图

其实在具体过程中第一二部是一起做的，我分开说是为了好理解。

下面给个例子，也是grads官方给的一个。我修改了一下：

*Here is a sample script written by M. Fiorino that uses these features:

*********************************************************************

* The following lines will display an arbitrary X section

* from one specified point to another.

*

* lon1 is the westernmost longitude point

* lon2 is the easternmost longitude point

* lat1 is the latitude that corresponds to lon1

* lat2 is the latitude that corresponds to lon2

*

* The loop is used to interpolate between points in

* the arbitrary cross section. This code will plot

* any cross section as long as you specify the points.

* My code plots cross sections of PV after I calculated

* PV on 11 pressure surfaces. I have another script

* that plots cross sections of potential temperature, and

* the code is very similar to this, except theta is substituted

* for PV.

*

* Many thanks to Brian Doty at COLA for his help with this code.

*

********************************************************************

'reinit'

'sdfopen f:\ncep\omega.2012.nc'

'set grads off'

'set zlog on'

'set csmooth on'

'set x 1'

'set y 1'

'set lev 1000 100'

lon1 = -95.0

lon2 = -90.0

lat1 = 55.0

lat2 = 15.0

lon = lon1

'collect 1 free'

while (lon <= lon2)

  lat = lat1 + (lat2-lat1)*(lon-lon1) / (lon2-lon1)

  'collect 1 gr2stn(omega,'lon','lat')'

  lon = lon + 1

endwhile

'set map 15 1 9'

*'set parea 1 10.2 1.5 8'

'set xlopts 1 6 0.14'

'set ylopts 1 6 0.14'

'set font 0'

'set x 14 16'

'set xaxis 'lon1' 'lon2

'set clab on'

'set gxout shaded'

'set cmax 0'

*'set ccols 0 0 7 0'

'd coll2gr(1,-u)'

'set gxout contour'

*'set cint .5'

'set cthick 11'

'd coll2gr(1,-u)'

'draw title omega along (55N,95W) to (15N,90W)'

'printim f:\test\1\section.gif white'

'c'

;

出来的图如下：这里需要注意的一点是，横轴所标示的值不是平时的经度值，他表示的是你作图所沿的那条直线上的点，也就是每一个经度对应了一个纬度。如果是站点资料表示的就是站点的经纬度。我们可以利用

set xlabs lab1 | lab2 | lab3 ...

set ylabs lab1 | lab2 | lab3 ...来强行改变它的表示。

下面再给出一个做流场图剖面的例子：

首先要想象一下，我们平时的流场图需要的平行于经度（纬度）和垂直于经度（纬度）的两个矢量。那同样的沿任意一条直线的流场图也需要两个量，一个量是平行于这条直线（可以利用U V风场投影在这条直线上），另一个量垂直于这条直线（也就是平行于高度z轴，那就是垂直速度了）

'reinit'

'sdfopen f:\2012\omega.2012.nc'

'sdfopen f:\2012\uwnd.2012.nc'

'sdfopen f:\2012\vwnd.2012.nc'

'sdfopen f:\2012\hgt.2012.nc'

'open f:\test\orog.ctl'

'set grads off'

'set grid off'

'set zlog on'

lon1=70.0

lon2=130.0

lat1=35.0

lat2=45.0

'define alfa=atan2('lat2-lat1','lon2-lon1')'

'set x 1'

'set y 1'

'set lev 1000 100'

'set t 1'

lon=lon1

'collect 1 free'

'collect 2 free'

'collect 3 free'

while(lon<=lon2)

lat=lat1+(lat2-lat1)*(lon-lon1)/(lon2-lon1)

'collect 1 gr2stn(-omega.1*100,'lon','lat')'

'collect 2 gr2stn(Uwnd.2*cos(alfa)+Vwnd.3*sin(alfa),'lon','lat')'

'collect 3 gr2stn(orog.5(t=1),'lon','lat')'

lon=lon+1

endwhile

'set lon 70 130'

'set lev 1000 100'

'set ylevs 1000 925 850 700 600 500 400 300 200 100'

'set xlabs 35N,70E|37N,82E|39N,94E|41N,106E|43N,118E|45N,130E'

'set gxout shaded'

'set csmooth on'

'set clab off'

'set clevs 0'

'set ccols 0 5 0'

'd coll2gr(3,-u)'

'set gxout stream'

'set clab on'

'd maskout(coll2gr(2,-u),-coll2gr(3,-u));maskout(coll2gr(1,-u),-coll2gr(3,-u));mag(coll2gr(2,-u),coll2gr(1,-u))'

'draw title along (35N,70E) to (45N,130E)'

'printim f:\test\1\anypoumian.gif white'

say 'ok'

;

出图如下

需要说明一下我利用NCEP的2.5°X2.5°的地形资料将地形加进去了，这个资料论坛上可以下载到的。（好吧，貌似有人一时找不到那个地形数据，我就传上来吧 NCEP地形.zip (432.6 KB, 下载次数: 2178)

2012-11-10 22:35 上传

点击文件名下载附件
阅读权限: 20 ）

另外还要注意omega负值表示上升，正值表示下降所以要记得乘以-1。还有就是这里omega的单位是Pa/s，为了和U V（单位是0.1m/s）的量级差不多就需要给omega乘以100，否则流线就会很平直。（垂直速度w和omega之间是有一个公式ω(omega)=﹣ρɡw，可以转化，那样会更好。我这个只是一个例子，有兴趣的同学自己转化一下吧）。

  综上就有了gs中的-omega.1*100

好了就这么多，一起学习一起进步。有问题的话请及时反馈，谢谢！

（2017年6月11日22时更新）以往大家使用fnl资料总是找不到合适分辨率的地形数据，大部分人都是只有NCEP的2.5*2.5的地形数据。其实利用OpenGrADS自带的re()函数很容易将其插值成1*1的。所以在此直接给大家上传吧。方便着急使用的，没有时间自己插值的同学们。

附件在此，小雨及以上级别可以下载 NCEP1x1地形.zip (3.12 MB, 下载次数: 663)

2017-9-18 08:50 上传

点击文件名下载附件
阅读权限: 20

再把插值的0.25分辨率地形传上来吧 NCEP地形0.25x0.25.zip (47.75 MB, 下载次数: 4849)

2017-9-18 08:47 上传

点击文件名下载附件
阅读权限: 20
插值使用的脚本也帖上来吧

1. 'reinit'
   
2. 'open g:/test/orog.ctl'
   
3. 'set gxout fwrite'
   
4. 'set fwrite g:/test/orog1.grd'
   
5. 'set t 1'
   
6. zz=1
   
7. while(zz<=17)
   
8. 'set z 'zz''
   
9. 'set x 1 144'
   
10. 'set y 1 73'
    
11. 'd re(orog,360,'linear',0,1,181,'linear',-90,1,bl)'
    
12. zz=zz+1
    
13. endwhile
    
14. 
    
15. 'disable fwrite'
    
16. ;

复制代码

发现的问题：贴子里提供的地形数据高度层次一共17层，这个和NCEP2.5*2.5数据的层次是完全一致的，但是和FNL资料的层次是不同的，所以你使用FNL资料来画任意直线剖面并叠加地形的时候就会报错（需要注意的是使用FNL资料画正常的经纬度剖面叠加这个地形的时候是可以用的，我觉得是这种剖面不涉及插值，所以可以用），没办法出图。但目前我还没有发现比较好的可以在垂直层次上给地形数据插值的程序，有人有的话是否可以提供一下。或者发帖分享一下！谢谢（已经解决，可以去看我的另一个贴【完美】解决FNL资料任意直线剖面叠加地形的问题）

cathy

river 发表于 2013-11-23 11:59

                               
登录/注册后可看大图

你去看一下coll2gr的官方文档，上面有详细用法和注意事项，我看你的资料好像不是等高面上的资料，后面-u的 ...

谢谢@river  ！问题解决了，只有等压面上的资料可以用-u，若等高面上的资料则需要改为所用等高面的层数，比如此处的26 。

YUNER

没什么技术含量,给有这方面疑问的提供个方便。

瑾珂

后面没敢弄地形的弄了一个斜剖，感觉画出来的图不完全，也不知道是什么原因，画MPV1的时候时能画出来的
'reinit'
'open L:\DZQDL\2014datafiles\fnl_2014.ctl'
i=1
while (i<=31)
'set t 'i''
'set lon 80 120'
'set lat 10 40'
'set z 1 21'
'define tmpz=TMPprs'
'define rhz=RHprs'
'define prsz=lev'
'define esz=(6.112*exp(17.67*(tmpz-273.15)/(tmpz-29.65)))'
'define qz=rhz*(0.62197*esz/(prsz-esz))/100.'
'define ez=prsz*qz/(0.62197+qz)+1e-10'
'define tlclz=55.0+2840.0/(3.5*log(tmpz)-log(ez)-4.805)'
'define thetaz=tmpz*pow((1000/prsz),(0.2854*(1.0-0.28*qz)))'
'define eqtz=thetaz*exp(((3376./tlclz)-2.54)*qz*(1.0+0.81*qz))'
'define dtx=cdiff(eqtz,x)'
'define dty=cdiff(eqtz,y)'
'define dx=cdiff(lon,x)*3.1416/180*cos(lat*3.1416/180)*6.371e6'
'define dy=cdiff(lat,y)*3.1416/180*6.371e6'
'define vy=VGRDprs'
'define uy=UGRDprs'

'define eqtxy=eqtz(z-1)-eqtz(z+1)'
'define vxy=vy(z-1)-vy(z+1)'
'define uxy=uy(z-1)-uy(z+1)'

'define dp=10000'
'define mp1=-9.8*absvprs*eqtxy/dp'
'define mp2=9.8*((vxy/dp)*(dtx/dx)-(uxy/dp)*(dty/dy))'
'define mp=mp1+mp2'
'define mpv=mp*1E6'
'define mpv1=mp1*1E6'
'define mpv2=mp2*1E6'

'set x 1'
'set y 1'
'set z 1 21'
lon1 = 96.0
lon2 = 108.0
lat1 = 21.5
lat2 = 28.6
lon3 = 101.12
lon4 = 102.04
lat3 = 26.64
lat4 = 25.09
lon = lon1
'collect 1 free'
while (lon <= lon2)
  ratio1 = (lat4-lat3)/ (lon4-lon3)
  ratio2 = -1/ratio1
  lat = lat4 + ratio2*(lon-lon4)
  'collect 1 gr2stn(mpv2,'lon','lat')'
  lon = lon + 1
endwhile

'set grads off'
'set zlog on'
'set grid off'
'set lon 96 108'
'set z 1 21'
'set ylab off'
'set xlab off'
'set csmooth on'
'set ylab on'
'set xlab on'
'set csmooth on'
*'set ylevs 1000 925 850 700 600 500 400 300 200 100'
'set xlabs 22N,96E|23N,99E|25N,102E|27N,105E|29N,108E'
'set csmooth on'
'set gxout stream'
'set ccolor 1'
'set clab on'
'd coll2gr(1,-u)'
'q dim'
it=sublin(result,5)
it=subwrd(it,6)
'draw title 'it' along (22N,96E) to (29N,108E)'
'printim L:\DZQDL\2014datafiles\mpv2_xpou\fnl_mpv2-xpou4_'it'.png white'
'c'
i=i+1
endwhile
'disable gxprint'
*

river

--惊鸿- 发表于 2020-3-1 18:47
啊啊啊啊把-u改成40一下子就出图了！可以了！
太感谢您了，这个斜剖面图困扰了我好久，真的太谢谢你的指 ...

没事儿多去GrADS官网逛逛，了解了解命令和和函数的用法，说不准就不要困扰喇嘛久了

river

爱喝柠檬汁 发表于 2017-12-29 22:00

不好意思，年底很忙。我没时间给你一一解释你问的这些问题。这些都是基本知识，你应该直接去看书和官方文档······

river

晚晴山 发表于 2012-11-10 21:41

                               
登录/注册后可看大图

你为何不把orog.ctl传上来~~~~

这个论坛的资料站有下载哦

蓝冰草

海蓝蓝 发表于 2014-11-28 12:47
把1000 925 ··· 500 ··· 10的图按照高度放起来，得到的三维数据不就是等压面的位势高度吗？
可是 ...

http://bbs.06climate.com/forum.php?mod=viewthread&tid=2663 来自luoziwuhui的贴子。我画了其中topo10.ctl所对应的数据，应该就是下垫面地形的高度。

freekiller

先赞一个再看   

随缘

神帖子！鉴定完毕！

晚晴山

恩，不错，我很喜欢~

晚晴山

你为何不把orog.ctl传上来~~~~

晚晴山

river 发表于 2012-11-10 21:57

                               
登录/注册后可看大图

这个论坛的资料站有下载哦

这................

liuyf2014

好~~~~~~~~~~~~~~~~~~

mofangbao

真的不错！会成为经典的！

river

mofangbao 发表于 2012-11-10 22:26

                               
登录/注册后可看大图

真的不错！会成为经典的！

还好你帮我移过来了，太感激了