爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 1578|回复: 0

wrfout用NCL画台风路径

[复制链接]
发表于 2020-12-10 10:30:43 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
程序如图怎么改?求助大佬
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
; ###To get lat/lon info.

  a = addfile("wrfout_d01_2017-08-19_00ls:00:00.nc","r")
  lat2d = a->XLAT(0,:,:)
  lon2d = a->XLONG(0,:,:)
  dimll = dimsizes(lat2d)
  nlat  = dimll(0)
  mlon  = dimll(1)
  t = wrf_user_getvar(a,"times",-1)
  dimt = dimsizes(t)
;print(t)
;print(dimt)
ndate = 45
; Sea Level Pressure
  slp = wrf_user_getvar(a,"slp",0)
  dims = dimsizes(slp)
; Array for track
  time = new(ndate,string)
  imin = new(ndate,integer)
  jmin = new(ndate,integer)
  smin = new(ndate,integer)
  slpmin = new(ndate,float)
  newlon = new(ndate,float)
  newlat = new(ndate,float)


;###########################
do fi=4,48
slp2d = wrf_user_getvar(a,"slp",fi)
maxlat = 23.5
minlat = 9.04
maxlon = 130.96
minlon = 109.00
opt = True
   loc  = wrf_user_ll_to_xy(a,(/minlon,maxlon/),(/minlat,maxlat/),opt)
   latlon = wrf_user_xy_to_ll(a,loc(0,:),loc(1,:),opt)
   loc = loc-1
  slpnew = slp2d(loc(1,0):loc(1,1),loc(0,0):loc(0,1))
  dimss = dimsizes(slpnew)

; We need to convert 2-D array to 1-D array to find the minima.
    slp1dnew     = ndtooned(slpnew)
        printVarSummary(slp1dnew)
    slp1d = ndtooned(slp2d)
        printVarSummary(slp1d)
        print(fi)
    smin(fi) = ind(slp1d.eq.min(slp1dnew))
    slpmin(fi) = min(slp1dnew)
; Convert the index for 1-D array back to the indeces for 2-D array.
    minij     = ind_resolve(ind(slp1d.eq.min(slpnew)),dims)
    imin(fi) = minij(0,0)
    jmin(fi) = minij(0,1)
    tlon = lon2d(minij(0,0),minij(0,1))
end do
;###########################

; Graphics section

  wks=gsn_open_wks("png","track")              ; Open PS file.
  gsn_define_colormap(wks,"BlGrYeOrReVi200")  ; Change color map.

  res                     = True
  res@gsnDraw             = False             ; Turn off draw.
  res@gsnFrame            = False             ; Turn off frame advance.
  res@gsnMaximize         = True              ; Maximize plot in frame.

  res@tiMainString = "HATO ";"Hurricane Isabel"       ; Main title

  WRF_map_c(a,res,0)                          ; Set up map resources
                                              ;    (plot options)
  plot = gsn_csm_map(wks,res)                 ; Create a map.

; Set up resources for polymarkers.
  gsres                = True
  gsres@gsMarkerIndex  = 16                  ; filled dot
  ;gsres@gsMarkerSizeF = 0.005               ; default - 0.007
  cols                  = (/5,160,40/)

; Set up resources for polylines.
  res_lines                      = True
  res_lines@gsLineThicknessF     = 3.           ; 3x as thick

  dot  = new(ndate,graphic)    ; Make sure each gsn_add_polyxxx call
  line = new(ndate,graphic)    ; is assigned to a unique variable.

; Loop through each date and add polylines to the plot.
  do i = 0,ndate-2
     res_lines@gsLineColor           = cols(0)
     xx=(/lon2d(imin(i),jmin(i)),lon2d(imin(i+1),jmin(i+1))/)
     yy=(/lat2d(imin(i),jmin(i)),lat2d(imin(i+1),jmin(i+1))/)
     line(i) = gsn_add_polyline(wks,plot,xx,yy,res_lines)
  end do
  do j = 0 ,ndate - 1
    newlon(j) = lon2d(imin(j),jmin(j))
    newlat(j) = lat2d(imin(j),jmin(j))
  end do
  asciiwrite ("track.txt" ,newlon +" "+ newlat)
  asciiwrite ("slp.txt" ,slpmin)
  lon1d = ndtooned(lon2d)
  lat1d = ndtooned(lat2d)

; Loop through each date and add polymarkers to the plot.
  do i = 0,ndate-1
     print("dot:"+lon1d(smin(i))+","+lat1d(smin(i)))
     gsres@gsMarkerColor  = cols(0)
     dot(i)=gsn_add_polymarker(wks,plot,lon1d(smin(i)),lat1d(smin(i)),gsres)
  end do

; Date (Legend)
  txres               = True
  txres@txFontHeightF = 0.015
  txres@txFontColor   = cols(0)
  txid1 = new(ndate,graphic)

  draw(plot)
  frame(wks)
end

路径这样怎么办

路径这样怎么办
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表