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

puck66

river 发表于 2020-9-15 23:21
遮地形的原理很简单啊，就是有地形的数据和变量的数据，画图的时候用地形的图把变量的图盖住。至于你说的 ...

嗯嗯，照你的方法，图到是都弄出来了

瑾珂

后面没敢弄地形的弄了一个斜剖，感觉画出来的图不完全，也不知道是什么原因，画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-9-23 16:36
后面没敢弄地形的弄了一个斜剖，感觉画出来的图不完全，也不知道是什么原因，画MPV1的时候时能画出来的
...

absvprs 这是哪来的变量？

瑾珂

river 发表于 2020-9-23 21:16
absvprs 这是哪来的变量？

虽然没有再次定义，但是是fnl里有的数据
* PROFILE hPa
zdef 26 levels 100000 97500 95000 92500 90000 85000 80000 75000 70000 65000 60000 55000 50000 45000 40000 35000 30000 25000 20000 15000 10000 7000 5000 3000 2000 1000
options pascals
vars 98
no4LFTXsfc   0,1,0   0,7,193 ** surface Best (4 layer) Lifted Index [K]
no5WAVA500mb   0,100,50000   0,3,197 ** 500 mb 5-Wave Geopotential Height Anomaly [gpm]
no5WAVH500mb   0,100,50000   0,3,193 ** 500 mb 5-Wave Geopotential Height [gpm]
ABSVprs    26,100  0,2,10 ** (1000 975 950 925 900.. 70 50 30 20 10) Absolute Vorticity [1/s]

river

瑾珂 发表于 2020-9-30 16:05
虽然没有再次定义，但是是fnl里有的数据
* PROFILE hPa
zdef 26 levels 100000 97500 95000 92500 90 ...

ratio1 = (lat4-lat3)/ (lon4-lon3)
  ratio2 = -1/ratio1
这两句应该不需要放到循环里面，你拿到循环外面再试一下

瑾珂

river 发表于 2020-10-1 09:32
ratio1 = (lat4-lat3)/ (lon4-lon3)
  ratio2 = -1/ratio1
这两句应该不需要放到循环里面，你拿到循环 ...

好的，O(∩_∩)O谢谢

liluos

谢谢分享~感谢楼主~

雨恬魅影

太牛了，非常感谢分享~

chuyingjia

你好，我按照您的方法用自己的数据画剖面，一直提示错误，自己实在找不出错在哪里，希望帮忙指教一下
我的数据ctl如下：
dset     F:\2016and2018\Radar_winddata\Radar_Wind_dBz_2015.grd

undef    999.000000   

title    radar vr

xdef     1100  linear    70.0000000       5.99999987E-02

ydef     640  linear    16.0000000       5.99999987E-02

zdef     15  linear   0.500000000      0.500000000   

tdef     24  linear  00z30Jul2015  1hr

vars     3

u        15  99 99

v        15  99 99

dbz      15  99 99

endvars
我的gs文件如下：
'reinit'
'open F:\2016and2018\Radar_winddata\Radar_Wind_dBz_2015.ctl'
'gxprint F:\2016and2018\Radar_winddata\radar_wind_poumian.png white'
'set grads off'
'set zlog on'
'set csmooth on'
'set x 1'
'set y 1'
'set lev 0.5 7'
'set t 9'

lon1 = 117.88
lon2 = 119.32
lat1 = 35.56
lat2 = 37
lon = lon1

'collect 1 free'
'collect 2 free'
while (lon <= lon2)
  lat = lat1 + (lat2-lat1)*(lon-lon1) / (lon2-lon1)
  'collect 1 gr2stn(u*2.5,'lon','lat')'
  'collect 2 gr2stn(v*2.5,'lon','lat')'
  lon = lon + 0.12
endwhile
'set lon 117.88 119.32'
'set lev 0.5 7'
'set t 9'
'set ylevs 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7'
'set gxout barb'
'set cthick 5'
'd coll2gr(1,-u);coll2gr(2,-u)'
'gxprint F:\2016and2018\Radar_winddata\radar_wind_poumian.png white'
;
提示错误为：cannot draw vectors/bars--U field all undefined

river

chuyingjia 发表于 2021-1-13 17:59
你好，我按照您的方法用自己的数据画剖面，一直提示错误，自己实在找不出错在哪里，希望帮忙指教一下{:eb50 ...

你的coll2gr函数的参数应该用错了，是你没有仔细研究过后面这个-u 参数的意义。建议以后多去官网看看

coll2gr()
This function creates a grid from a collection of station data. The syntax is:

coll2gr(cnum <,num>)

cnum     collection number
num      number of vertical levels in the output grid (default is 10)
If num is given as "-u", then the vertical levels will be a union of all the levels in the collection of station data.

Usage Notes
1.Before using coll2gr, you must create the collection of stations using the GrADS command collect.

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

3.coll2gr produces an output grid that varies in X and Z; the dimension environment used when coll2gr is invoked must also be X and Z varying. The X axis of the output grid will contain the equally spaced station profiles and will span the range of the current X dimension environment. The Z axis of the output grid will span the range of the current Z dimension environment and will have either the specified number of levels or a union of the levels. Data points outside of the range of levels will be used for interpolating to within the range if appropriate.


4.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.


See the section of the User's Guide on Arbitrary Cross Sections for more information.