站点eof特征向量作图问题

letring

站点资料，数据读进去且能够用set gxout print 打印在屏幕上，可就是作不出图[用到插值函数，作空间分布图]，可能是什么问题呢？求解

river

少于两个站点，无法插值啊。看看是不是你的fortran程序里面经纬度写反了，应该纬度在前

letring

river 发表于 2013-5-20 20:41

                               
登录/注册后可看大图

少于两个站点，无法插值啊。看看是不是你的fortran程序里面经纬度写反了，应该纬度在前

并没有写反，能把数据打印在屏幕上都是对的，而且个数也没有少。并且第一个时次的能话出来，第二个时次以后的就都不能画出来了，都是这种错误提示。。。。数据ctl和格点ctl的时间设置是一致的呀，怎么会这样，愁死了

river

letring 发表于 2013-5-21 11:35

                               
登录/注册后可看大图

并没有写反，能把数据打印在屏幕上都是对的，而且个数也没有少。并且第一个时次的能话出来，第二个时次以 ...

那就是你多时次资料转成grd过程种可能有问题，把程序帖上来吧。并且简单介绍下你用的资料。

蔚蓝的天

资料转成grd过程种可能有问题

nuistajun

有没有站点的EOF程序啊，跪求大神给一个

hxyj

请问能否提供一个站点EOF的程序

g18289769965

我想问一下EOF的特征向量图是用什么软件来做的？？

陶蜜云

      program main_stationEOF
        parameter(m=160,n=33,mnl=min(m,n),np=min(20,mnl),ll=m)
      parameter(ks=1,undef=32767.0)
!-----input array
      dimension f0(m,n),F(M,N)
!-----work array
      dimension gvt(ll,mnl),rgvt(ll,np),cof(mnl,n),rcof(np,n),
     &sh(n,12,m)
!-----output arrays
      dimension er(mnl,4),egvt(m,mnl),ecof(mnl,n),rain(m,n)
      dimension rer(np,2),regvt(m,np),recof(np,n)
       integer sta(m)
         real jd(m),wd(m),aaa
        character aa
c-----read data
      open(1,file='160-rain.txt',status='old')


        f0=0.0
        read(1,*)
        do I=1,M
        read(1,*) sta(i),jd(i),wd(i),(rain(i,j),j=1,n)
          do j=1,n
             f0(i,j)=rain(i,j)
          enddo
        end do
        print*,f0(1,1),f0(1,2)
!101  FORMAT(33F10.9)
        close(1)
c-----Remove the terrain or missing value.
      l=0
        do j=1,m
          if(f0(j,1).ne.undef)then
          l=l+1
            do k=1,n
            f(l,k)=f0(j,k)
            enddo
          endif
        enddo
      write(*,*)'ll=',l
c-----call the subroutine
      call reof(ll,n,mnl,np,f,ks,er,gvt,ecof,rer,rgvt,recof)
c-----Add the terrain or missing value.
      l=0
        do i=1,m
          if(f0(i,1).ne.undef)then
            l=l+1
            do k=1,mnl
              egvt(i,k)=gvt(l,k)
            end do
                do k=1,np
              regvt(i,k)=rgvt(l,k)
            end do
          else
            do k=1,mnl
              egvt(i,k)=undef
            end do
            do k=1,np
              regvt(i,k)=undef
            end do
          endif
        enddo
!-----output the result        
c-----output the error
      OPEN(10,file='er.dat')
    !  write(10,*)'error=',sum-tr
      write(10,*)'EOF lanmda (eigenvalues) from big to small'
        write(10,*)(er(i,1),i=1,mnl)
        write(*,*)'lamda='
      write(*,*)(er(i,1),i=1,mnl)
        write(10,*)'EOF accumulated eigenvalues from big to small'
        write(10,*)(er(i,2),i=1,mnl)
        write(10,*)'EOF explained variances'
        write(10,*)(er(i,3),i=1,mnl)
        write(10,*)'EOF accumulated explained variances'
        write(10,*)(er(i,4),i=1,mnl)
        write(10,*)'REOF explained variances'
        write(10,*)(rer(i,1),i=1,np)
        write(10,*)'REOF accumulated explained variances'
        write(10,*)(rer(i,2),i=1,np)
c-----output eigenvactors matrix of EOF
      OPEN(11,FILE='egvt.dat',form='binary')       
        do i=1,m
        write(11)(egvt(i,j),j=1,np)
        end do
        close(11)

      OPEN(11,FILE='egvt.txt')       
        do i=1,m
        write(11,102)sta(i),jd(i),wd(i),(egvt(i,j),j=1,np)
        enddo
102   format(i7,2x,2(f7.2,2x),20f8.2)
        close(11)
c-----output time coefficients matrix of EOF
      open(14,file='ecof.dat',form='binary')
      do i=1,n
        write(14)(ecof(j,i),j=1,np)
        end do
        close(14)

      open(14,file='ecof.txt')
      do i=1,n
        write(14,103)(ecof(j,i),j=1,np)
        end do
103        format(20f8.2)
        close(14)
c-----output loading vectors of REOF
      OPEN(21,FILE='regvt.dat',form='binary')
      do i=1,m
           write(21)(regvt(i,j),j=1,np)
        end do
        close(21)

      OPEN(21,FILE='regvt.txt')
      do i=1,m
        write(21,102)sta(i),jd(i),wd(i),(regvt(i,j),j=1,np)
      end do
!103   format(i7,2x,2(f7.2,2x),20f10.2)
        close(21)
c-----output time coefficients matrix of REOF
      OPEN(22,FILE='recof.dat',form='binary')
      do j=1,n
      write(22)(recof(i,j),i=1,np)
        end do
        close(22)

      OPEN(22,FILE='recof.txt')
      do j=1,n
      write(22,103)(recof(i,j),i=1,np)
        end do
!104   format(20f10.4)
        close(22)
      stop
        end

c-----
      subroutine reof(m,n,mnl,np,f,ks,er,egvt,ecof,rer,regvt,recof)
c-----input array
      dimension f(m,n)
c-----work arrays
      dimension cov(mnl,mnl),s(mnl,mnl),d(mnl)
c-----output arrays
      dimension er(mnl,4),egvt(m,mnl),ecof(mnl,n)
      dimension rer(np,2),regvt(m,np),recof(np,n)

c---- Preprocessing data
      call transf(m,n,f,ks)       
c---- Crossed product matrix of the data f
      call crossproduct(m,n,mnl,f,cov,sum)       
c---- Eigenvalues and eigenvectors by the Jacobi method
      eps=1e-7
      call jcb(mnl,cov,s,d,eps)
c---- Specified eigenvalues
      call arrang(m,n,mnl,d,s,er,tr)
        write(*,*)'error=',sum-tr       
c---- Normalized eignvectors and their time coefficients
      call tcoeff(m,n,mnl,f,s,er,egvt,ecof)
        write(*,*)'EOF accomplished'
c-----linear transposed (rotate)
      call rotaor(m,n,mnl,np,egvt,ecof,tr,rer,regvt,recof)
        write(*,*)'REOF accomplished'
c---- Specified eignvectors and time coefficients with
c     explained variances of REOF (if you need)
      call arrange2(m,n,np,regvt,recof,rer)

      return
      end

      subroutine transf(m,n,f,ks)
!-----input array
      dimension f(m,n)
!-----work array
      dimension fw(n)
      if(ks.eq.-1)return
      if(ks.eq.0)then                !anomaly of f
        do i=1,m
          do j=1,n
            fw(j)=f(i,j)
          enddo
          call meanvar(n,fw,af,sf,vf)
          do j=1,n
            f(i,j)=f(i,j)-af
          enddo  
        enddo
        return
      endif
      if(ks.eq.1)then                 !normalizing f
        do i=1,m
          do j=1,n
            fw(j)=f(i,j)
          enddo
          call meanvar(n,fw,af,sf,vf)
          do j=1,n
            f(i,j)=(f(i,j)-af)/sf
          enddo
        enddo
      endif
      return
      end

c-----
      subroutine crossproduct(m,n,mnl,f,cov,sum)
!-----input array
      dimension f(m,n)
!-----work array
        dimension cov(mnl,mnl)
      
      if(n-m) 10,50,50
  10  do 20 i=1,mnl
      do 20 j=i,mnl
        cov(i,j)=0.0
        do is=1,m
          cov(i,j)=cov(i,j)+f(is,i)*f(is,j)
        enddo
        cov(j,i)=cov(i,j)
  20  continue
      sum=0.
        do i=1,mnl
        sum=sum+cov(i,i)
        end do
        sum=sum/(mnl*1.)
        print*,'bb'
      return
       
  50  do 60 i=1,mnl
      do 60 j=i,mnl
        cov(i,j)=0.0
        do js=1,n
          cov(i,j)=cov(i,j)+f(i,js)*f(j,js)
          
        enddo
        cov(j,i)=cov(i,j)
  60  continue
       print*,'aa'
      sum=0.
        do i=1,mnl
        sum=sum+cov(i,i)
        end do
        sum=sum/(mnl*1.)
      return
      end

c-----
      subroutine jcb(m,a,s,d,eps)
!-----input array
      dimension a(m,m)
!-----output arrays
        dimension s(m,m),d(m)

      do 30 i=1,m
      do 30 j=1,i
        if(i-j) 20,10,20
  10    s(i,j)=1.
        go to 30
  20    s(i,j)=0.
        s(j,i)=0.
  30  continue
      g=0.
      do 40 i=2,m
        i1=i-1
        do 40 j=1,i1
  40      g=g+2.*a(i,j)*a(i,j)
      s1=sqrt(g)
      s2=eps/float(m)*s1
      s3=s1
      l=0
  50  s3=s3/float(m)
  60  do 130 iq=2,m
        iq1=iq-1
        do 130 ip=1,iq1
        if(abs(a(ip,iq)).lt.s3) goto 130
        l=1
        v1=a(ip,ip)
        v2=a(ip,iq)
        v3=a(iq,iq)
        u=0.5*(v1-v3)
        if(u.eq.0.0) g=1.
        if(abs(u).ge.1e-10) g=-sign(1.,u)*v2/sqrt(v2*v2+u*u)
        st=g/sqrt(2.*(1.+sqrt(1.-g*g)))
        ct=sqrt(1.-st*st)
        do 110 i=1,m
          g=a(i,ip)*ct-a(i,iq)*st
          a(i,iq)=a(i,ip)*st+a(i,iq)*ct
          a(i,ip)=g
          g=s(i,ip)*ct-s(i,iq)*st
          s(i,iq)=s(i,ip)*st+s(i,iq)*ct
  110     s(i,ip)=g
        do 120 i=1,m
          a(ip,i)=a(i,ip)
  120     a(iq,i)=a(i,iq)
        g=2.*v2*st*ct
        a(ip,ip)=v1*ct*ct+v3*st*st-g
        a(iq,iq)=v1*st*st+v3*ct*ct+g
        a(ip,iq)=(v1-v3)*st*ct+v2*(ct*ct-st*st)
        a(iq,ip)=a(ip,iq)
  130 continue
      if(l-1) 150,140,150
  140 l=0
      go to 60
  150 if(s3.gt.s2) goto 50
      do 160 i=1,m
        d(i)=a(i,i)
  160 continue
      return
      end

c-----
      subroutine arrang(m,n,mnl,d,s,er,tr)
!-----input arrays
      dimension d(mnl),s(mnl,mnl)
!-----output array
        dimension er(mnl,4)

        tr=0.0
      do 10 i=1,mnl
        tr=tr+d(i)/(mnl*1.)
        er(i,1)=d(i)/(mnl*1.)
  10  continue
      mnl1=mnl-1
      do 20 k1=mnl1,1,-1
      do 20 k2=k1,mnl1
        if(er(k2,1).lt.er(k2+1,1)) then
          c=er(k2+1,1)
          er(k2+1,1)=er(k2,1)
          er(k2,1)=c
          do 15 i=1,mnl
            c=s(i,k2+1)
            s(i,k2+1)=s(i,k2)
            s(i,k2)=c
  15      continue
        endif
  20  continue
      er(1,2)=er(1,1)
      do 30 i=2,mnl
        er(i,2)=er(i-1,2)+er(i,1)
  30  continue
      do 40 i=1,mnl
        er(i,3)=er(i,1)*100./tr
        er(i,4)=er(i,2)*100./tr
  40  continue
      return
      end

c-----
      subroutine tcoeff(m,n,mnl,f,s,er,egvt,ecof)
!-----input arrays
      dimension f(m,n),s(mnl,mnl),er(mnl,4)
!-----work arrays
      dimension v(m),v1(n)
!-----output arrays
        dimension egvt(m,mnl),ecof(mnl,n)

      do j=1,mnl
        do i=1,m
          egvt(i,j)=0.
        enddo
        do i=1,n
          ecof(j,i)=0.
        enddo
      enddo
c-----Normalizing the input eignvectors (c=sqrt(lamda))
      do 10 j=1,mnl
        c=0.
        do i=1,mnl
          c=c+s(i,j)*s(i,j)
        enddo
        c=sqrt(c)
        do i=1,mnl
          s(i,j)=s(i,j)/c
        enddo
  10  continue
c-----(m<n)
      if(m.le.n) then      
          do js=1,mnl
          do i=1,m
            egvt(i,js)=s(i,js)
          enddo
        enddo
        do 30 j=1,n
          do i=1,m
            v(i)=f(i,j)
          enddo
          do is=1,mnl
            do i=1,m
              ecof(is,j)=ecof(is,j)+v(i)*egvt(i,is)
            enddo
          enddo
  30    continue
      else
c-----(m>n)
        do 40 i=1,m
          do j=1,n
            v1(j)=f(i,j)
          enddo
          do js=1,mnl
          do j=1,n
            egvt(i,js)=egvt(i,js)+v1(j)*s(j,js)
          enddo
          enddo
  40    continue

        do 50 js=1,mnl
          do j=1,n
            ecof(js,j)=s(j,js)*sqrt(abs(er(js,1)))
          enddo
          do i=1,m
            egvt(i,js)=egvt(i,js)/sqrt(abs(er(js,1)))
          enddo
  50    continue       

      do js=1,mnl
          do i=1,m
          egvt(i,js)=(sqrt(abs(er(js,1)/(mnl*1.0))))*egvt(i,js)
          end do
          do i=1,n
            ecof(js,i)=(ecof(js,i))/(sqrt(abs(er(js,1)/(mnl*1.0))))
          end do
        end do
        end if
      return
      end

c-----
        subroutine rotaor(m,n,mnl,np,egvt,ecof,tr,rer,regvt,recof)
c-----input arrays
        dimension egvt(m,mnl),ecof(mnl,n)         
c-----work arrays       
        dimension h(m),st(np),vrlv(50),er(mnl,4),s(m,np)
c-----output arrays         
      dimension regvt(m,np),recof(np,n),rer(np,2)
      integer t,k

      write(*,*)'rotate beginning'
c-----Take out forward NP eigenvectors to rotate
      do j=1,np
        do i=1,m
            regvt(i,j)=egvt(i,j)
          end do
        end do
c-----Take out forward NP time coefficients to rotate
      do i=1,np
        do j=1,n
                recof(i,j)=ecof(i,j)
        end do
        end do
c-----compute the common degree
      do i=1,m
          h(i)=0.0
        do j=1,np
                h(i)=h(i)+regvt(i,j)**2
        end do
        end do
      do i=1,m
          h(i)=sqrt(h(i))
      end do
c-----compute the variance of the variance contribution by per common degree
c-----(equal to the ratated one)
      vrlv0=0.0
      do k=1,np
          g1=0.0
        g2=0.0
        do i=1,m
                g1=g1+(regvt(i,k)**2/h(i)**2)**2/real(m)
                g2=g2+(regvt(i,k)**2/h(i)**2)/real(m)
        end do
          g2=g2**2
          vrlv0=vrlv0+g1-g2
      end do
c-----rotated
      lll=0
        write(*,*)'rotate times'
800  continue
      do 505 t=1,np
      do 505 k=1,np
        if(t.ge.k) goto 505
                call rot(regvt,recof,h,t,k,m,n,np)
505  continue
      lll=lll+1
        write(*,*)'LLL=',LLL
      vrlv(lll)=0.0
      do k=1,np
          g1=0.0
          g2=0.0
        do i=1,m
            g1=g1+(regvt(i,k)**2/h(i)**2)**2/real(m)
            g2=g2+(regvt(i,k)**2/h(i)**2)/real(m)
        end do
          g2=g2**2
          vrlv(lll)=vrlv(lll)+g1-g2
      end do
      if(lll.lt.2) goto 800
        ci=(vrlv(lll)-vrlv(lll-1))/vrlv0
        if(ci.lt.0.001) goto 600
        if(lll.ge.40) goto 600
      goto 800
600  continue
c-----compute the explained variances
        do j=1,np
          st(j)=0.0
        do i=1,m
                st(j)=st(j)+regvt(i,j)**2
        end do
      end do
      do j=1,np
          rer(j,1)=st(j)*100/tr
        end do
c-----compute the accumulated explained variances
      ddd=0.0
      do k=1,np
        ddd=ddd+st(k)/tr
        rer(k,2)=ddd*100.
      enddo
      return
      end

c-----
      subroutine rot(a,b,h,t,k,m,n,np)
c-----input arrays
      dimension a(m,np),b(np,n),h(m)
c-----work arrays
      dimension u(m),vr(m),wt(m),wk(m),wbt(n),wbk(n)
c-----output arrays
c      dimension a(m,np),b(np,n)
      integer t,k

c-----compute fi
      do 25 i=1,m
      u(i)=(a(i,t)/h(i))**2-(a(i,k)/h(i))**2
25   vr(i)=2.0*(a(i,t)/h(i))*(a(i,k)/h(i))
      c=0.0
      d=0.0
      s=0.0
      g=0.0
      do 30 i=1,m
      c=c+(u(i)**2-vr(i)**2)
      d=d+2.0*u(i)*vr(i)
      s=s+u(i)
30   g=g+vr(i)
      tg1=d-2.0*s*g/(m*1.)
      tg2=c-(s**2-g**2)/(m*1.)
      fi=atan2(tg1,tg2)/4.0
c-----compute new a and b with fi
      do 75 i=1,m
      wt(i)=a(i,t)
75   wk(i)=a(i,k)
      do 84 kk=1,n
      wbt(kk)=b(t,kk)
84   wbk(kk)=b(k,kk)
      do 78 i=1,m
      a(i,t)=wt(i)*cos(fi)+wk(i)*sin(fi)
78   a(i,k)=wt(i)*(-sin(fi))+wk(i)*cos(fi)
      do 89 it=1,n
      b(t,it)=wbt(it)*cos(fi)+wbk(it)*sin(fi)
89   b(k,it)=wbt(it)*(-sin(fi))+wbk(it)*cos(fi)
      return
      end

c-----
      subroutine arrange2(m,n,np,regvt,recof,rer)
c-----input and output arrays         
      dimension regvt(m,np),recof(np,n),rer(np,2)

      mnl1=np-1
      do 20 k1=mnl1,1,-1
      do 20 k2=k1,mnl1
        if(rer(k2,1).lt.rer(k2+1,1)) then
          c=rer(k2+1,1)
          rer(k2+1,1)=rer(k2,1)
          rer(k2,1)=c
          do i=1,m
            d=regvt(i,k2+1)
            regvt(i,k2+1)=regvt(i,k2)
            regvt(i,k2)=d
          end do
          do i=1,n
            e=recof(k2+1,i)
            recof(k2+1,i)=recof(k2,i)
            recof(k2,i)=e
          end do
        endif
  20  continue
      do i=1,np
          rer(i,2)=0.0
        end do
      rer(1,2)=rer(1,1)
      do 30 i=2,np
        rer(i,2)=rer(i-1,2)+rer(i,1)
  30  continue
      return
      end


      subroutine meanvar(n,x,ax,sx,vx)
!-----input array
      dimension x(n)

      ax=0.
      vx=0.
      sx=0.
      do 10 i=1,n
        ax=ax+x(i)
  10  continue
      ax=ax/float(n)
      do 20 i=1,n
        vx=vx+(x(i)-ax)**2
  20  continue
      vx=vx/float(n)
      sx=sqrt(vx)
      return
      end

陶蜜云

我的天哪，好长啊