- 积分
- 1274
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-4-29
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
此程序为网格点数据插值到站点,一个时次。
请大神指点在哪里加循环,改为4个时次的循环
- ! 双线性插值,格点到站点,考虑权重 !
- ! modifioed from mofangbao(清风) by djj!
- !===================================================!
- program bilineargird2stn
- integer,parameter::oldx=97,oldy=61 !格点数据xy向网格点数
- integer,parameter::newxy=1308 !站点数
- integer,parameter::time=4 !时间层次
- real,parameter::oldLon_start=100.0,oldLat_start=20.0,oldInterpx=0.25,oldInterpy=0.25 !格点数据起始经纬度和格距
- 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='C:\PCGrADS\win32\st1308-2.txt')
- do i=1,newxy
- read(27,*) stid,newPoint(i,2),newPoint(i,1) !先经度后纬度
- enddo
- close(27)
- !!读取网格点数据
- open(29,file='C:\PCGrADS\win32\24cma2015081619.grib',form='binary')
- read(29) ((oldVal(i,j),i=1,oldx),j=1,oldy)
- close(29)
- !do k=1,4
- !!调用插值程序,这里只哟一个时次,因此没有循环
- call getNewVal(oldVal,newVal,oldLon_start,oldLat_start,oldInterpx,oldInterpy,newPoint,isFirst)
- open(25,file='C:\PCGrADS\win32\grid2stn.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=97,oldy=61,newxy=1308
- 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,1308
- 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(25,file='C:\PCGrADS\win32\weight.txt')
- do i=1,1308
- write(25,*) weight(i,:)
- enddo
- close(25)
- 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
复制代码
|
|