- 积分
- 2729
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2015-10-7
- 最后登录
- 1970-1-1
![未绑定新浪微博用户 新浪微博达人勋](source/plugin/sina_login/img/gray.png)
|
![](static/image/common//ico_lz.png)
楼主 |
发表于 2016-7-21 19:10:10
|
显示全部楼层
! 双线性插值,格点到站点,考虑权重 !
! modifioed from mofangbao(清风) by djj!
!===================================================!
program bilineargird2stn
integer,parameter::oldx=29,oldy=17 !格点数据xy向网格点数
integer,parameter::newxy=192 !站点数
integer,parameter::time=1 !时间层次
real,parameter::oldLon_start=70.0,oldLat_start=15.0,oldInterpx=2.5,oldInterpy=2.5 !格点数据起始经纬度和格距
real oldVal(oldx,oldy),newVal(newxy),newPoint(newxy,2) !分别是:格点数据,站点数据,站点数据经纬度,newVal的点顺序和nowPoint中的对应
character*8 ::stid
integer::i,j,k
logical isFirst !判断是否第一个时次
isFirst=.true. !后面无需手动修改此值,程序会自动修改
!!读取站点经纬度
open(27,file='E:\data\ob\194yuzhi\lon_lat.txt')
do i=1,newxy
read(27,*) stid,newPoint(i,2),newPoint(i,1) !先经度后纬度
enddo
close(27)
!!读取网格点数据
open(29,file='E:\data\NCEP-1\1\yuzhi\nc1tmin2.5.grd',form='binary')
read(29) ((oldVal(i,j),i=1,oldx),j=1,oldy)
close(29)
!!调用插值程序,这里只有一个时次,因此没有循环
call getNewVal(oldVal,newVal,oldLon_start,oldLat_start,oldInterpx,oldInterpy,newPoint,isFirst)
open(25,file='nc1grid2stn.txt')
do i=1,newxy
write(25,*) newVal(i)
enddo
end
subroutine getNewVal(oldVal,newVal,oldLon_start,oldLat_start,oldInterpx,oldInterpy,newPoint,isFirst)
integer,parameter::oldx=29,oldy=17,newxy=192
real::oldLon_strat,oldLat_strat
real,parameter::undef=-999.0
real::oldInterpx,oldInterpy
integer pos(newxy,4),nstart,nend !weight 1,2是左侧点和右侧点的分别权重,3,4是下面点和上面点分别的权重
real tmpVal_up,tmpVal_down !本次纬向下方和纬向上方的线性插值结果以及权重
real::weight(newxy,4)
real newLon,newLat,weight_end !本次插值的新点位置,以及返回的权重
real oldVal(oldx,oldy),newVal(newxy),newPoint(newxy,2) !旧值,新值,新点经纬度
logical isfirst !是否第一个时次,只有第一个时次需要计算四周格点分布以及权重,后面的套用即可
if(isfirst)then
!获取新点的四周的点以及每个点的权重
do i=1,newxy
newLon=newPoint(i,1)
! print*,newLon,oldLon_start
call getPosAndWeight(newLon,oldLon_start,nstart,nend,weight_end,oldInterpx,isFirst)
!print*,newLon,nstart,nend,weight_end
if(nend>oldx)then
if(nstart==oldx)then
pos(i,1)=nstart
pos(i,2)=nstart
weight(i,1)=1
weight(i,2)=0
else
pos(i,1)=-1
pos(i,2)=-1 !超出了原来格点的范围,标记为无效值
weight(i,1)=-1
weight(i,2)=-1
endif
else
pos(i,1)=nstart
pos(i,2)=nend
weight(i,1)=1-weight_end
weight(i,2)=weight_end
endif
!print*, pos(i,1),pos(i,2)
!write(*,*) weight(i,1),weight(i,2)
newLat=newPoint(i,2)
! write(*,*) newLat,oldLat_start
call getPosAndWeight(newLat,oldLat_start,nstart,nend,weight_end,oldInterpy,isFirst)
if(nend>oldy)then
if(nstart==oldy)then
pos(i,3)=nstart
pos(i,4)=nstart
weight(i,3)=1
weight(i,4)=0
else
pos(i,3)=-1
pos(i,4)=-1 !超出了原来格点的范围,标记为无效值
weight(i,3)=-1
weight(i,4)=-1
endif
else
pos(i,3)=nstart
pos(i,4)=nend
weight(i,3)=1-weight_end
weight(i,4)=weight_end
endif
!print*, pos(i,3),pos(i,4)
!write(*,*) weight(i,3),weight(i,4)
!pause
enddo
endif
!开始插值
do i=1,newxy
if(pos(i,1)<0 .or.pos(i,3)<0)then
newval(i)=-999.0
else
!先做纬向插值
if(pos(i,1)==undef .or. pos(i,2)==undef .or. pos(i,3)==undef .or. pos(i,4)==undef)then
newVal(i)=-9999.0
else
if(oldVal(pos(i,1),pos(i,3))>0.and.oldVal(pos(i,2),pos(i,3))>0) tmpVal_down=weight(i,1)*oldVal(pos(i,1),pos(i,3))+weight(i,2)*oldVal(pos(i,2),pos(i,3))
if(oldVal(pos(i,1),pos(i,3))>0.and.oldVal(pos(i,2),pos(i,3))<0) tmpVal_down=oldVal(pos(i,1),pos(i,3))
if(oldVal(pos(i,1),pos(i,3))<0.and.oldVal(pos(i,2),pos(i,3))>0) tmpVal_down=oldVal(pos(i,2),pos(i,3))
if(oldVal(pos(i,1),pos(i,3))<0.and.oldVal(pos(i,2),pos(i,3))<0) tmpVal_down=-999.0
if(oldVal(pos(i,1),pos(i,4))>0.and.oldVal(pos(i,2),pos(i,4))>0) tmpVal_up=weight(i,1)*oldVal(pos(i,1),pos(i,4))+weight(i,2)*oldVal(pos(i,2),pos(i,4))
if(oldVal(pos(i,1),pos(i,4))<0.and.oldVal(pos(i,2),pos(i,4))>0) tmpVal_up=oldVal(pos(i,2),pos(i,4))
if(oldVal(pos(i,1),pos(i,4))>0.and.oldVal(pos(i,2),pos(i,4))<0) tmpVal_up=oldVal(pos(i,1),pos(i,4))
if(oldVal(pos(i,1),pos(i,4))<0.and.oldVal(pos(i,2),pos(i,4))<0) tmpVal_up=-999.0
if(tmpVal_down>0.and.tmpVal_up>0) then
newVal(i)=weight(i,3)*tmpVal_down+weight(i,4)*tmpVal_up
else if(tmpVal_down>0.and.tmpVal_up<0) then
newVal(i)=tmpVal_down
else if(tmpVal_down<0.and.tmpVal_up>0) then
newVal(i)=tmpVal_up
else
newVal(i)=-999.0
endif
endif
endif
enddo
open(30,file='weight.txt')
do i=1,newxy
write(30,*) weight(i,:)
enddo
close(30)
endsubroutine
!计算某个方向的两侧的点以及后侧点(上方点)的权重
subroutine getPosAndWeight(newPos,oldstart,nstart,nend,weight,oldInterp,isfirst)
integer nstart,nend !两侧点
real newPos,weight,oldInterp
real oldstart
logical isfirst
!离的近的点权重大,但是占的距离比重小,所以权重=1-距离比重,经纬度从0开始,但是标号从1开始,所以加1
!print*,newPos,oldstart,oldInterp
nstart=int((newPos-oldstart)/oldInterp)+1
nend=nstart+1
weight=(newPos-(nstart-1)*oldInterp-oldstart)/oldInterp ! |---|-------| 此表达式计算的是后一段的权重(前一段的距离比重),该段总长为interp
isfirst=.false.
endsubroutine
|
|