按照对OLR进行EOF分析后PC1正负位相筛选台风问题

什夜0.0

请教各位大佬，想按照OLR进行EOF分析后PC1的正负位相筛选台风，实心是正位相台风位置，空心是负位相台风位置，但在１５０E以东几乎没有台风生成，这是为什么呢？导师说应该在０－１８０E区间都有台风生成，原始数据没有问题可能是我进行分类时的程序有问题，请教应该如何改进程序呢？如有回答感激不尽！我的程序如下：
begin
  neof   = 1
;*******************************************
;  Read data
;*******************************************

  f    = addfile ("/home/mars/Liubs/30-60olr.nc", "r")
  time = (f->time)

   time@units  = "hours since 1800-01-01 00:00:0.0"
   time_ut  = cd_calendar(time,0)

   month   = time_ut(:,1)  ; Convert to integer for

   indtime0 = ind((month.ge.6).and.(month.le.11))

  x    = f->olr(indtime0,:,:)           ;
  ;printVarSummary(x)

;*******************************************
;  EOF
;*******************************************
  X = x(lat|:,lon|:,time|:)              ; Space x Time

  optEof = True
  eof    = eofunc_Wrap( X, neof, optEof)
  eof_ts = eofunc_ts_Wrap( X, eof, False)

  printVarSummary(eof)
  printVarSummary(eof_ts)


f1  = addfile ("/home/mars/Liubs/MJO_PC_INDEX.nc", "r")

time1 = (f1->time)

   time1@units  = "hours since 1800-01-01 00:00:0.0"
   time_ut1  = cd_calendar(time1,0)
   year1   = time_ut1(:,0)  ; Convert to integer for
   month1  = time_ut1(:,1)  ; Convert to integer for
   day1    = time_ut1(:,2)

  PC1 = (f1->PC1)
  f2 = addfile("/home/mars/Liubs/lunwen/data/IBTrACS.WP.v04r00.nc","r")

   time2 = (f2->time)
   time2@units  = "days since 1858-11-17 00:00:00"
   time_ut2  = cd_calendar(time2,0)  
   year2   = time_ut2(:,0,0)  ; Convert to integer for
   ;print(year2)
   month2  = time_ut2(:,0,1)  ; Convert to integer for
   day2    = time_ut2(:,0,2)

  indyear2 = ind((year2.ge.1979).and.(year2.le.2020).and.(month2.ge.6).and.(month2.le.11))
  indmonth2=ind((year2.ge.1979).and.(year2.le.2020).and.(month2.ge.6).and.(month2.le.11))
  indday2  = ind((year2.ge.1979).and.(year2.le.2020).and.(month2.ge.6).and.(month2.le.11))
  printVarSummary(indyear2)
  year2  := year2(indyear2)
  month2 := month2(indmonth2)
  day2   := day2(indday2)


    indTC1 = new((/1220/),float)
do nTC = 0,1219
    anyTC1 = any((year1.eq.year2(nTC)).and.(month1.eq.month2(nTC)).and.(day1.eq.day2(nTC)).and.(PC1.gt.0))

    if (anyTC1) then
      indTC1(nTC) = nTC
   end if
  delete([/anyTC1/])
end do
indTC1  := indTC1(ind(.not.ismissing(indTC1)))

   indTC1!0    = "time"

   indTC1&time = time2(indyear2(ind(.not.ismissing(indTC1))),0)
   indtime1 = toint(indTC1)
   ;printVarSummary(indTC1)
   lat1 = (f2->lat(indtime1,0))
   lon1 = (f2->lon(indtime1,0))
   ;print(lon1)

   indTC2 = new((/1220/),float)
do nTC = 0,1219
    anyTC2 = any((year1.eq.year2(nTC)).and.(month1.eq.month2(nTC)).and.(day1.eq.day2(nTC)).and.(PC1.lt.0))

    if (anyTC2) then
      indTC2(nTC) = nTC
   end if
  delete([/anyTC2/])
end do
indTC2  := indTC2(ind(.not.ismissing(indTC2)))

   indTC2!0    = "time"

   indTC2&time = time2(indyear2(ind(.not.ismissing(indTC2))),0)
   indtime2 = toint(indTC2)
   ;printVarSummary(indTC2)
   lat2 = (f2->lat(indtime2,0))
   lon2 = (f2->lon(indtime2,0))
   ;print(lon2)
;*******************************************
;  plots
;*******************************************
  ;---TC
  wks  = gsn_open_wks("png","MJO-pc1+-&TC3 ")  

        res                = True
        res@gsnDraw        = False
        res@gsnFrame       = False
        res@gsnMaximize    = True
        res@tmXTOn         = False
        res@tmYROn         = False
        res@gsnLeftString  = "MJO"
        res@gsnRightString = ""

;;set map;;

        res@gsnDraw        = False
        res@gsnFrame       = False
        res@gsnAddCyclic                = False
        res@mpOutlineOn                 = True
        res@mpGeophysicalLineThicknessF = 2
        res@mpNationalLineThicknessF    = 2
        res@mpFillDrawOrder             = "PostDraw"
        res@mpFillOn                    = False


;;set area;;

        res@mpMinLatF                   = 0
        res@mpMaxLatF                   = 50
        res@mpMinLonF                   = 100
        res@mpMaxLonF                   = 180

;;set contour;;

        res@cnFillDrawOrder             = "PreDraw"
        res@cnFillOn                    = True
        res@cnLinesOn                   = False
        res@pmLabelBarWidthF            = 0.6
        res@pmLabelBarHeightF           = 0.08
        res@pmLabelBarOrthogonalPosF    = 0.1
        res@lbLabelFontHeightF          = 0.01
        res@lbLabelAngleF               = 0
        res@cnLevelSelectionMode   =  "ManualLevels"
        res@cnMinLevelValF         =  -0.1
        res@cnMaxLevelValF         =  0.12
        res@cnLevelSpacingF        =  0.02
        res@tiMainString         =  ""
; Older way to subset a color map
       res@cnFillPalette               = "BlueDarkRed18"
       res@gsnLeftString               = ""

      txres                       = True
      txres@txFontHeightF         = 0.02
      txres@txFontColor           = "black"
      txres@txPerimOn               = True
      txres@txPerimColor            = "black"
      txres@txBackgroundFillColor = "White"
      txres@txFontOpacityF        = 0.8
      txres@txFontThicknessF      = 5.0

   contour = gsn_csm_contour_map(wks,eof(0,:,:),res)
  rts                        = True
  rts@gsMarkerIndex = 16
  rts@gsMarkerColor = "purple4"
  rts@gsMarkerSizeF = 0.01            
  rts@gsMarkerThicknessF = 3
  rts@tfPolyDrawOrder =  "Draw"
  rps                        = True
  rps@gsMarkerIndex = 4
  rps@gsMarkerColor = "purple4"
  rps@gsMarkerSizeF = 0.01            
  rps@gsMarkerThicknessF = 3
  rps@tfPolyDrawOrder =  "Draw"


   dum = gsn_add_text(wks,(/contour/),"(a)PC1&TC",107.5,47,txres)


   points1 = gsn_add_polymarker(wks,contour,lon1,lat1,rts)
   points2 = gsn_add_polymarker(wks,contour,lon2,lat2,rps)
    draw(contour)
    frame(wks)

end

JGD飘

似乎你的indtime1，并不是正负位相的storm，还需要再ind一次，拙见！

什夜0.0

JGD飘 发表于 2022-6-10 16:13
似乎你的indtime1，并不是正负位相的storm，还需要再ind一次，拙见！

确实确实是这里筛选有问题，已经解决啦 感谢回复！