file:///C:/Users/hp/AppData/Local/Temp/msohtml1/01/clip_image002.gif
首先感谢武大毕业的尹雄锐博士提供的程序.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 本程序用于2001年气象数据空间插值
c 2006-9-28
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
parameter(n_sta=35,n=365,nn=1,n_pix=584915)
c n－－－一年的总天数 nn－－－需要插值的天数
dimension elev_sta(n_sta),temp_sta(n,n_sta)
real mxtemp_sta(n,n_sta),mntemp_sta(n,n_sta)
dimension vpre_sta(n,n_sta),suntime_sta(n,n_sta)
dimension wind_sta(n,n_sta),p_sta(n,n_sta)
character*50 No_sta(n_sta),Name_sta(n_sta),Name_prov,temps
c n_sta----------------------气象站站点数
c NO_sta---------------------气象站号
c elev_sta-------------------气象站高程
c temp_sta-------------------气象站气温
c mxtemp_sta,mntemp_sta------气象站最高、最低气温
c vpre_sta-------------------气象站水气压
c suntime_sta----------------气象站日照时数
c wind_sta-------------------气象站风速
c P_sta----------------------气象站降雨
c Name_sta,Name_prov---------气象站站名和所在省名
dimension cor_x(n_sta),cor_y(n_sta)
integer year(n,n_sta),month(n,n_sta),day(n,n_sta)
c cor_x,cor_y ---------气象站经纬度坐标
c
dimension x(n_pix),y(n_pix),dem(n_pix)
c x，y---------网格投影坐标（单位：米）
c dem----------网格高程
integer ye(nn),mo(nn),da(nn)
c 插值日期
dimension x0(n_sta),y0(n_sta)
c x0，y0---------气象站投影坐标(Alberes,WGS84)
real temperature(n_sta),mxtemp(n_sta),mntemp(n_sta)
dimension vpre(n_sta),suntime(n_sta),wind(n_sta),p(n_sta)
dimension Vpre_out(n_pix),wind_out(n_pix),suntime_out(n_pix)
dimension temperature_out(n_pix),P_out(n_pix)
real lanbuda,Ning
c 读入气象站点各项数据
open(1,file='2001.txt',status='old')
read(1,*)
do i=1,n_sta
do j=1,n
read(1,*)NO_sta(i),cor_x(i),cor_y(i),Name_sta(i),Name_prov,temp1,
$ elev_sta(i),temp2,temps,temp3,temp4,temp5,temp6,temp7,temp8,
$ year(j,i),month(j,i),day(j,i),temp9,temp_sta(j,i),
$ mxtemp_sta(j,i),mntemp_sta(j,i),
$vpre_sta(j,i),temp10,temp11,wind_sta(j,i),p_sta(j,i),
$ temp12,suntime_sta(j,i)
enddo
enddo
close(1)
c do i=1,n_sta
c write(*,*)cor_x(i),cor_y(i)
c enddo
c 读入栅格坐标
open(2,file='coordinate.txt',status='old')
do i=1,n_pix
read(2,*) x(i),y(i),dem(i)
enddo
close(2)
c 输入日期
open(3,file='input_data.txt',status='old')
do i=1,nn
read(3,*)ye(i),mo(i),da(i)
enddo
close(3)
c open(211,file='wind.txt')
c open(333,file='relation.txt')
c 找出该天的所有气象站点的数据并进行插值
do ii=1,1 !nn
!挑出当天各站的气象数据
do i=1,n_sta
do j=1,n
if(year(j,i).eq.ye(ii).and.month(j,i).eq.mo(ii).
$ and.day(j,i).eq.da(ii))then
temperature(i)=temp_sta(j,i)*0.1 !0.1--scale factor
mxtemp(i)=mxtemp_sta(j,i)*0.1
mntemp(i)=mntemp_sta(j,i)*0.1
vpre(i)=vpre_sta(j,i)*0.1
suntime(i)=suntime_sta(j,i)*0.1
wind(i)=wind_sta(j,i)*0.1
c 降雨对前五天累计进行插值,要保证j>5，否则不行
p(i)=0.0
do kk=1,5
if(P_sta(j-kk,i).gt.3000) P_sta(j-kk,i)=0.0
p(i)=p(i)+P_sta(j-kk,i)*0.1
enddo
endif
enddo
write(*,*)temperature(i),vpre(i),suntime(i),P(i)
c write(333,*)elev_sta(i),suntime(i)*0.1
! 找出各气象站点的风速，不插值
lanbuda=cor_x(i)
fi=cor_y(i)
write(*,*)lanbuda,fi,wind(i)
call proj(lanbuda,fi,Eing,Ning)
x0(i)=Eing
y0(i)=Ning
c write(211,212)lanbuda*180/3.141593,fi*180/3.141593,x0(i)
c $,y0(i),wind(i),name_sta(i)
c212 format(2(1x,f10.4),2(1x,f10.1),3x,f4.1,2x,a20)
enddo
c close(211)
do i=1,n_pix
wind_out(i)=-10.0
do j=1,n_sta
if(abs(x(i)-x0(j)).le.500.and.abs(y(i)-y0(j)).le.500) then
wind_out(i)=wind(j)
write(*,*) 'find',j,wind_out(i)
endif
enddo
enddo
!------------------------------------------------------------------
! 对网格点插值
! 水汽压插值，对高程有明显的趋势，随高程指数下降
do i=1,n_sta
vpre(i)=vpre(i)*10**(elev_sta(i)/5000.0)
enddo
callDDWA(x0,y0,vpre,n_sta,x,y,n_pix,vpre_out)
do i=1,n_pix
if(dem(i).eq.-32768)then
vpre_out(i)=-999
else
vpre_out(i)=vpre_out(i)*10**(-dem(i)/5000.0)
endif
enddo
!------------------------------------------------------------------
! 风速插值，对比高程和风速没有明显的趋势
c call DDWA(x0,y0,wind,n_sta,x,y,n_pix,wind_out)
!-------------------------------------------------------------------
! 太阳日照时数插值, 对高程没有趋势
call DDWA(x0,y0,suntime,n_sta,x,y,n_pix,suntime_out)
!-------------------------------------------------------------------
! 气温插值，需要将气象站高度订正到同一平面再插值，然后网格再算回去
do i=1,n_sta
offset=elev_sta(i)/100.0*0.6 !海拔每上升100米，气温平均下降0.6度
temperature(i)=temperature(i)+offset
enddo
callDDWA(x0,y0,temperature,n_sta,x,y,n_pix,temperature_out)
do i=1,n_pix
if(dem(i).eq.-32768)then
temperature_out(i)=-999
else
offset=dem(i)/100*0.6
temperature_out(i)=temperature_out(i)-offset
endif
enddo
c
!------------------------------------------------------------------------
! 降雨插值
call DDWA(x0,y0,P,n_sta,x,y,n_pix,P_out)
do i=1,n_pix
if(dem(i).eq.-32768)then
p_out(i)=-999
else
p_out(i)=p_out(i)
endif
enddo
!------------------------------------------------------------------------
open(4,file='Vpre.txt')
do i=1,n_pix
write(4,111)x(i),y(i),Vpre_out(i)
111 format(1x,f10.1,1x,f10.1,1x,f10.3)
enddo
close(4)
open(5,file='Wind.txt')
do i=1,n_pix
write(5,112)x(i),y(i),wind_out(i)
112 format(1x,f10.1,1x,f10.1,1x,f10.3)
enddo
close(5)
open(6,file='suntime.txt')
do i=1,n_pix
if(dem(i).eq.-32768) then
suntime_out(i)=-999
else
suntime_out(i)=suntime_out(i)
endif
write(6,113)x(i),y(i),suntime_out(i)
113 format(1x,f10.1,1x,f10.1,1x,f10.3)
enddo
close(6)
open(7,file='Temperature.txt')
do i=1,n_pix
write(7,114)x(i),y(i),Temperature_out(i)
114 format(1x,f10.1,1x,f10.1,1x,f10.3)
enddo
close(7)
open(8,file='precipitation.txt')
do i=1,n_pix
write(8,115)x(i),y(i),P_out(i)
115 format(1x,f10.1,1x,f10.1,1x,f10.3)
enddo
close(8)
enddo
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 该子程序用于将点的经纬度坐标转化为Albers的投影坐标,用于将气象站的经纬度转化为米
c 经纬度为度的形式，非度分秒形式
c 取中国区的参数为：fi1-------第一标准纬度25
c fi2-------第二标准纬度47
c lanbuda0--中央经线105
c fi0-------中央纬线0
c EF--------东偏0,单位：米
c NF--------北偏0,单位：米
c WGS84椭球体参数 长半轴a=6378137
c 短半轴b=6356752.3142
c Karasovsky椭球体参数 长半轴a=6378245
c 短半轴b=6356863.0188
c 偏心率e=SQRT(2f-f^2),1/f=a/(a-b)
cc 这里采用WGS84
subroutineproj(lanbuda,fi,Eing,Ning)
c 输入（地理经纬度）lanbuda------精度
c fi-----------纬度
c 输出（投影坐标m）Eing---------东的坐标
c Ning---------北的坐标
real lanbuda,lanbuda0,NF,m1,m2,n,Ning
pi=4*atan(1.0)
fi0=0*pi/180
lanbuda0=105*pi/180
fi1=25*pi/180
fi2=47*pi/180
EF=0
NF=0
fi=fi*pi/180
lanbuda=lanbuda*pi/180
a=6378137
b=6356752.3142
e1=SQRT(1-b/a*b/a)
e2=1-e1*e1
alfa=e2*(sin(fi)/(1-e1*e1*sin(fi)*sin(fi))-0.5/e1
$ *log((1-e1*sin(fi))/(1+e1*sin(fi))))
alfa0=e2*(sin(fi0)/(1-e1*e1*sin(fi0)*sin(fi0))-0.5/e1
$ *log((1-e1*sin(fi0))/(1+e1*sin(fi0))))
alfa1=e2*(sin(fi1)/(1-e1*e1*sin(fi1)*sin(fi1))
$ -0.5/e1*log((1-e1*sin(fi1))/(1+e1*sin(fi1))))
alfa2=e2*(sin(fi2)/(1-e1*e1*sin(fi2)*sin(fi2))
$ -0.5/e1*log((1-e1*sin(fi2))/(1+e1*sin(fi2))))
m1=cos(fi1)/sqrt(1-e1*e1*sin(fi1)*sin(fi1))
m2=cos(fi2)/sqrt(1-e1*e1*sin(fi2)*sin(fi2))
n=(m1*m1-m2*m2)/(alfa2-alfa1)
C=m1*m1+n*alfa1
rou0=a*sqrt(C-n*alfa0)/n
rou=a*sqrt(C-n*alfa)/n
sita=n*(lanbuda-lanbuda0)
Eing=EF+rou*sin(sita)
Ning=NF+rou0-rou*cos(sita)
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 该子程序采用距离方向加权平均法进行空间插值
c x0,y0---------气象站点坐标,n0-----站数（站号）
c x1,y1---------待求点坐标
c P-------------气象站实测值
c P-out---------点插值结果
c dis(j)--------与第j个站的距离
cc
subroutine DDWA(x0,y0,p,n0,x1,y1,n1,P_out)
dimension x0(n0),y0(n0),x1(n1),y1(n1)
dimension P(n0),P_out(n1)
dimension dis(n0),num(n0),w(n0),alf(n0)
c 通过最近8个站计算每个目标点的值
do i=1,n1
do j=1,n0
dis(j)=sqrt((x1(i)-x0(j))*(x1(i)-x0(j))
$ +(y1(i)-y0(j))*(y1(i)-y0(j)))
if(dis(j).eq.0) then
P_out(i)=p(j)
goto 33
endif
enddo
!找出最近的8个站及其与待求点的距离
do j=1,n0
num(j)=j
enddo
do ii=1,8
s=1000000000
do j=ii+1,n0
if(dis(j).lt.s) then
s=dis(j)
k=j
endif
enddo
if(s.lt.dis(ii))then
dis(k)=dis(ii)
dis(ii)=s
tempk=num(ii)
num(ii)=num(k)
num(k)=tempk
endif
enddo
!计算权重距离
dis0=dis(8) !取8个最远的点的距离（经验衰减距离）
do j=1,8
w(j)=(exp(-dis(j)/dis0))**4.0
enddo
!计算方向修正系数
do j=1,8
s1=0
s2=0
do l=1,8
s1=s1+w(l)
d=SQRT((x0(num(j))-x0(num(l)))**2
$ +(y0(num(j))-y0(num(l)))**2)
cos=(dis(j)*dis(j)+dis(l)*dis(l)-d*d)/(2*dis(j)*dis(l))
s2=s2+w(l)-w(l)*cos
enddo
alf(j)=s2/s1
enddo
!计算总修正权重
do j=1,8
w(j)=w(j)*(1+alf(j))
enddo
!计算目标点值
s1=0
s2=0
do j=1,8
s1=s1+w(j)*P(num(j))
s2=s2+w(j)
enddo
P_out(i)=s1/s2
33 enddo
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 该子程序采用修正距离倒数平方法
c x0,y0---------气象站点坐标,n0-----站数（站号）
c x1,y1---------待求点坐标
c dem-----------待求点高程
c P-------------气象站实测值
c P-out---------点插值结果
c dis(j)--------与第j个站的距离
cc
subroutine MIDW(x0,y0,elev_sta,p,n0,x1,y1,dem,n1,P_out)
dimension x0(n0),y0(n0),elev_sta(n0),x1(n1),y1(n1),dem(n1)
dimension P(n0),P_out(n1)
dimension dis(n0),alf(n0)
b=-2 ! b是权重指数，一般情况下取负值
do i=1,n1
ss=0.0
do j=1,n0
de=abs(dem(i)-elev_sta(j))
if(de.eq.0) de=1
c de=1 !如果de＝1，则变为距离倒数平方法
dis(j)=SQRT((x1(i)-x0(j))**2+(y1(i)-y0(j))**2)
if(dis(j).eq.0)then
P_out(i)=p(j)
goto 440
endif
dis(j)=(dis(j)/de)**b
ss=ss+dis(j)
enddo
do j=1,n0
alf(j)=dis(j)/ss
write(*,*)alf(j)
p_out(i)=p_out(i)+alf(j)*p(j)
enddo
440 enddo
return
end
窥视卡
雷达卡
发表于 2013-4-7 15:11:51
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
发表于 2013-4-7 15:14:07
楼主