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

river

李小毛123 发表于 2020-8-8 20:41
'reinit'
'open F:\fnl\fnl072000.ctl'
'set grads off'

你这个问题应该是因为在出图之前少设置一个经度范围，你不设置他默认是最前面的set x 1 和 set y 1，只有一个点，所以函数就报错了。你需要在出图之前重新设置，比如 set lon 113.8 118.2，应该就可以了

river

puck66 发表于 2020-9-13 13:28
一直没搞明白遮地形的原理，不知道论坛哪里是否有介绍？另外就是，我用gtopo30_10min.nc的地面高程数据，是 ...

遮地形的原理很简单啊，就是有地形的数据和变量的数据，画图的时候用地形的图把变量的图盖住。至于你说的那个地形资料，我不了解它的格式，不清楚可不可以用

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

谢谢分享~感谢楼主~

雨恬魅影

太牛了，非常感谢分享~