关于ncl添加地形

karenlk

1649518749 发表于 2016-8-18 14:16
不行啊

不会啊，什么问题

紫儿

楼主，最近正在画剖面图，我试过一些方法添加不了，所以想看看你的脚本，对比一下，但贡献被我用完了，麻烦能把etopo-xz.ncl这个脚本发一下我吗，853619839@qq.com,谢谢！

刘振

不用res@trYReverse        = True就可以画成楼主的效果的吗，我的默认1000hPa都是在上面的，overlay之后就没地形了

沐..

karenlk 发表于 2016-6-7 19:17
改进剖面地形添加

楼主能否把你这个给我发一下啊   唯一的一个贡献下载了你前面5-4的那个脚本 但是运行错误哎

沐..

楼主那沿着任意两点的垂直剖面图如何画呢，以下是我的脚本   单独可以画出地形高度，但是没有没办法overlay~你知道是什么原因吗

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"


undef("read_elev_data")
function read_elev_data(topo_file)
local nlat, nlon, topo_file, lat, lon
begin
;---Read data as a straight binary file
  nlat = 2160
  nlon = 4320
  setfileoption("bin","ReadByteOrder","BigEndian")
  elev = cbinread(topo_file,(/nlat,nlon/),"short")

;---Create 1D coordinate arrays
  lat       = fspan(90,-90,nlat)
  lon       = fspan(0,360,nlon)
  lat!0     = "lat"
  lon!0     = "lon"
  lat@units = "degrees_north"
  lon@units = "degrees_east"
  lat&lat   = lat
  lon&lon   = lon

;---Attach the coordinate arrays
  elev!0    = "lat"
  elev!1    = "lon"
  elev&lat  = lat
  elev&lon  = lon

  return(elev)
end

;----------------------------------------------------------------------
; Main code
;

begin
lon       = fspan(0,360,4320)
dir  ="/mnt/hgfs/res_gong/ec/2009/ec/0712/data/"
dir2  ="/mnt/hgfs/res_gong/ec/2009/ec/0712/"
dir3  ="/mnt/hgfs/res_gong/ec/2009/ec/0712/surface/"
dir4="/mnt/hgfs/res_gong/"

  f0=addfile(dir+"4.grib2","r")
  lattc     =f0->lat_0             ; coordinates
  lontc     =f0->lon_0
  p       =f0->lv_ISBL0

  mm=dimsizes(lattc)
  nn=dimsizes(lontc)
  

  setfileoption("nc","MissingToFillValue",False)


elev= read_elev_data(dir4+"ETOPO5.dat")
elev0 = 1013.25*(1-elev*0.0065/288.15)^5.25145        ;½«µØÐÎÊý¾ÝµÄµ¥Î»£ºmת»¯ÎªhPa
copy_VarCoords(elev,elev0)

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

  ngr=6
  gr= new(ngr, "integer")
  gr=stringtointeger(str_get_field(datafenzu,1, " "))

  ;i=1
  itime=0
  do i=0,6
  wks  = gsn_open_wks("x11","g+_sub_ty_volo_"+itime*6+"time"+gr(i))  
  ;;;;;;;;读入差值点信息;;;;;;;;;;;;;;
  ftc=dir3+gr(i)+".txt"
  datatc=asciiread(ftc,-1, "string")
  lat1=stringtofloat(str_get_field(datatc,1, ","))
  lon1=stringtofloat(str_get_field(datatc,2, ","))
  lat2=stringtofloat(str_get_field(datatc,3, ","))
  lon2=stringtofloat(str_get_field(datatc,4, ","))


  f=addfile(dir+gr(i)+".grib2","r")     ;time,level,lat,lon
  u=f->u_P1_L100_GLL0(itime,:,mm-1:0,:)
  v=f->v_P1_L100_GLL0(itime,:,mm-1:0,:)
  scale = 1.e05                                 
  t  = u                                      ; retain coordinates
  t=uv2vrF_Wrap(u,v)*scale
;************************************************
;;;;;;;;开始计算差值经纬度;;;;;;;;;;;;;;;;;;;;;;;;;;
;************************************


  k=(lat1(itime)-lat2(itime))/(lon1(itime)-lon2(itime))
  b=lat1(itime)-k*lon1(itime)
  leftlon  =  118
  leftlat  =k*leftlon+b
  rightlon =  124
  rightlat =k*rightlon+b
  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
  trans   = linint2_points(t&lon_0,t&lat_0,t,True,dist@gclon,dist@gclat,2)
  copy_VarAtts(t,trans)          ; copy attributes
  trans!0      = "p"           ; create named dimension and assign
  trans&p   = t&lv_ISBL0   
;;;;;;开始画图;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

elon=dist@gclon
elat=dist@gclat
print(elat)

elev2=new(13,"float")
do ie=0,12
elev2(ie)=elev0({elat(ie)},{elon(ie)})        ;¹Ì¶¨Ä³Ò»Î³¶ÈµÄµØÐÎÊý¾Ý
print(elev2(ie))
end do
elev2@_FillValue= -9.96921e+36
print(elon)



  res                   = True  
  res@gsnDraw           = False
  res@gsnFrame          = False                 ; plot mods desired
  res@cnFillOn          = 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@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@tmXBMode   ="Explicit"
  res@tmXBTickSpacingF=2

  res@tiMainString="forecast_time="+itime*6
  res@gsnRightString="TC lon="+lon1(itime)+" SC lont="+lon2(itime)
; res@cnFillPalette              
  res@cnSmoothingOn               =True
  res@cnSmoothingDistanceF        =0.0005
  res@cnLevelSelectionMode = "ManualLevels"       ; set manual contour levels
  res@cnMinLevelValF       =0                ; set min contour level
  res@cnMaxLevelValF       =80                ; set max contour level
  res@cnLevelSpacingF      =5                ; set contour spacing
  res@cnFillPalette               = "precip3_16lev";"GMT_polar"
    res@cnFillDrawOrder      = "Predraw"
  plot =gsn_csm_pres_hgt(wks,trans,res)  


res2=True
res2@gsnMaximize      = True
res2@trYMaxF=1000
res2@gsnYRefLine      = 1000.          ; create a reference line
res2@gsnBelowYRefLineColor= "black"              ; above ref line fill red
res2@gsnDraw        = False             ; don't draw the plots yet
res2@gsnFrame       = False
    res2@cnFillDrawOrder      = "Predraw"
plot_hgt = gsn_csm_xy(wks,elon,elev2,res2)  
overlay(plot,plot_hgt)

;************************************************
  draw(wks)
  frame(wks)
  delete([/f,u,v/])
  delete([/ftc,datatc,lat1,lat2,lon1,lon2/])
end do
end

karenlk

沐.. 发表于 2016-10-25 15:55
楼主能否把你这个给我发一下啊   唯一的一个贡献下载了你前面5-4的那个脚本 但是运行错误哎

邮箱说一下

沐..

karenlk 发表于 2016-10-25 18:39
邮箱说一下

342082431@qq.com   谢谢~

蚂蚁

赞一个，积蓄分享!!!!!

Eric.Yang

{:eb502:}

蚂蚁

不知道为什么无法下载（自动回复：请不要使用迅雷等下载工具，点我查看下载帮助）楼主挂出来的文件