【求助】多时次站点资料作图生产映射文件发生错误

shirlyqiu · 发表于 2017-12-22 22:13:18

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

根据清风大大的《GrADS站点文件作图详细解决方案》，生成了.grd文件，但是在生成.map文件时出现如下错误：Invalid station hdr found in station binary file.

QQ截图20171222215807.png

各种看解决方案，翻了一整天的论坛，大家说的解决方案都试过了还是不行，实在无解了，po出来请大神瞅瞅，在此先行谢过！！！

Fortran program（为了看着简单，已经把实际文件名改简单）：
Program sta2grd
Implicit none
!!!!!Define variables
character*30 filename1,filename2,filename3,filename4
Character*8 stid(1000), sta, staid(1000) !station id
Real lon(1000),lat(1000),tim,tmax(366,1000),tmin(366,1000)
real temp_lon,temp_lat,t1,t2
integer nlev,flag,tt,i,j,m,n,ii,jj,day,l,a,k,z,day0 !nlev 表示总层次，地面变量算1层，有高空层次在加上；flag表示有无地面数据，有为1无为0
integer sta_num,line !sta_num表示每一年的站点总数，line表示每一年数据的总行数
integer year,month,dat,temp_year,temp_mon,temp_dat !判断时间是否连续是否需要加缺测
real undef
undef = -999.0
nlev=1
flag=1

open(14,file='station_id.txt',status='old')
      read(14,*)(staid(k),k=1,sta_num)
sta_num=130
line=28234
year=1980
if(mod(year,4)==0)then !闰年和非闰年的天数设定
tt=366
else
tt=365
endif
!!!!Read input data
Open(1,file='1980.txt',status='old')
read(1,*)!!!!跳过第一行无效数据
i=1
m=1 !站点
n=1 !时间
do while(i<line)
read(1,*)sta,temp_lat,temp_lon,temp_year,temp_mon,temp_dat,t1,t2 !读入临时变量中看是否需要加缺测值
if(mod(year,4)==0)then
if(temp_mon==1)then
day=temp_dat !得到是year的第几天
else if(temp_mon==2)then
day=31+temp_dat
else if(temp_mon==3)then
day=60+temp_dat
else if(temp_mon==4)then
day=91+temp_dat
else if(temp_mon==5)then
day=121+temp_dat
else if(temp_mon==6)then
day=152+temp_dat
else if(temp_mon==7)then
day=182+temp_dat
else if(temp_mon==8)then
day=213+temp_dat
else if(temp_mon==9)then
day=244+temp_dat
else if(temp_mon==10)then
day=274+temp_dat
else if(temp_mon==11)then
day=305+temp_dat
else if(temp_mon==12)then
day=335+temp_dat
endif
else
if(temp_mon==1)then
day=temp_dat !得到是year的第几天
else if(temp_mon==2)then
day=31+temp_dat
else if(temp_mon==3)then
day=59+temp_dat
else if(temp_mon==4)then
day=90+temp_dat
else if(temp_mon==5)then
day=120+temp_dat
else if(temp_mon==6)then
day=151+temp_dat
else if(temp_mon==7)then
day=181+temp_dat
else if(temp_mon==8)then
day=212+temp_dat
else if(temp_mon==9)then
day=243+temp_dat
else if(temp_mon==10)then
day=273+temp_dat
else if(temp_mon==11)then
day=304+temp_dat
else if(temp_mon==12)then
day=334+temp_dat
endif
endif
if(m==1 .and. n==1)then !第一行有效值全部赋值
stid(m)=sta
lat(m)=temp_lat
lon(m)=temp_lon
if(n==day)then !如果时间是连续的，那直接给tmax和tmin赋值
tmax(n,m)=t1
tmin(n,m)=t2
!write(*,*) day,n,m,stid(m),lon(m),lat(m),tmax(n,m),tmin(n,m)
else !如果时间不连续，那就把中间不连续的tmax和tmin都赋值为undef，并把n连接到更新的天数
do jj=n,(day-1)
tmax(jj,m)=-999.0
tmin(jj,m)=-999.0
!write(*,*) day,n,m,stid(m),lon(m),lat(m),tmax(n,m),tmin(n,m)
enddo
n=day
tmax(n,m)=t1
tmin(n,m)=t2
endif
n=n+1
day0=day-1
write(*,*)m,day0
else !第二行以后的情况
if(sta/=stid(m))then !如果是新站点，则站点数+1，时间n重新从1开始
if(day0 < tt)then!如果是新站点，且前一天的天总数不等于year的总天数，那说明前一个站点不是一年连续，需要把剩余天数设为undef
do jj=n,tt
tmax(jj,m)=-999.0
tmin(jj,m)=-999.0
enddo
day0=0
endif
m=m+1
n=1
stid(m)=sta
lat(m)=temp_lat
lon(m)=temp_lon
if(n==day)then !如果时间是连续的，那直接给tmax和tmin赋值
tmax(n,m)=t1
tmin(n,m)=t2
else !如果时间不连续，那就把中间不连续的tmax和tmin都赋值为undef，并把n连接到更新的天数
do jj=n,(day-1)
tmax(jj,m)=-999.0
tmin(jj,m)=-999.0

enddo
n=day
tmax(n,m)=t1
tmin(n,m)=t2
write(*,*) n
endif
n=n+1
day0=day-1
else !如果不是新站点，则不需要再给站点的经纬度等赋值，直接判断tmax和tmin
if(n==day)then
tmax(n,m)=t1
tmin(n,m)=t2
else
do jj=n,(day-1)
tmax(jj,m)=-999.0
tmin(jj,m)=-999.0
enddo
n=day
tmax(n,m)=t1
tmin(n,m)=t2
endif
!write(*,*) day,n,m,stid(m),lon(m),lat(m),tmax(n,m),tmin(n,m)
n=n+1
day0=day-1
if(m==sta_num)then
write(*,*)m,n-1,tmax(n-1,m),tmin(n-1,m)
endif
endif
endif
i=i+1
enddo
if(m==sta_num .and. day0<tt)then
write(*,*)m,n,tmax(n-1,m),tmin(n-1,m)
do jj=n,tt
tmax(jj,m)=-999.0
tmin(jj,m)=-999.0
write(*,*)m,n,tmax(jj,m),tmin(jj,m)
enddo
day0=0
endif
close(1)
!!!replace sta_id
do z=1,sta_num
stid(z)=staid(z)
enddo
            !!!!以上可以忽略不看，确认过数据处理没有问题
!!!write date
open(2,file='1980.grd',status='replace',form='binary')
do l=1,tt
tim=0.0
nlev=1
flag=1
do a=1,sta_num
write(2)stid(a),lat(a),lon(a),tim,nlev,flag,tmax(l,a),tmin(l,a)
enddo
nlev=0
write(2)stid(sta_num),lat(sta_num),lon(sta_num),tim,nlev,flag
enddo
close(2)
open(4,file="1980_test.txt",status='replace')
do l=1,tt
tim=0.0
nlev=1
flag=1
do a=1,sta_num
write(4,*)stid(a),lat(a),lon(a),tim,nlev,flag,tmax(l,a),tmin(l,a)
enddo
nlev=0
write(4,*)stid(sta_num),lat(sta_num),lon(sta_num),tim,nlev,flag
enddo
close(4)
!!!write ctl file
close(12)
close(13)
end

已结仔仔细细各种write验证了数据没问题，把grd文件格式按TXT文件打印出来也没有问题。
ctl：
dset ^1980.grd
dtype station
stnmap ^1980.map
* OPTIONS big_endian
undef -999.0
title station tmax tmin
tdef 366 linear 01jan1980 1dy
vars 2
tmax 0 99 daily maximum temperature
tmin 0 99 daily minimum temperature
endvars

求问哪里出了问题，太感谢了！

蜻蜓队长 · 发表于 2017-12-23 09:55:32

请问你的代码为何排列如此整齐?是用的什么编辑器吗？已帮忙看完程序，未发现问题所在

shirlyqiu · 发表于 2017-12-25 10:13:29

蜻蜓队长发表于 2017-12-23 09:55
请问你的代码为何排列如此整齐?是用的什么编辑器吗？已帮忙看完程序，未发现问题所在

就是用的txt编辑器呀，就是习惯比较好

蜻蜓队长 · 发表于 2017-12-25 10:17:38

shirlyqiu 发表于 2017-12-25 10:13
就是用的txt编辑器呀，就是习惯比较好

打call，太棒了这个习惯

shirlyqiu · 发表于 2017-12-25 10:48:34

大神们是否能来帮忙看看是什么问题，实在找不到原因了

shirlyqiu · 发表于 2017-12-25 20:48:44

最后找到的问题是不小心在批处理的时候把ctl文件的dset ^....ctl，但是第一次检查的没有发现，所以悲剧了

所以以上的代码和描述文件都是没有问题的。。。

		自动登录	找回密码
密码			立即注册

[脚本编辑] 【求助】多时次站点资料作图生产映射文件发生错误

登录后查看更多精彩内容~