idl 的EOF分析

liang0816

最近在学些EOF分析，并且用idl做了一个程序，虽然出结果了，但是有问题，不知什么原因。
请各位给指导一下，谢谢！
pro eof_sample

    path =file_dirname(routine_filepath('eof_sample'))
  print,path
  cd,path
     file = file_search('spei.tif',count =n)
        print,file
        data =read_tiff(file,geotiff =geographicgeotypekey)      
        size =size(data)  
      print,size      
       ts =size[1]
       nc =size[2]
       nl =size[3]
       type =size[4]
    print,ts,nc,nl ,type
      data =reform(data,ts,nc*nl,/overwrite)
      data =reform(data,nc*nl,ts,/overwrite)     
      data_anomalies =float(data)
     for j =0,nc*nl-1 do begin
       data_anomalies[j,*]= (data[j,*]-mean(float(data[j,*])))/10000       ;距平处理
     endfor   
;
          matrix =(1.0/(ts-1))*(data_anomalies##transpose(data_anomalies)) ;协方差矩阵
     la_svd,matrix,W,U,V
     ; Calculate the eigenvectors and eigenvalues of the covariance matrix.
     ; W: ntime          奇异值
     ; U: ntime x ntime  左矩阵
     ; V: ntime x ntime  右矩阵
  help,W ,U,V  

    ;The first value IN the singular value VECTOR corresponds to the first mode, and,
    ;as a percentage OF the TOTAL values, tells you how much the first mode VECTOR explains the data VARIANCE.
    ;We sometimes compute AND display this as a percent VARIANCE.
     percent_variance =W^2/total(W^2)*100            ;w矩阵的值和矩阵的特征值λ 之间有个换算关系，w=sqrt（nλ )
print,percent_variance                              ;贡献
;
   dims = Size(data_anomalies, /Dimensions)
   eof = FltArr(dims[1], dims[0])
FOR j=0,dims[1]-1 DO BEGIN
     t = Transpose(data_anomalies) ## U[j,*]                                           ;矩阵乘,得到第j个模态向量
     eof[j,*]  = t / SQRT(Total(t^2))
ENDFOR

;
mode = 1                                       ;第一主成分
theEOF = eof[mode-1,*]
theEOF = Reform(theEOF,nc, nl,/overwrite)
help,theEOF
write_tiff,'theEOF.tif',theEOF,/float;,geotiff =geographicgeotypekey,/float
write_tiff,'data_anomalies.tif',data_anomalies,geotiff =geographicgeotypekey,/float
      pc = FltArr(dims[1], dims[1])
   FOR j=0,dims[1]-1 DO BEGIN
      pc[j,*]  = data_anomalies ## eof[j,*]
   ENDFOR     
      plot,pc[0,*]         
  end