grads绘制-180，180范围图

峰少

大佬们打扰，求教一下，grads绘制-180 180 范围的波通量-180 0 范围内的波通量无法显示出来，怎么解决，谢谢啦，新人只有两个金币，望大佬帮忙

伽蓝鸟

你代码里有一句set lon -5 365，改成set lon -185 185，再试试？

峰少

伽蓝鸟 发表于 2021-4-10 14:54
你代码里有一句set lon -5 365，改成set lon -185 185，再试试？

您好，还是不行啊

伽蓝鸟

峰少 发表于 2021-4-10 21:08
您好，还是不行啊

那你好歹把代码贴出来吧。。这样是打算让我们盲猜嘛。。

峰少

'reinit'

* change the directory name
* u
'sdfopen F:\climatedata\uwnd.mon.mean.nc'
* v
'sdfopen F:\climatedata\vwnd.mon.mean.nc'
* height
'sdfopen F:\climatedata\hgt.mon.mean.nc'

*  gas constant
'define Ra=290'
* earth radius
'define a=6400000'
*角度制转弧度制
'define dlat = cdiff(lat,y)*3.1416/180'
'define dlon = cdiff(lon,x)*3.1416/180'

'define coslat = cos(lat*3.1416/180)'
'define sinlat = sin(lat*3.1416/180)'
* Coriolis parameter
'define f = 2*7.24/100000*sinlat'
'define f0 = 2*7.24/100000*sin(43*3.1416/180)'
'define g=9.8'

* unit [hPa] 200
'set lev 200'

* For drawing polar projection map  
'set lon 0  360'
'set lat 0 90'

* making basic state (climatology)
'define uclm = smth9((ave(u.1,time=jun1981,time=jul2010,1yr)+ave(u.1,time=jul1981,time=jun2010,1yr))/2)'
'define vclm = smth9((ave(v.2,time=jun1981,time=jul2010,1yr)+ave(v.2,time=jul1981,time=jun2010,1yr))/2)'
'define zclm =  smth9((ave(z.3,time=jun1981,time=jul2010,1yr)+ave(z.3,time=jul1981,time=jun2010,1yr))/2/g)'
* anomalies
'define zano=smth9(ave(z.3,time=jun2020,time=jul2020)/g))-zclm'
*去除纬向
'define za=zano-ave(zano,lon=0,lon=360,-b)'

* QG stream function
'define psia=g*za/f'


*psia lat mean
'define psiam=ave(psia,lon=0,lon=360,-b)'

* magnitude of basic state wind speed
'define magU = mag(uclm,vclm)'

'define dpsidlon = cdiff(psia,x)/dlon'
'define ddpsidlonlon = cdiff(dpsidlon,x)/dlon'

'define dpsidlat = cdiff(psia,y)/dlat'
'define ddpsidlatlat = cdiff(dpsidlat,y)/dlat'
'define ddpsidlatlon = cdiff(dpsidlat,x)/dlon'

'define termxu = dpsidlon*dpsidlon-psia*ddpsidlonlon'
'define termxv = dpsidlon*dpsidlat-psia*ddpsidlatlon'
'define termyv = dpsidlat*dpsidlat-psia*ddpsidlatlat'

* "p" is normalized by 1000hPa
'define coeff=coslat*(lev/1000)/(2*magU)'
*x-component
*'define px = coeff/(a*a*coslat)*( uclm*termxu+ vclm/coslat*termxv)'
* corrected on 14th Sep. 2015: Thanks to Yan Jin.
'define px = coeff/(a*a*coslat)*( uclm*termxu/coslat + vclm*termxv)'
*y-component
'define py = coeff/(a*a)*( uclm/coslat*termxv + vclm*termyv)'

'set parea 0.5 5.5  0 8.5'
'set lon  180 360'
'set lat 0 90'
'set ccolor 1'
'set cint 10'
'set gxout contour'
'set map 0'

* QG stream-function
*'set black -0.1 0.1'
*'set cthick 20'
*'set cint 2000000'
*'d maskout( (psia-psiam),  abs(lat)-10)'
'd za'


* horizontal wave-activity flux
'set arrscl 0.6 10 '
'set arrowhead -0.3'
* maskout regions where
* (i) QG approximation is supposed to be invalid (lower latitudes)
* (ii) westerly wind speed is weak or negative
*'d skip(px,2,2);maskout(maskout( py , abs(lat)-10),magU-5)'
'd skip(px,30,20);maskout(maskout( py , abs(lat)-10),magU-5)'

'set parea 5.5 10.5  0 8.5'
'set lon  0 180'
'set lat 0 90'
'set gxout contour'
'set ccolor 1'
'set mpdset lowres'
'set map 1'
'set ylab off'

* QG stream-function
*'set black -0.1 0.1'
*'set cthick 20'
*'set cint 2000000'
*'d maskout( (psia-psiam),  abs(lat)-10)'
*'d za'


* horizontal wave-activity fluza
'set arrscl 0.6 10 '
'set arrowhead -0.3'
* maskout regions where
* (i) QG approximation is supposed to be invalid (lower latitudes)
* (ii) westerly wind speed is weak or negative
*'d skip(px,2,2);maskout(maskout( py , abs(lat)-10),magU-5)'
'd skip(px,30,20);maskout(maskout( py , abs(lat)-10),magU-5)'
'set parea 0.5 10.5  0 8.5'
'set lon -180 180'
'set ccolor 1 '
'set xlab off'
'd za'


'draw title  WAF in Jun and July,meanly,2020'
'printim F:\200hPajunjulqvchuweixiang2020-TN01.png white'

峰少

我d过psia，其实在psia计算的时候就会产生问题显示不全，目前的代码是通过虚页设置画了两张图拼了起来，但是0度有坐标轴线

峰少

F:\200hPajunjulqvchuweixiang2020-TN01.png

峰少

伽蓝鸟 发表于 2021-4-10 23:39
那你好歹把代码贴出来吧。。这样是打算让我们盲猜嘛。。

麻烦大佬了

峰少

已经自己解决，问题出在了有几个变量定义之前未选择-180 180 范围，系统默认是0 360 ，自己又将问题复杂化，谢谢论坛里的大佬