【求助】利用气象局日尺度数据做站点对比时不能顺利生成sta.map

shangbn305

想做站点对比分析，我在国家气象局下载的12年6月份天尺度气温文件，也就是每个站点都有30天的数据，我想对比每天的平均气温。

第一步：生成sta.grd文件，我的程序如下，生成sta,grd

      program sta2grd
      implicit none
      character*8 stid
      real H,avg,max,min,na,nb,nc,tim,lat,lon
      integer nlev,flag,year,month,day
      tim=0.0
      nlev=1
      flag=1
      open(1,file='TEM-201206.TXT',status='old')
      open(2,file='sta.grd',status='replace',form='binary')
10    read(1,*,end=100)stid,lat,lon,H,year,month,day,avg,max,
     &min,na,nb,nc
      lat=lat/100+(lat-(lat/100)*100)*0.0166666666666667
      lon=lon/100+(lon-(lon/100)*100)*0.0166666666666667
      avg=avg/10
      max=max/10
      min=min/10
      write(2)stid,lat,lon,tim,nlev,flag,avg
      goto 10100 continue
      close(1)
      nlev=0
      write(2)stid,lat,lon,tim,nlev,flag
      close(2)
      end
第二步：为生成的sta.grd文件编写sta.ctl文件，如下：对于T我使用两种方法定义，第一种方法可以顺利进行第三步，第二种方法出错DSET /public/home/shangbingnan/station/Temperature/sta.grd
DTYPE station
STNMAP /public/home/shangbingnan/station/Temperature/sta.map
UNDEF 32766
TITLE Temperature Data Sample
TDEF 1 linear 00jun2012 1dy         
TDEF 30 linear 00jun2012 1dy
VARS 1
avg 0 99 Temperature Data
ENDVARS
第三步：使用grads为站点的二进制文件生成站点映射文件sta.map!stnmap -i ./sta.ctlsta.ctl中第一个T定义可以顺利生成sta.map，信息如下：


第二个T定义出错，错误信息如下：


疑问：时间T应该是三十天啊，因为我的数据是一个月的，为什么会出错呢？

shangbn305

最后那个图和前边的图重复了

lqouc

在人人上貌似回复不了这么多，你在这里看看吧，不一定能用，可能还要你稍微改改。
      do i=1,99999
      read(1,*,iostat=it)stid,lat,lon,H,year,month,day,avg,max,
      &min,na,nb,nc
      if(it=0) exit
       lat=lat/100+(lat-(lat/100)*100)*0.0166666666666667
       lon=lon/100+(lon-(lon/100)*100)*0.0166666666666667
       avg=avg/10
       max=max/10
       min=min/10
       write(2)stid,lat,lon,tim,nlev,flag,avg
       goto 10100 continue
       close(1)
       end do

shangbn305

lqouc 发表于 2013-3-6 16:20

                               
登录/注册后可看大图

在人人上貌似回复不了这么多，你在这里看看吧，不一定能用，可能还要你稍微改改。
      do i=1,99999
   ...

我算了一下，因为站台是194个，每个站台是30个即30天的数据，所以总行数是194*30=5820行。修改后重新运行后还是同样的错误。
如果可以将错就错就好了可是这个气温是每天2:00,8:00,14:00,20:00测的，取平均。我的wrf结果使用ARWpost处理，输出的结果和这个时间相对应的，所以wrf结果的ctl文件内的T定义是tdef 5 linear 02Z01JUN2012 360MN，这个定义和grid.ctl的不统一，还是会出错

mofangbao

请把站点和时间的循环分开写，所有站点的第N个时次结束后都要nlev=0 再write一个标记的，你的程序只最后write了一次标记，所以生成map文件的时候会被认为是1个时次，5820个站点，还不明白就仔细看一下我那个站点数据作图的帖子

river

楼上正解

shangbn305

mofangbao 发表于 2013-3-6 23:10

                               
登录/注册后可看大图

请把站点和时间的循环分开写，所有站点的第N个时次结束后都要nlev=0 再write一个标记的，你的程序只最后wri ...

明白了，恍然大悟，怪我自己不够认真，谢谢

mofangbao

shangbn305 发表于 2013-3-7 08:19

                               
登录/注册后可看大图

明白了，恍然大悟，怪我自己不够认真，谢谢

把你改完的也顺便贴过来吧

shangbn305

mofangbao 发表于 2013-3-6 23:10

                               
登录/注册后可看大图

请把站点和时间的循环分开写，所有站点的第N个时次结束后都要nlev=0 再write一个标记的，你的程序只最后wri ...

  下边是我修改后的程序，以及sta.ctl，可以生成30个时次的sta.map，但是作图的时候出问题了，
作图程序sta.gs:
'open c:\grads20\win32\sta.ctl'
'set lat 33 60'
'set lon 105 148'
'set t 1'
'd avg'
>run c:\grads20\win32\sta.gs
错误信息：
No hardcopy metafile open
all filese closedl all defined objects released;
all grads attributes have been reinitialized

fortran程序：
      program sta2grd
      implicit none
      character*8 stid
      real H,avg,max,min,na,nb,nc,tim,lat,lon
      integer nlev,flag,year,month,day
      integer i,j
      open(1,file='TEM-201206.TXT',status='old')
      open(2,file='sta.grd',status='replace',form='binary')
      
      do i=1,194
      tim=0.0
      flag=1   
      do j=1,30
      nlev=1   
      read(1,*)stid,lat,lon,H,year,month,day,avg,max,
     &min,na,nb,nc
      lat=lat/100+(lat-(lat/100)*100)*0.0166666666666667
      lon=lon/100+(lon-(lon/100)*100)*0.0166666666666667
      avg=avg/10
      max=max/10
      min=min/10
      write(2)stid,lat,lon,tim,nlev,flag,avg   
      end do
      nlev=0
      write(2)stid,lat,lon,tim,nlev,flag
      end do
      
      close(1)
      
      close(2)
      end

这是我的sta.ctl
DSET C:/GrADS20/win32/sta.grd
DTYPE station
STNMAP C:/GrADS20/win32/sta.map
UNDEF 32766
TITLE Temperature Data Sample
TDEF 30 linear 01jun2012 1dy
VARS 1
avg 0 99 Temperature Data
ENDVARS

mofangbao

抱歉 那个不是错误信息，正常运行都会有的，你翻译一下看看？