清风的双线性插值方法，为什么得到的TXT全是-999.00的缺测值呢？

jolincai

请帮忙看一下程序哪里有问题吗 ？我只改了一些参数  主要的程序 计算步骤没有改

lqouc

jolincai 发表于 2016-7-28 22:07
啊.大神..我都不好意思啦 但是又得学会...不知道怎么上传附件可以不用金钱就下载...楼下两个附件就要10金 ...

我大概看了下，是子程序写的有点问题，计算过程没错。具体程序的问题出在哪里我就不解释了，告诉你解决办法好了。
在主程序里面找到下面这个call，然后再上下分别加入我给的这两句计算，再运行就好了。

1. oldVal=oldVal+100
   
2. !!调用插值程序,这里只有一个时次，因此没有循环
   
3. call getNewVal(oldVal,newVal,oldLon_start,oldLat_start,oldInterpx,oldInterpy,newPoint,isFirst)
   
4. newVal=newVal-100

复制代码

程序出错的原因在于你的网格数据存在负值，清风的这个程序只能对正值进行差值，负值会被认为是缺省（我理解的他的程序是这样）。所以以后用到这个程序的时候就先在原来数值的基础上加一个足够大的数，让你的数据全部变成正值，插值之后再减去还原为原来的数值进行输出。
ps：在Fortran板块我下载东西是不需要金钱的，如果以后要让其他会员下载你的程序帮忙修改，直接给人家评分补偿就好了。当然也可以去免积分提问板块，那里下载附件貌似不需要金钱。

lqouc

直接把代码贴上来，讲清楚具体改了什么。

jolincai

!   双线性插值，格点到站点，考虑权重  !
!      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

jolincai

lqouc 发表于 2016-7-21 18:45
直接把代码贴上来，讲清楚具体改了什么。

不知道怎么贴代码 我就粘贴过来的，就读入了我的数据 改了格点数站点数之类的 计算方法什么的都没有改过

jolincai

lqouc 发表于 2016-7-21 18:45
直接把代码贴上来，讲清楚具体改了什么。

麻烦大神看看哪里有错误

lqouc

jolincai 发表于 2016-7-21 19:13
麻烦大神看看哪里有错误

看不出什么问题，确定站点的经纬度顺序对么？

jolincai

lqouc 发表于 2016-7-21 20:00
看不出什么问题，确定站点的经纬度顺序对么？

恩 就是我194个站的经纬度 应该是没问题的 有什么顺序的要求吗？

lqouc

jolincai 发表于 2016-7-21 20:19
恩 就是我194个站的经纬度 应该是没问题的 有什么顺序的要求吗？

显然有顺序好不好，你试试反过来。

jolincai

lqouc 发表于 2016-7-22 00:28
显然有顺序好不好，你试试反过来。

反过来是什么意思 经纬度反过来吗？还是站号的顺序？

kaiwang44

好高深，看不懂