也发一个格点插值为站点数据的程序

MeteoInfo · 发表于 2013-7-23 16:52:06

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

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

x

本帖最后由 MeteoInfo 于 2013-7-31 13:52 编辑

之前已经有一些帖子讨论了格点插值为站点数据的算法了，通常会用双线性插值，也有用最近距离方法。MeteoInfo C# 和 Java版都有双线性插值的算法（其实挺简单的），因为工作需要最近用Fortran写了一个格点插值为站点数据的程序，希望对大家有用。程序支持两种插值方法：双线性和最近距离。

格点转站点的fortran程序见附件，编译后（比如a.out）用法：a.out -i -o -s -x -m

-i后是输入文件名（micaps 4格点数据文件，必须有），-o后面是输出文件名（micaps 3站点数据文件，必须有）， -s后面是站点列表数据（stations.txt，必须有，每一行包含站点号、纬度、经度、高度，用空格分隔），-x后面是n（将最近格点的值赋给站点）或者i（双线性插值，用站点周边4个格点数据进行双线性插值获得站点值）（该参数是可选的，缺省是双线性插值），-m后面是最大站点数（该参数是可选的，缺省为500）。

例如：a.out -i13071408.003 -otest.003 -sstations.txt
a.out -i13071408.003 -otest_n.003 -sstations.txt -xn

这个程序做得比较灵活，可以根据需要进行使用。如果输入和输出数据的格式不同可以修改输入输出部分的代码。
7月31日更新了程序，能够处理不同文件头写法的micaps 4数据。

tostation.F (5.55 KB, 下载次数: 355)

2013-7-31 13:50 上传

点击文件名下载附件

站点列表文件内容示例如下：
50136       53.47       122.37       297.3
50353       51.72       126.65       178.2
50434       50.48       121.68       733.4
50442       50.40       124.12       372.1
50468       50.25       127.45       166.9
50514       49.57       117.43       662.0
50527       49.22       119.75       612.5
50557       49.17       125.23       243.0
50603       48.67       116.82       555.6

代码也贴出来吧：

PROGRAM GRIDTOSTATION
!-----------------------------------------------
! Author: Yaqiang Wang
! Date: Jul. 15, 2013
! Purpose: Interpolate grid data to station data
! Method: Neighbor or Bi-linear interpolation
! Note: Input is MICAPS 4 data; Output is MICAPS 3 data;
! Station data file is space delimited
! (stid latitude longitude altitude)
!-------------------------------------------------
integer :: xnum,ynum,error,nstn,mstn,xidx,yidx
integer :: iyear,imonth,iday,ihour,phour,level
integer, allocatable :: sid(:)
real :: mvalue,x1val,x2val,sum,a,b,c,d
real :: xmin,xmax,ymin,ymax,xdelt,ydelt,xi,yi
real, allocatable :: xarray(:),yarray(:),data(:,:)
real, allocatable :: conc(:),slat(:),slon(:),salt(:)
character(80) :: finp,fout,fstn
character(200) :: ctemp,desc
character(500):: label
logical :: ftest,linear
!--------------------------------------------------
! Command line argument list
!--------------------------------------------------
karg=iargc()
if(karg.eq.0)then
write(*,*)'USAGE: tostation -i -o -s -x -m]'
write(*,*)' -i[Input file name (MICAPS 4 data)]'
write(*,*)' -o[Output file name (MICAPS 3 data)]'
write(*,*)' -s[Station list file name]'
write(*,*)' -x[n(neighbor) or i(interpolation)]'
write(*,*)' -m[Maximum number of stations (500)]'
write(*,*)' '
write(*,*)'NOTE: leave no space between option and value'
stop
end if
! Set default values
mvalue=9999.0
nstn=0
mstn=500
linear=.true.
! Go through each argument
do while (karg.gt.0)
call getarg(karg,label)
select case (label(1:2))
case ('-i','-I') ! Input grid data file
klen=index(label,' ')-1
finp=label(3:klen)
case ('-o','-O') ! Output station data file
klen=index(label,' ')-1
fout=label(3:klen)
case ('-s','-S') ! Station lat/lon file
klen=index(label,' ')-1
fstn=label(3:klen)
case ('-m','-M') ! Maximum station number
klen=index(label,' ')-1
read(label(3:),'(I10)') mstn
case ('-x','-X') ! Interpolation methods
klen=index(label,' ')-1
ctemp=label(3:klen)
if(ctemp(1:1).eq.'n'.or.ctemp(1:1).eq.'N')linear=.false.
end select
karg=karg-1
end do
write(*,*)linear
allocate (conc(mstn),sid(mstn),slat(mstn),slon(mstn),salt(mstn))
!----------------------------------------------------
! Open station file and read station id and lon/lat
!----------------------------------------------------
inquire (file=fstn,exist=ftest)
if(ftest)then
open(15,file=fstn)
else
write(*,*)'Station list file not found: ',fstn
stop
end if
do while(nstn.lt.mstn)
nstn=nstn+1
read(15,*,end=100)sid(nstn),slat(nstn),slon(nstn),salt(nstn)
! write(*,*)sid(nstn),slat(nstn),slon(nstn)
end do
100 close(15)
nstn=nstn-1
!----------------------------------------------------
! Open data file and allocate variables
!----------------------------------------------------
inquire (file=finp,exist=ftest)
if(ftest)then
open(10,file=finp,status='old',iostat=error)
else
write(*,*)'Input file not found: ',finp
stop
end if
read(10,*) ctemp,ctemp,desc
write(*,'(A)') desc
read(10,*) iyear,imonth,iday,ihour,phour,level,xdelt,ydelt,xmin,xmax,ymin,ymax,xnum,ynum&
,ctemp,ctemp,ctemp,ctemp,ctemp
write(*,*) iyear,imonth,iday,ihour,phour,level,xdelt,ydelt,xmin,xmax,ymin,ymax,xnum,ynum
allocate (xarray(xnum),yarray(ynum),data(ynum,xnum))
do i=1,xnum
xarray(i)=xmin+xdelt*(i-1)
end do
do i=1,ynum
yarray(i)=ymin+ydelt*(i-1)
end do
read(10,*)((data(j,i),i=1,xnum),j=1,ynum)
close(10)
!--------------------------------------------------------------------------
! Interpolation
!--------------------------------------------------------------------------
do i=1,nstn
conc(i)=mvalue
xi = (slon(i)-xarray(1))/xdelt+1.0
yi = (slat(i)-yarray(1))/ydelt+1.0
xidx = INT(xi)
yidx = INT(yi)
! write(*,*)xi,yi
if(xidx.ge.1.and.xidx.lt.xnum.and.yidx.ge.1.and.ydix.lt.ynum)then
if(linear)then
i1=INT(yi)
j1=INT(xi)
i2=i1+1
j2=j1+1
a=data(i1,j1)
b=data(i1,j2)
c=data(i2,j1)
d=data(i2,j2)
! write(*,'(I3,I4,I4,I4,F8.2)')i1,i2,j1,j2,data(i1,j1)
n=0
sum=0.0
if(a.ne.mvalue)then
sum=sum+a
n=n+1
end if
if(b.ne.mvalue)then
sum=sum+b
n=n+1
end if
if(c.ne.mvalue)then
sum=sum+c
n=n+1
end if
if(d.ne.mvalue)then
sum=sum+d
n=n+1
end if
if(n.eq.0)then
cycle
else if(n.le.3)then
conc(i)=sum/n
else
x1val=a+(c-a)*(slat(i)-yarray(i1))/ydelt
x2val=b+(d-b)*(slat(i)-yarray(i1))/ydelt
conc(i)=x1val+(x2val-x1val)*(slon(i)-xarray(j1))/xdelt
write(*,'(I3,F8.2,F8.2,F8.2,F8.2,F8.2,F8.2,F8.2)')i,a,b,c,d,x1val,x2val,conc(i)
end if
else
i1=NINT(yi)
j1=NINT(xi)
conc(i)=data(i1,j1)
write(*,'(I3,I4,I4,F8.2)')i,i1,j1,conc(i)
end if
end if
end do
!----------------------------------------------------------------------------------
! Write output station file
!----------------------------------------------------------------------------------
open(20,file=fout)
ctemp = 'diamond 3 ' // desc
write(20,'(A80)') ctemp
ctemp=' 0 0 0 0 1 '
write(20,'(I4,I3,I3,I3,I5,A11,I3)') iyear,imonth,iday,ihour,level,ctemp,nstn
do i=1, nstn
write(20,'(I5,F8.2,F7.2,F8.2,F8.2)') sid(i),slon(i),slat(i),salt(i),conc(i)
end do
close(20)
stop
end

PROGRAM GRIDTOSTATION 
!----------------------------------------------- 
! Author: Yaqiang Wang 
! Date: Jul. 15, 2013 
! Purpose: Interpolate grid data to station data 
! Method: Neighbor or Bi-linear interpolation 
! Note: Input is MICAPS 4 data; Output is MICAPS 3 data; 
!       Station data file is space delimited 
!       (stid  latitude  longitude altitude) 
!------------------------------------------------- 
 
integer       :: xnum,ynum,error,nstn,mstn,xidx,yidx 
integer       :: iyear,imonth,iday,ihour,phour,level 
integer,      allocatable :: sid(:) 
real          :: mvalue,x1val,x2val,sum,a,b,c,d 
real          :: xmin,xmax,ymin,ymax,xdelt,ydelt,xi,yi 
real,         allocatable :: xarray(:),yarray(:),data(:,:) 
real,         allocatable :: conc(:),slat(:),slon(:),salt(:) 
character(80) :: finp,fout,fstn 
character(200) :: ctemp,desc 
character(500):: label 
logical       :: ftest,linear 
 
!-------------------------------------------------- 
! Command line argument list 
!-------------------------------------------------- 
 
karg=iargc() 
 
if(karg.eq.0)then 
  write(*,*)'USAGE: tostation -i -o -s -x -m]' 
  write(*,*)'  -i[Input file name (MICAPS 4 data)]' 
  write(*,*)'  -o[Output file name (MICAPS 3 data)]' 
  write(*,*)'  -s[Station list file name]' 
  write(*,*)'  -x[n(neighbor) or i(interpolation)]' 
  write(*,*)'  -m[Maximum number of stations (500)]' 
  write(*,*)'  ' 
  write(*,*)'NOTE: leave no space between option and value' 
  stop 
end if 
 
! Set default values 
mvalue=9999.0 
nstn=0 
mstn=500 
linear=.true. 
 
! Go through each argument 
do while (karg.gt.0) 
  call getarg(karg,label) 
  select case (label(1:2)) 
    case ('-i','-I')    ! Input grid data file 
      klen=index(label,' ')-1 
      finp=label(3:klen) 
    case ('-o','-O')    ! Output station data file 
      klen=index(label,' ')-1 
      fout=label(3:klen) 
    case ('-s','-S')    ! Station lat/lon file 
      klen=index(label,' ')-1 
      fstn=label(3:klen) 
    case ('-m','-M')    ! Maximum station number 
      klen=index(label,' ')-1 
      read(label(3:),'(I10)') mstn 
    case ('-x','-X')    ! Interpolation methods 
      klen=index(label,' ')-1 
      ctemp=label(3:klen) 
      if(ctemp(1:1).eq.'n'.or.ctemp(1:1).eq.'N')linear=.false.       
  end select 
  karg=karg-1 
end do 
write(*,*)linear 
allocate (conc(mstn),sid(mstn),slat(mstn),slon(mstn),salt(mstn)) 
 
!---------------------------------------------------- 
! Open station file and read station id and lon/lat 
!---------------------------------------------------- 
 
inquire (file=fstn,exist=ftest) 
if(ftest)then 
  open(15,file=fstn) 
else 
  write(*,*)'Station list file not found: ',fstn 
  stop 
end if 
 
do while(nstn.lt.mstn) 
  nstn=nstn+1 
  read(15,*,end=100)sid(nstn),slat(nstn),slon(nstn),salt(nstn) 
!  write(*,*)sid(nstn),slat(nstn),slon(nstn) 
end do 
100 close(15) 
nstn=nstn-1 
 
!---------------------------------------------------- 
! Open data file and allocate variables 
!---------------------------------------------------- 
 
inquire (file=finp,exist=ftest) 
if(ftest)then 
  open(10,file=finp,status='old',iostat=error) 
else 
  write(*,*)'Input file not found: ',finp 
  stop 
end if 
 
read(10,*) ctemp,ctemp,desc 
write(*,'(A)') desc 
read(10,*) iyear,imonth,iday,ihour,phour,level,xdelt,ydelt,xmin,xmax,ymin,ymax,xnum,ynum& 
  ,ctemp,ctemp,ctemp,ctemp,ctemp 
write(*,*) iyear,imonth,iday,ihour,phour,level,xdelt,ydelt,xmin,xmax,ymin,ymax,xnum,ynum 
allocate (xarray(xnum),yarray(ynum),data(ynum,xnum)) 
do i=1,xnum 
  xarray(i)=xmin+xdelt*(i-1) 
end do 
do i=1,ynum 
  yarray(i)=ymin+ydelt*(i-1) 
end do 
read(10,*)((data(j,i),i=1,xnum),j=1,ynum) 
close(10) 
 
!-------------------------------------------------------------------------- 
! Interpolation 
!-------------------------------------------------------------------------- 
 
do i=1,nstn 
  conc(i)=mvalue 
  xi = (slon(i)-xarray(1))/xdelt+1.0   
  yi = (slat(i)-yarray(1))/ydelt+1.0 
  xidx = INT(xi) 
  yidx = INT(yi) 
!  write(*,*)xi,yi 
  if(xidx.ge.1.and.xidx.lt.xnum.and.yidx.ge.1.and.ydix.lt.ynum)then 
    if(linear)then 
      i1=INT(yi) 
      j1=INT(xi) 
      i2=i1+1 
      j2=j1+1 
      a=data(i1,j1) 
      b=data(i1,j2) 
      c=data(i2,j1) 
      d=data(i2,j2) 
!      write(*,'(I3,I4,I4,I4,F8.2)')i1,i2,j1,j2,data(i1,j1) 
      n=0    
      sum=0.0 
      if(a.ne.mvalue)then 
        sum=sum+a 
        n=n+1 
      end if 
      if(b.ne.mvalue)then 
        sum=sum+b 
        n=n+1 
      end if 
      if(c.ne.mvalue)then 
        sum=sum+c 
        n=n+1 
      end if 
      if(d.ne.mvalue)then 
        sum=sum+d 
        n=n+1 
      end if 
      if(n.eq.0)then 
        cycle 
      else if(n.le.3)then 
        conc(i)=sum/n 
      else 
        x1val=a+(c-a)*(slat(i)-yarray(i1))/ydelt 
        x2val=b+(d-b)*(slat(i)-yarray(i1))/ydelt 
        conc(i)=x1val+(x2val-x1val)*(slon(i)-xarray(j1))/xdelt 
        write(*,'(I3,F8.2,F8.2,F8.2,F8.2,F8.2,F8.2,F8.2)')i,a,b,c,d,x1val,x2val,conc(i) 
      end if 
    else 
      i1=NINT(yi) 
      j1=NINT(xi) 
      conc(i)=data(i1,j1) 
      write(*,'(I3,I4,I4,F8.2)')i,i1,j1,conc(i) 
    end if 
  end if 
end do 
 
!---------------------------------------------------------------------------------- 
! Write output station file 
!---------------------------------------------------------------------------------- 
open(20,file=fout) 
ctemp = 'diamond 3 ' // desc 
write(20,'(A80)') ctemp 
ctemp=' 0 0 0 0 1 ' 
write(20,'(I4,I3,I3,I3,I5,A11,I3)') iyear,imonth,iday,ihour,level,ctemp,nstn 
do i=1, nstn 
  write(20,'(I5,F8.2,F7.2,F8.2,F8.2)') sid(i),slon(i),slat(i),salt(i),conc(i) 
end do 
close(20) 
 
stop 
end

О.o贰沫 · 发表于 2013-7-23 16:59:18

谢谢王老师分享！！

又是那隻貓 · 发表于 2013-7-23 16:59:28

太棒了，！！

lqouc · 发表于 2013-7-23 18:28:50

感谢王老师分享

topmad · 发表于 2013-7-23 19:25:00

看王老师的代码能学到好多东西的说

ps 王老师你究竟会多少种语言....

java C# fortran ironPython grads ......

做个霸气的木头 · 发表于 2013-7-23 21:03:45

王老师一直是我的偶像~下下来学习学习

MeteoInfo · 发表于 2013-7-23 21:33:28

topmad 发表于 2013-7-23 19:25

登录/注册后可看大图

看王老师的代码能学到好多东西的说

ps 王老师你究竟会多少种语言....

用过的语言比较多了，还有VB, VB.Net, C++, C, Groovy, Avenue，其实都只是知道皮毛而已，只是在用的时候会学习一点，毕竟编程语言很多东西是想通的。精通一门语言更重要，不过每一门语言都博大精深，我们非软件开发为职业的人来说很难真正做到精通，基本够用就行了。

MeteoInfo · 发表于 2013-7-23 21:43:15

发现一个bug，

allocate (conc(mstn),sid(mstn),slat(mstn),slon(mstn),salt(mstn))

应该放在从参数中读取mstn变量之后，已经在一楼中进行了更新。

張顧.煒Zane · 发表于 2013-7-23 21:44:14

收藏了！！

ぁ無聲か · 发表于 2013-7-23 23:15:32

谢谢王老师分享！！

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

[源代码] 也发一个格点插值为站点数据的程序

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

评分

点评