【完美】解决FNL资料任意两点间斜剖面叠加地形的问题

river

    之前发过的任意两点间斜剖面做法总结【利用GRADS自带函数】贴子有很多人发现了问题：原贴使用的资料是NCEP2.5*2.5的数据叠加NCEP2.5*2.5地形数据，很方便就能出图。因为数据本身和地形数据的水平分辨率和垂直分辨率完美的一致。    但现在大家经常使用的都是高分辨率的FNL资料，一般是1*1的。虽然利用GrADS的re()函数可以方便的将地形数据也插值成相同分辨率的地形，但是垂直层次上却还是只能保持17层（从1000hPa到100hPa）。目前最新版的FNL资料很多变量都达到了31层，可以直达1hPa的高度。这样就出现利用FNL资料绘制任意斜剖面的过程中由于垂直分辨率不同导致gr2stn函数运行过程中报错，提示 不兼容的网格。
    之前想到的是直接利用插值程序将地形资料垂直层次插值到需要的层数，但是个人能力有限，脑容量也有限，在紧张的工作之余想了半天也没有好方案。我想了一下那个NCEP资料的地形数据是用等压面去水平切割地形得到的数据，这个和平时的高空图一个道理，有地形的地方是正值，没有的地方赋负值。而且地形基本都在500hPa以下，所以理论上只需要将FNL资料多出来的等压面层次转换成相应的高度数据增加到地形数据中就可以了。说是这么说，我编程能力、时间都有限，主要还是智商跟不上，所以提供思路，这个工作只能让高人来帮忙了！
    现在说我想到的一个方案。其实比较简单，不知道大家有没有注意到：虽然FNL资料和地形资料垂直层次不同，但画正常的经度--高度或者纬度--高度剖面图的时候，叠加这个地形是没有任何问题的。也就是说只要水平分辨率相同，无论层次怎么样，正常的经纬度剖面图都是可以和这个地形无缝叠加的。所以出现上述的网格兼容问题是出在gr2stn函数本身，因为它涉及到格点插值到站点，插值过程中变量和地形垂直层次不同，造成出现了错误。那么我们换一个角度，把gr2stn的插值过程分开来，地形插值地形的，变量插值变量的，然后将插值以后的变量和地形叠加起来。这样不就和前面说的正常经纬度剖面一样了吗。
   可能有点抽象，下面上脚本：

*垂直环流图，要求有两个变量，一个是平行于高度(Z轴)的omega场，另一个是平行于斜线的风场分量。后者需要用u,v投影到斜线上。具体gs文件如下：
*顺便将地形(orog1.ctl)加入其中
'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=70.0
lon2=130.0
lat1=35.0
lat2=45.0

lon=lon1
'collect 3 free'
while(lon<=lon2)
lat=lat1+(lat2-lat1)*(lon-lon1)/(lon2-lon1)
'collect 3 gr2stn(orog,'lon','lat')'
lon=lon+1
endwhile

*再插值变量到所需斜线
'set dfile 1'
'set x 1 360'
'set y 1 181'
'set z 1 21'
'set t 3'
*批量描述FNL，多时次
'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=70.0
lon2=130.0
lat1=35.0
lat2=45.0
lon=lon1
'collect 1 free'
'collect 2 free'
while(lon<=lon2)
lat=lat1+(lat2-lat1)*(lon-lon1)/(lon2-lon1)
'collect 1 gr2stn(-VVELprs*100,'lon','lat')'
'collect 2 gr2stn(uv,'lon','lat')'
lon=lon+1
endwhile

'set grads off'
'set zlog on'
'set grid off'
'set lon 70 130'
'set z 1 21'
'set ylab off'
'set xlab off'
'set csmooth on'
*画地形
'set clevs 0'
'set ccols 0 1 0'
'set gxout shaded'
'd coll2gr(3,-u)'
*画变量
'set ylab on'
'set xlab on'
'set csmooth on'
*'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 xaxis 'lon1' 'lon2''

'set gxout stream'
'set clab on'
'd coll2gr(2,-u);coll2gr(1,-u)'
'draw title along (35N,70E) to (45N,130E)'
'gxprint g:/test/anypoumian3.png white'
say 'ok'
;



过程中没有任何报错，出来的图如下：

FNL斜剖面叠加地形

需要注意的一点就是，地形的插值过程一定要放在最前面，否则一些设置会导致报错，能出图，但是叠加起来会错位。
好啦，至此这个问题比较完美的解决了！如有任何问题，请不吝指出！



    后面有人问，是不是只能让等值线或者流线之类的和地形相同的颜色，不然的话地形的遮挡作用就没有了。这里要跟大家说一下2.0开始的grads是支持透明色的，也就是说shaded图形不一定非要放在其他图形前面画的。这个问题grads官网有很详细的说明，有能力有精力有想法的人可以去好好研究下。我今天就简单普及一下。
  这里主要用到的两个命令就是 set clevs 和 set ccols ，下面是官网介绍的一部分：

Omitting Colors

The default behavior of GrADS when plotting filled contours or shaded grid cells is to colorize all areas. To omit a particular color (or contour level) from the plot, simply assign the background color. For example:

set gxout shaded
set clevs -5 -4 -3 -2 -1 1 2 3 4 5
set ccols 0 17 18 19 20 0 21 22 23 24 0

This example is similar to the one given above, but notice where some of the ccols have been set to "0" (the background color). The first, last, and middle colors have been omitted. These commands set up a plot that will only shade areas where the anomalies are between 1 and 5 and -1 and -5. The remaining areas will be black (or white, depending on what the background color is set to be.)

When using version 2.0.0+ and 'gxout shade2' or 'gxout shade2b', if any of the color numbers is < 0, the contour is not drawn at all (i.e., it is effectively transparent（透明的）). In version 2.1+, the 'set gxout' option 'shaded' is an alias for 'shade2'.

所以利用这两个命令，把你不想被阴影图覆盖的地方的颜色设置成负数，它就会被忽略，相当于透明的效果。还得注意新版grads其实已经没有了shaded这个图形，他已经升级成了三个命令：shade1,shade2,shade2b。shade1和原来的shaded是一样的，后面的两个有什么不同请大家到官网看文档吧，但是要用透明色，必须设置后面两个才行。
下面先上图：

看到了吧，地形和流线是不同的颜色，但是覆盖的很好。
脚本在这：
'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=70.0
lon2=130.0
lat1=35.0
lat2=45.0
lon=lon1
'collect 3 free'
while(lon<=lon2)
lat=lat1+(lat2-lat1)*(lon-lon1)/(lon2-lon1)
'collect 3 gr2stn(orog,'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=70.0
lon2=130.0
lat1=35.0
lat2=45.0
lon=lon1
'collect 1 free'
'collect 2 free'
while(lon<=lon2)
lat=lat1+(lat2-lat1)*(lon-lon1)/(lon2-lon1)
'collect 1 gr2stn(-VVELprs*100,'lon','lat')'
'collect 2 gr2stn(uv,'lon','lat')'
lon=lon+1
endwhile

'set grads off'
'set zlog on'
'set grid off'
'set lon 70 130'
'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 35N,70E|37N,82E|39N,94E|41N,106E|43N,118E|45N,130E'
*'set xaxis 'lon1' 'lon2''
'set gxout stream'
'set ccolor 5'
'set clab on'
'd coll2gr(2,-u);coll2gr(1,-u)'
'draw title along (35N,70E) to (45N,130E)'

'set lon 70 130'
'set z 1 21'
'set ylab off'
'set xlab off'
'set csmooth on'
'set gxout shade2'
'set clevs 0'
'set ccols -1 3'
'd coll2gr(3,-u)'

'gxprint g:/test/anypoumian4.png white'
say 'ok'
;
看到了吧，阴影图是画在最后的，但是没有覆盖所有的流线，只是覆盖了地形的一部分。我简单说下怎么靠那两个命令实现的。
'set clevs 0'
这句是特指出一条线，0线。相当于分界线，大于这个线是一种颜色，小于的是另一个颜色。

'set ccols -1 3'
这个命令来设置颜色，根据上面说的就知道这个命令设置的颜色必须要比上面的命令设置的等值线数量多一个才行。在这里它的意思是小于0的一边全部忽略不画，或者说透明。而大于零的全部画成3这个颜色，grads里3就是绿色。

更具体的关于这两个命令的说明请大家移步官网看官方文档吧。
就到这，有问题，请不吝指出！

kmdqlb

谢谢楼主分享

独奏古筝到天明

赞一个，完美

Flash

多谢多谢

iam428jll

楼主很厉害！！！！

tulalang

效果拔群！楼主果然厉害，很久没发帖，都是在酝酿！

扫眉才子

确实完美！river大神很久没发帖了，一发帖就是好东西，感谢感谢！

migrant

这个问题比较完美的解决了！

micolian

谢谢楼主在繁忙的工作中，抽出时间帮我解决问题，在此万分感谢！

river

micolian 发表于 2017-9-19 15:13
谢谢楼主在繁忙的工作中，抽出时间帮我解决问题，在此万分感谢！

不客气，有用就好