任意垂直剖面的地形画法

黎大页子

GrADs自带函数可以完成任意剖面的作图，主要依靠 collect 和 gr2stn 函数。
具体方法 river 已经有非常好的总结：http://bbs.06climate.com/forum.php?mod=viewthread&tid=10966

在垂直剖面有地形的时候，需要自己添加。river用的是下载的地形资料。但是在没有额外地形资料，只能靠wrfoutput的HGT变量画地形。
网上有一种办法：沿切面逐点读出地形高度，用 sublin 和 subwrd 读取，然后利用 q w2xy 查找该高度在虚页上对应的(x,y)坐标。最后用多边形（draw polyf）连接相邻两点的山高的(x,y)点。即可以做出连续的遮蔽的地形。
http://cookbooks.opengrads.org/i ... rary_Cross_Sections
详见链接页面里的fuction函数。

但是以上方法需要在做变量和地形时统一固定虚页的大小，否则两张图不能对应。并且wrfoutput里面的地形变量是HGT，只有一层（即一定要set z 1才能读出）。所以我延续以上思路，稍作修改。利用draw line 函数，在线段足够粗，步长足够小的情况下，足以画出连续的地形。并且不需要再人为设定虚页大小，改为利用q w2xy自动查找。
以下是脚本中画好变量后画地形的部分。代码很简单，希望能给有需要的人一些参考。

*说明：
*lon1，lat1---剖面起始经纬度
*lon2，lat2---剖面结束经纬度
*HGT--wrfoutput的地形
*0.04858 是z=1时的高度，可根据自己的资料修改。
*0.002是画地形时的步长，可自行调整。
*---------------------------------------------
  ln=lon1
  i=1
  while(ln<=lon2)
*------read terrain height--------------------
    lt=lat1+(lat2-lat1)*(ln-lon1)/(lon2-lon1)
    'set lon 'ln
    'set lat 'lt
    'set z 1'
    'd HGT*0.001'
    res1=sublin(result,2)
    top=subwrd(res1,4)

    'q w2xy 'ln' 0.04858'
    res2=sublin(result,1)
    tx=subwrd(res2,3)
    ty=subwrd(res2,6)

    'q w2xy 'ln' 'top
    res3=sublin(result,1)
    mt=subwrd(res3,6)

    if(mt>ty)
      'set line 1 1 6'
      'draw line 'tx' 'ty' 'tx' 'mt
    endif

    ln=ln+0.002
    i=i+1
  endwhile
*---------------------------------------------------------------
成图如下：

river

大叶子又来发帖，赶紧前来支持！

sjysjxyf66

无私分享，大赞！学习了，谢谢啊！

尽头的尽头

赞一个

topmad

来顶大叶子....未曾谋面  久闻大名

传说中的谁

在满屏的提问求助帖中，楼主的作品有如出水之芙蓉，如此清新脱俗，沁人心脾

kongfeng0824

谢谢分享                             

mofangbao

大页子比大A厉害啊

萝卜和白菜

学习了~谢谢论坛的大神们

黎大页子

river 发表于 2014-4-22 06:30
大叶子又来发帖，赶紧前来支持！

哈哈，多谢鼓励。你的那份总结帮了我不少，现在的年轻人都这么厉害了啊， 我等真是要被拍死在沙滩上了。