【求助】用grads画站点风，怎么也画不出来

瑾珂

站点数据为风速和风向，我用Fortran转为二进制：
Program sta2grd
Implicit none
Character*8 stid
Real lon,lat,wd,ws,u,v,tim !这里是程序的变量声明
integer::nlev,flag
!变量声明结束

!程序开始
tim=0.0
nlev=1
flag=1
Open(1,file='L:\DJBY\feng\2017010308w1.txt',status='old')
open(2,file='L:\DJBY\feng\2017010308w1.dat',status='replace',form='binary')

10 Read(1,*,end=100)lon,lat,stid,wd,ws
u=(-1)*sin(wd*3.14/180)*ws
v=(-1)*cos(wd*3.14/180)*ws
Write(2)stid,lon,lat,tim,nlev,flag,u,v
Goto 10
100 continue
Close(1)
nlev=0
flag=0
Write(2)stid,lon,lat,tim,nlev,flag
close(2)
!程序结束
End

再编写ctl:
dset L:\DJBY\feng\2017010308w1.grd
dtype station
stnmap L:\DJBY\feng\2017010308w1.map
undef 999.9
title wind
tdef 1 linear JAN2017 1mo
vars 2
u  0  99  U-Component of Wind [m/s]
v  0  99  V-Component of Wind [m/s]
ENDVARS

编写gs:
'reinit'
'open L:/DJBY/feng/2017010308.ctl'


'set lat 20 30'
'set lon 96 108'
'set ylint 2'
'set xlint 2'
'set ylopts 1 4 0.18'
'set xlopts 1 4 0.18'
'set mpdset cnworld'
'set map 1 1 9'
'set csmooth on'

'set gxout barb'
'set ccolor 2'
'd u;v'

'printim L:\DJBY\feng\windrain2017010308.png white'
'disable gxprint'
；

始终没有结果，不知道是为什么

edwardli

还是我来吧。“格点”和“站点”是不同的。常规熟悉命令都是针对格点的，站点是离散的。
set   gxout   model
可以绘制站点填图。

灵心流萤

没结果是啥意思呀？也没报错？disable gxprint这是啥命令啊？

瑾珂

edwardli 发表于 2022-3-10 16:18
还是我来吧。“格点”和“站点”是不同的。常规熟悉命令都是针对格点的，站点是离散的。
set   gxout   mod ...

还是不能出图，是不是我转二进制的那个程序错了

瑾珂

灵心流萤 发表于 2022-3-11 08:18
没结果是啥意思呀？也没报错？disable gxprint这是啥命令啊？

就是画出来的图是空白的

灵心流萤

瑾珂 发表于 2022-3-11 09:19
就是画出来的图是空白的

你试试直接open ctl文件看能不能画出图，不能的话就是grd有问题，或者ctl没写对。或者再试试先把计算好的站点数据生成txt文件，检查下没问题再写一个程序生成grd，一步一步检查看看哪里有问题

栤唥de嶶笑

1、gxout barb肯定不行，barb是针对格点的。

2、可以把gxout后面的改成以下任一试试看，能不能出东西。例如：set gxout stnmark

For station data, these additional graphics output types are also available:

findstn    Finds nearest station
model      Plots station model
stnmark    Plots a mark at station location that is colorized by data value
tserbarb   Time series of wind barbs at a point (1-D)
tserwx     Time series of weather symbols at a point (1-D)
value      Plots station values
wxsym      Plots weather symbols at station, controlled by set wxopt

3、如果你非要画格点的barb，可以试试分别给U风和V风进行插值运算，即用Grads自带函数oacres()，最后display中， d oacres(grd.1,u.2);oacres(grd.1,v.3)，这个grd.1是需要自己构造的网格（分辨率自定），说来话长。。。

瑾珂

灵心流萤 发表于 2022-3-14 08:54
你试试直接open ctl文件看能不能画出图，不能的话就是grd有问题，或者ctl没写对。或者再试试先把计算好的 ...

生成的txt，打开是乱码
program micaps_diamond

implicit none

integer::time,nlev,nflag   ! nb is the number of stations of file, time, nlev,
                               ! and nflag are the ones that will be used while writing grd
integer::i


character*8,allocatable::id(:)
character*8,allocatable::sta(:)
real,allocatable::dir(:)
real,allocatable::speed(:)
real,allocatable::lat(:)
real,allocatable::lon(:)
real,allocatable::u(:)
real,allocatable::v(:)
real,allocatable::ncld(:)




open(11,file="L:/DJBY/feng/wind2017010308.txt",status="old",form="formatted")

   read(11,*)
   read(11,*)
   read(11,*)
  
   
   allocate(id(426))
   allocate(lat(426))
   allocate(lon(426))
   allocate(dir(426))
   allocate(speed(426))
   allocate(u(426))
   allocate(v(426))
   allocate(sta(426))
   allocate(ncld(426))


  
   do i=1,426
      read(11,*) lon(i),lat(i),id(i),sta(i),dir(i),speed(i),ncld(i)
          u(i)=(-1)*sin(dir(i)*3.14/180)*speed(i)
          v(i)=(-1)*cos(dir(i)*3.14/180)*speed(i)
   enddo
  
      
close(11)

open(21,file="L:\DJBY\feng\2017010308wind2.txt",status="replace",form="binary")


   time=1
   nlev=1
   nflag=1
   do i=1,426
      write(21) id(i),lon(i),lat(i),time,nlev,nflag,u(i),v(i)
   enddo
   nlev=0
   write(21) lon(i-1),lat(i-1),id(i-1),time,nlev,nflag
      
close(21)

end

灵心流萤

瑾珂 发表于 2022-3-16 09:35
生成的txt，打开是乱码
program micaps_diamond

你生成的txt是二进制的，当然不能直接查看呀，我的意思是你输出一个十进制的u.txt和v.txt，看看有没有问题