使用NCL沿着任意两点做垂直剖面图

BLUEANGEL2012

请问楼主解决了吗

xuebiz

仅供参考~~~~~~

也是抄的官网例子transect_1.ncl改的，注意看后面注释部分。。。

前面读nc文件，读温度数据请自动忽略。。。反正你也没俺用的数据。。。
数据是这样子的：    short sz_h(time, dep, lat, lon) ;

;************************************
; CSM_Graphics: transect_1.ncl
;************************************
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"
;************************************
begin
;  in = addfile("h_avg_Y0191_D000.00.nc","r")
  in = addfile("../data/G2016042812.nc", "r")
  t  = in->sz_h(0,:,:,:)

;************************************
; calculate great circle along transect
;************************************
  leftlat  =  10.0
  rightlat =  30.0

  leftlon  =  90.0
  rightlon =  150.0

  npts     =   10                    ; number of points in resulting transect

  dist     = gc_latlon(leftlat,leftlon,rightlat,rightlon,npts,2)
  points   = ispan(0,npts-1,1)*1.0
;printVarSummary(t)
;********************************
; interpolate data to great circle
;********************************
;  trans   = linint2_points(t&lon_t,t&lat_t,t,True,dist@gclon,dist@gclat,2)
  trans   = linint2_points(t&lon,t&lat,t,True,dist@gclon,dist@gclat,2)
  copy_VarAtts(t,trans)          ; copy attributes

;  trans!0      = "z_t"           ; create named dimension and assign
;  trans&z_t    = t&z_t           ; coordinate variable for 0th dimension only
  trans!0      = "dep"
  trans&dep    = t&dep
  trans@_FillValue = t@missing_value
  trans@scale_factor = t@scale_factor
  trans = trans*trans@scale_factor
;printVarSummary(trans)
;********************************
; create plot
;********************************
  wks = gsn_open_wks("png","trans_salt")       ; open a ps file
;  gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; choose colormap
  gsn_define_colormap(wks,"rainbow")
  res                     = True          ; plot mods desired
  res@tmXBMode            = "Explicit"    ; explicitly label x-axis
  res@tmXBValues          = (/points(0),points(npts-1)/) ; points to label
; label values
  res@tmXBLabels          = (/leftlat +", "+leftlon,rightlat+", "+rightlon/)

  res@cnFillOn            = True         ; turn on color
  res@lbLabelAutoStride   = True         ; nice label bar label stride
  res@gsnSpreadColors     = True         ; use full range of colormap
  res@cnLinesOn           = True        ; turn off countour lines
  res@lbOrientation       = "vertical"   ; vertical label bar
  res@pmLabelBarOrthogonalPosF = -0.05        ; move label bar closer to plot

  res@tiMainString        = "Transect"   ; add title
  res@tiXAxisString       = "lat/lon along transect"
  res@trYReverse          = True         ; reverse y axis
;  res@trXReverse          = True         ; reverse x axis (neg longitudes)
;  res@cnLevelSpacingF     = 1.0          ; set contour spacing
  
  plot = gsn_csm_contour(wks,trans,res)  ; create plot

end

吴阿好大人

xuebiz 发表于 2016-9-21 08:27
仅供参考~~~~~~

也是抄的官网例子transect_1.ncl改的，注意看后面注释部分。。。

谢谢！让我研究研究

xuebiz

吴阿好大人 发表于 2016-9-21 08:29
谢谢！让我研究研究

http://www.ncl.ucar.edu/Applications/transect.shtml

这是我用的NCL官网的例子，看看适用你的情况不？

沐..

不知道楼主解决没 这个是我的部分脚本  沿着特定两点剖面图   
http://www.ncl.ucar.edu/Applications/transect.shtml    这个是官网参考例子
  f0=addfile(dir+"4.grib2","r")
  lat     =f0->lat_0
  lat=-lat              ; coordinates
  lon     =f0->lon_0
  p       =f0->lv_ISBL0

  mm=dimsizes(lat)
  nn=dimsizes(lon)
  

  fenzu=dir2+"fenzu.txt"
  nrowfenzu=numAsciiRow(fenzu)
  datafenzu=asciiread(fenzu, -1, "string")

  ngr=7
  gr= new(ngr, "integer")
  gr=stringtointeger(str_get_field(datafenzu,2, " "))
   
  itime=1
  wks  = gsn_open_wks("ps","sub_ty_poor")  
  i=0

  f=addfile(dir+gr(i)+".grib2","r")     ;time,level,lat,lon
  u=f->u_P1_L100_GLL0(itime,:,:,:)
  v=f->v_P1_L100_GLL0(itime,:,:,:)

  scale = 1.e05                                 
  vor   = u                                      ; retain coordinates
  vor   = uv2vrG_Wrap(u,v) * scale
  
  mn=(/8,mm,nn/)
  t=new(mn,"float")
  
  t(:,:,:)=vor(:,360:0,:)
  lat@units="degrees_north"
  lon@units="degrees_east"
  lon@units="level"

  t!2="lon"
  t!1="lat"
  t!0="p"


  t&lat=lat
  t&lon=lon
  t&p=p

  t@long_name = "vorticity"
  t@units     = "1.e05"
;************************************************
;;;;;;;;开始沿着副中心到台风中心差值;;;;;;;;;;;;;;;;;;;;;;;;;;
;************************************
  leftlat  =  25.3
  rightlat =  22.5

  leftlon  =  118
  rightlon =  124
  npts     =  13            ; number of points in resulting transect
  dist     = gc_latlon(leftlat,leftlon,rightlat,rightlon,npts,2)
  points   = ispan(0,npts-1,1)*1.0
;********************************
; interpolate data to great circle
;********************************
  trans   = linint2_points(t&lon,t&lat,t,True,dist@gclon,dist@gclat,2)
  copy_VarAtts(t,trans)          ; copy attributes

  trans!0      = "p"           ; create named dimension and assign
  trans&p   = t&p   


  res                   = True  
  res@gsnDraw = False
  res@gsnFrame = False                 ; plot mods desired
  res@cnFillOn          = True  

  res                     = True          ; plot mods desired
  res@tmXBMode            = "Explicit"    ; explicitly label x-axis
  res@tmXBValues          = (/points(0),points(npts-1)/) ; points to label
; label values
  res@tmXBLabels          = (/leftlat +", "+leftlon,rightlat+", "+rightlon/)

                 ; turn on color               ; use full range of color map
  res@lbLabelAutoStride = True                   ; nice label bar labels
  res@cnLinesOn         = False                  ; no contour lines
  res@mpMinLonF =118
  res@mpMaxLonF = 124
  res@gsnAddCyclic = False
  res@mpGridAndLimbOn =True

  ;res@tmYLMode   ="Explicit"
  ;res@tmYLTickSpacingF=2
  res@gsnLeftString="forecast_time="+itime*6+gr(i)
  res@gsnRightString="lon ="+lontc(itime)
; res@cnFillPalette              
  res@cnSmoothingOn               =True
  res@cnSmoothingDistanceF        =0.0005
  res@cnLevelSelectionMode = "ManualLevels"       ; set manual contour levels
  res@cnMinLevelValF       = -30                 ; set min contour level
  res@cnMaxLevelValF       =  30                  ; set max contour level
  res@cnLevelSpacingF      =  5                 ; set contour spacing
  res@cnFillPalette               = "GMT_polar"


  plot =gsn_csm_pres_hgt(wks,trans,res)  

   draw(wks)
   frame(wks)
  delete([/f,u,v,ftc,lattc,lontc/])
end

yeah...

沐.. 发表于 2016-10-12 16:27
不知道楼主解决没 这个是我的部分脚本  沿着特定两点剖面图   
http://www.ncl.ucar.edu/Applications/tra ...

请教一下那个npts=13是如何计算得到的呢   还有我看你说开始沿着副中心到台风中心插值  我不太懂 求教了

yeah...

沐.. 发表于 2016-10-12 16:27
不知道楼主解决没 这个是我的部分脚本  沿着特定两点剖面图   
http://www.ncl.ucar.edu/Applications/tra ...

请教一下那个npts=13是如何计算得到的呢   还有我看你说开始沿着副中心到台风中心插值  我不太懂 求教了

沐..

yeah... 发表于 2016-10-12 21:07
请教一下那个npts=13是如何计算得到的呢   还有我看你说开始沿着副中心到台风中心插值  我不太懂 求教了

13是自己订的，根据你的网格分辨率和计算范围来定，一般不要多于本身的差值格点吧。我是沿着两个台风中心求画的坡面图  

yeah...

沐.. 发表于 2016-10-16 14:08
13是自己订的，根据你的网格分辨率和计算范围来定，一般不要多于本身的差值格点吧。我是沿着两个台风中心 ...

谢谢您的回复。我利用npts设置了不同的数值画了一下，发现那个npts越大极值中心越明显呀  还想问您一下  在画剖面图的时候  您考虑了地形吗？

沐..

这个没有哦    不过差值的补数应该最好不要跟原来的分辨率相差太多