NCL绘制台风路径密度（track density）——比想象的简单

Mid_Student

最近要计算台风路径密度，网上没有找到现成的程序，只能写一个，比想象的要简单。
使用的数据是中国气象局台风中心登陆中国的台风路径，关键是提取经纬度坐标。还请各位大佬批评指正。

begin

  lat0  = ispan(0,50,2)            ; 2度 格点
  lon0  = ispan(90,180,2)          ; 空间范围必须要覆盖所有台风路径经纬度，否则后面会出错
  
  nlat  = dimsizes(lat0)
  nlon  = dimsizes(lon0)  
   
num_tc = new((/nlat,nlon/),"float",-999)  

  num_tc!0    = "lat"
  num_tc!1    = "lon"
  num_tc&lat  = lat0
  num_tc&lon  = lon0   
  num_tc&lat@units  = "degrees_north"
  num_tc&lon@units  = "degrees_east"  
  
files   = "/mnt/e/work/TC_hazard/Track_density/tc_track_2000-2019.txt"

  var    = asciiread (files,-1,"integer")
  
     nt = num(var.eq.66666)        ; 台风个数
     nm = ind(var.eq.66666)        ; 定位每个台风
         
   num_tc = 0
           
do ii = 0,nt-1
          ; print(ii)
    if (ii.lt.nt-1) then
     var0 := var(nm(ii)+8:nm(ii+1)-1)  ;定位每个台风初始位置
    else
     var0 := var(nm(ii)+8::)
    end if
            
     nrow := dimsizes(var0)/6        
      Ntc := var0(1::6)
      Lat := var0(2::6)/10.
      Lon := var0(3::6)/10.
      
    do nn = 0,nrow-1
      if (.not.ismissing(Lat(nn)).and..not.ismissing(Lon(nn))) then
       num_tc({Lat(nn)},{Lon(nn)}) = (/num_tc({Lat(nn)},{Lon(nn)})+1/)   ；关键判断语句
      end if
    end do
end do
           num_tc = num_tc/20.                      ； 年平均
           num_tc = where(num_tc.le.0, num_tc@_FillValue, num_tc)
wks = gsn_open_wks("pdf", "track_density_2000-2019")
gsn_define_colormap(wks,"prcp_1")     

  res                = True


  res@gsnDraw        = False
  res@gsnFrame       = False
  res@cnFillOn       = True
  res@cnLinesOn      = False
  res@gsnAddCyclic   = False
  res@cnLineLabelsOn = False  
  res@tmXTOn         = False
  res@tmYROn         = False   
  res@mpMinLatF      = 10
  res@mpMaxLatF      = 40
  res@mpMinLonF      = 100
  res@mpMaxLonF      = 130
  res@tmXBLabelFontHeightF  = 0.02
  res@tmYLLabelFontHeightF  = 0.02
  res@mpOutlineOn           = True
  res@mpOutlineBoundarySets = "National"
  res@mpDataBaseVersion     = "MediumRes"
  res@mpDataSetName         = "/database/Earth..4"
  res@mpOutlineSpecifiers   = (/"China:States"/)
  res@mpFillOn              = False
  res@pmTickMarkDisplayMode = "Always"  


  res@cnFillMode            = "RasterFill"
  res@cnLevelSelectionMode  = "ExplicitLevels"   ; set explicit contour levels
  res@cnLevels              = (/0,1,2,3,4,5,6,7,8,9,10/)
  res@cnFillColors          = (/2,3,6,7,8,9,10,11,12,13,14,15/)
  res@gsnRightStringOrthogonalPosF  = 0.005
  res@gsnLeftStringOrthogonalPosF   = 0.005
  res@gsnLeftString         = " TC Track Density "   
  res@gsnRightString        = " 2000-2019 "  
      
  plot = gsn_csm_contour_map_ce(wks,num_tc,res)
  
    draw(plot)
    frame(wks)


end

折磨死

你好同学，首先感谢分享！我想问一下那个关键判断语句是什么意思呢？我没有看懂，而且用了您的程序和数据之后，也没有画出来图，输出了一下关键判断语句的左右变量，出来的值还是令我比较困惑。期待您的回复

Mid_Student

折磨死 发表于 2021-12-5 22:52
你好同学，首先感谢分享！我想问一下那个关键判断语句是什么意思呢？我没有看懂，而且用了您的程序和数据之 ...

程序我又执行了一次，没发现什么问题。判断语句是指单位格点内台风点（经纬度）的个数累加，一共20年的台风路径，最后除以20，得到平均每年的。

折磨死

Mid_Student 发表于 2021-12-7 15:04
程序我又执行了一次，没发现什么问题。判断语句是指单位格点内台风点（经纬度）的个数累加，一共20年的台 ...

好哒，谢谢您，是因为我没有注意到您提醒的设置格点的范围需要覆盖所有的路径，所以没有画出来。使用您的程序之后出的图没有问题，再次感谢！

DeQiang

学到了，谢谢楼主的分享