【利用GRADS自带函数】制作沿任意“折线”的剖面图，并叠加地形的做法

river

说正题之前先推荐自己的另两个帖子，因为下面要用到的地形资料在这两个帖子里有：①任意两点间斜剖面做法总结【利用GRADS自带函数】②【完美】解决FNL资料任意两点间斜剖面叠加地形的问题

       正题：利用GrADS自带函数制作沿任意折线剖面图，并叠加地形的做法
    近期在瞎捉摸一些事情，发现GrADS自带的三个函数（collect 、coll2gr 、gr2stn）可以比较简单的完成任意折线的剖面图。
    collect 函数可以把站点数据和不同的时间序列的数据暂时存放在内存中作为可以调用的一个数据集，用法比较简单：collect cnum expr ，cnum就是个编号，为0-31之间任意一个值，expr就是grads中任意一个表达或者表达式。
   
    gr2stn 函数就是把格点资料插值到你要的经纬度上或者说站点上，用法：gr2stn(grid_expr, stn_expr, <-n>) or gr2stn(grid_expr, lon, lat, <-n>)

  

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

       其他的不多说了，可以去看我另两个帖子或者去看官方文档。

具体思路：

一，利用gr2stn将需要的数据插值到任意直线上，重复多次可以插值到多条直线上。

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

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

在具体过程中第一二部是一起做的，gr2stn一边插值数据到站点，collect一边把插值好的数据存放在内存中作为一个数据集。

其中一个要点就是为了防止其他的干扰，一般在正式插值建立数据集前会使用 'collect num(0-31任意整数） free '  来建立一个空数据集，然后再往里写数据。为了将不同直线上的数据放在同一个数据集中，这个命令在最早一次插值前使用，后面再插值到其他直线上时就不需要了。最后给个例子大家来体会吧。

脚本如下：

*垂直环流图，要求有两个变量，一个是平行于高度(Z轴)的omega场，另一个是平行于斜线的风场分量。后者需要用u,v投影到斜线上。

*顺便将地形(orog1.ctl)加入其中,具体gs文件如下：

'reinit'

'open G:\ncep\2017/201702.ctl'

'open G:\test\orog1.ctl'

'set dfile 2'

'set t 1'

'set x 1'

'set y 1'

'set z 1 17'

*先插值第一条斜线上的地形数据

lon1=140.0

lon2=160.0

lat1=60.0

lat2=70.0

lon=lon1

'collect 3 free'

while(lon<=lon2)

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

'collect 3 gr2stn(orog.2,'lon','lat')'

lon=lon+1

endwhile

*再插值第二条斜线上的数据

lon1=75.0

lon2=100.0

lat1=30.0

lat2=40.0

lon=lon1

*'collect 3 free'  注意这句已经注释掉了，这样就把第二条斜线上的数据增加到数据集3 中了

while(lon<=lon2)

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

'collect 3 gr2stn(orog.2,'lon','lat')'

lon=lon+1

endwhile

'set dfile 1'

'set x 1 360'

'set y 1 181'

'set z 1 21'

'set t 3'

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

*atan2返回的结果单位是弧度

'define uv=UGRDprs*cos(alfa)+VGRDprs*sin(alfa)'

'set x 1'

'set y 1'

'set z 1 21'

lon1=140.0

lon2=160.0

lat1=60.0

lat2=70.0

lon=lon1

'collect 1 free'

'collect 2 free'

while(lon<=lon2)

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

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

'collect 2 gr2stn(uv,'lon','lat')'

lon=lon+1

endwhile

lon1=75.0

lon2=100.0

lat1=30.0

lat2=40.0

lon=lon1

*'collect 1 free'

*'collect 2 free'

while(lon<=lon2)

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

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

'collect 2 gr2stn(uv,'lon','lat')'

lon=lon+1

endwhile

'set parea 1.5 10.5 0.5 7.5'

'set grads off'

'set zlog on'

'set grid off'

'set lon 70 180'

'set z 1 21'

*'set xlab off'

'set csmooth on'

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

'set xlabs 60N,140E||||70N,160E(35N,70E)|||||40N,100E'

*先画地形阴影图

'set clevs 0'

'set ccols 0 1 0'

'set gxout shade2b'

'd coll2gr(3,-u)'

*再画折线上的全风速

'set gxout contour'

'set cint 30'

*'set clab on'

*'set strmden 5 0.5 0.05 1'

'd smth9(mag(coll2gr(2,-u),coll2gr(1,-u)))'

*标记两条斜线的分界

'drawline1 Lon 120'

'draw title along (60N,140E) to (70N,160E)\and (35N,75E) to (40N,100E)'

'gxprint g:/test/anypoumian3(3).png white'

say 'ok'

;

我单独画了两条斜线上的剖面图来和上面合成一起的图做对比，验证了相关函数功能确实可以做到任意折线剖面图的合成。图如下：

来自群组: 南气院校友总会

ayzqs

强！鼓掌鼓掌

carlye2002

太强了，学习啦

Zerlinda

多谢大大分享

hm_others

请教一下，地形数据是什么样子的，如果是海拔高度，怎么解决纵坐标是百帕和米的关系呢？

river

hm_others 发表于 2021-11-27 17:04
请教一下，地形数据是什么样子的，如果是海拔高度，怎么解决纵坐标是百帕和米的关系呢？

没懂你的意思。你是说你的数据资料不是等压面的资料，是等高面上的？

puck66

这次一定要学会

puck66

river 发表于 2021-11-27 21:12
没懂你的意思。你是说你的数据资料不是等压面的资料，是等高面上的？

楼上似乎是想问这里地形的单位是啥子？如果地形是海拔高度，单位是米或者千米，绘图的时候，变量本身的纵坐标是hPa，那地形和变量在纵坐标方向的单位怎么协调？我看过一个帖子是把海拔高度换算成等压面hPa去匹配，还有一种处理办法是在图形右侧给海拔高度增加一个纵坐标，就相当于是双Y轴，后者要简便一点

river

puck66 发表于 2022-1-20 10:29
楼上似乎是想问这里地形的单位是啥子？如果地形是海拔高度，单位是米或者千米，绘图的时候，变量本身的纵 ...

这个地形资料本身就已经是等压面的了，所以直接用就好了

蒙同学_/╱°.

给力！最终还是grads承担了所有，希望有一天python版本的也有大神分享一下！