爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
楼主: 小其其格

NCL绘制wrfout风场垂直剖面并叠加地形

  [复制链接]

新浪微博达人勋

 楼主| 发表于 2019-4-28 08:59:44 | 显示全部楼层
雨萌萌 发表于 2019-4-27 15:53
楼主您好!我用了你的程序试着画了一下,可是地形与物理量对不上,如下图是怎么回事呢

在插值的时候,如果给定了线段AB,则opts要设置为True

opts = True
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-4-29 13:35:01 | 显示全部楼层
小其其格 发表于 2019-4-28 08:59
在插值的时候,如果给定了线段AB,则opts要设置为True

opts = True

感谢楼主回复!我就是设为的True,换了个经纬度范围,还是这样,是窜位的C:\Users\kristina\Desktop\plt_CrossSection4.000005
plt_CrossSection4.000005.png
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2019-4-30 19:37:37 | 显示全部楼层
本帖最后由 小其其格 于 2019-4-30 19:39 编辑
雨萌萌 发表于 2019-4-27 15:53
楼主您好!我用了你的程序试着画了一下,可是地形与物理量对不上,如下图是怎么回事呢


用以下代码绘制试试看,我的出图没问题: plt_UW_mdBz_08.ab3.ncl.000003.png

  1. ;   Example script to produce plots for a WRF real-data run,
  2. ;   with the ARW coordinate dynamics option.
  3. ;   Plot data on a cross section
  4. ;   This script will plot data at a set angle through a specified point
  5. ;   This script adds lon/lat info along X-axis

  6. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
  7. load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

  8. begin
  9. ;
  10. ; The WRF ARW input file.  
  11. ; This needs to have a ".nc" appended, so just do it.
  12. a = addfile("../wrfout_d03_2016-05-05_08_00_00.nc","r")


  13. ; We generate plots, but what kind do we prefer?
  14. ; type = "x11"
  15. ; type = "pdf"
  16. ; type = "ps"
  17. ; type = "ncgm"
  18. type = "png"
  19. ;wks = gsn_open_wks(type,"plt_UW_mdBz_08.ncl")
  20. ;wks = gsn_open_wks(type,"plt_UW_mdBz_08.ab.2.ncl")
  21. wks = gsn_open_wks(type,"plt_UW_mdBz_08.ab3.ncl")
  22. ;gsn_define_colormap(wks,"WhViBlGrYeOrReWh")       ; Overwrite the standard color map
  23. ;gsn_define_colormap(wks, "Radar2")  ;加载自己定义的色标
  24. gsn_define_colormap(wks, "Radar.CINRAD")  ;加载自己定义的色标


  25. ; Set some basic resources
  26. res = True
  27. res@MainTitle = "REAL-TIME WRF"
  28. res@Footer = False


  29. pltres = True
  30. mpres = True
  31. mpres@mpGeophysicalLineColor = "Black"
  32. mpres@mpNationalLineColor    = "Black"
  33. mpres@mpUSStateLineColor     = "Black"
  34. mpres@mpGridLineColor        = "Black"
  35. mpres@mpLimbLineColor        = "Black"
  36. mpres@mpPerimLineColor       = "Black"


  37. mpres@mpDataSetName               = "Earth..4"
  38. mpres@mpDataBaseVersion           = "MediumRes"
  39. mpres@mpOutlineOn                 = True
  40. mpres@mpOutlineSpecifiers         = (/"China:states","Taiwan","Zhejiang"/)
  41. mpres@mpGeophysicalLineThicknessF = 2
  42. mpres@mpNationalLineThicknessF    = 2
  43. mpres@mpNationalLineColor         = "black"
  44. mpres@mpGeophysicalLineColor      = "black"
  45. mpres@mpProvincialLineColor       = "black"
  46. mpres@mpProvincialLineThicknessF  = 2

  47. mpres@mpFillOn                    = True
  48. mpres@mpFillAreaSpecifiers        = (/"water",       "land" /)
  49. mpres@mpSpecifiedFillColors       = (/"skyblue1","white"/)
  50. mpres@mpMaskAreaSpecifiers        = (/"China:states","Taiwan","Zhejiang"/)


  51. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  52. FirstTime = True
  53. FirstTimeMap = True
  54. times  = wrf_user_getvar(a,"times",-1) ; get times in the file
  55. ntimes = dimsizes(times)         ; number of times in the file

  56. mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file
  57. dNames = getfilevardims(a,"P")    ; NCL 6.4.0 and earlier
  58. nd = dimsizes(mdims)              ; 查看mdims为几围数组
  59. print(dNames + " : " + mdims)     ; 查看各个变量维数

  60. xlat = wrf_user_getvar(a, "XLAT",0)
  61. xlon = wrf_user_getvar(a, "XLONG",0)
  62. ter = wrf_user_getvar(a, "HGT",0)


  63. u10 = wrf_user_getvar(a,"U10",0)    ; u at 10 m, mass point
  64. v10 = wrf_user_getvar(a,"V10",0)    ; v at 10 m, mass point

  65. u10@units = "m/s"
  66. v10@units = "m/s"

  67. printVarSummary(xlat)



  68. ;---------------------------------------------------------------

  69. do it = 0,ntimes-1             ; TIME LOOP

  70.     print("Working on time: " + times(it) )
  71.     res@TimeLabel = times(it)   ; Set Valid time to use on plots


  72.     z    = wrf_user_getvar(a, "z",it)      ; grid point height
  73.     dbz = wrf_user_getvar(a,"dbz",it)    ;dBz
  74.     mdbz = wrf_user_getvar(a,"mdbz",it)    ;dBz
  75.     tc   = wrf_user_getvar(a, "tc", it)
  76.     p  = wrf_user_getvar(a, "pressure",it) ;
  77.     rh = wrf_user_getvar(a,"rh",it)        ; relative humidity;
  78.     printVarSummary(tc)

  79.     u = wrf_user_getvar(a,"ua",it)
  80.     v = wrf_user_getvar(a,"va",it)
  81.     w = wrf_user_getvar(a,"wa",it)

  82.     if ( FirstTime ) then                ; get height info for labels
  83.       zmin = 0.
  84.       zmax = 15.                          ; We are only interested in the first 6km
  85.       nz   = floattoint(zmax + 1)
  86.   end if


  87.   ;---------------------------------------------------------------

  88.   do ip = 1, 1,1        ; we are doing 3 plots;我们这里只绘制2张图
  89.     ; all with the pivot point (plane) in the center of the domain
  90.     ; at angles 0, 45 and 90
  91.     ;  
  92.     ;                   |
  93.     ;       angle=0 is  |          ;垂直剖面直线AB的形状
  94.     ;                   |
  95.     ;

  96.     ;线段AB的起点和终点经纬度设置

  97.    


  98.    
  99.     ;08.ab UTC
  100.     startlon = 119.0
  101.     startlat = 28.6

  102.     endlon   = 120.0
  103.     endlat   = 29.35

  104.     ;startlat = 29.35
  105.     ;endlat   = 28.6

  106.   

  107.     PI=3.1415926
  108.     ;--------------------------------------------------------------------------
  109.     ;经纬度转为网格点后对应的最小最大值
  110.     ;--------------------------------------------------------------------------
  111.     ; Xstart = xy(0,0)         ;经度最小
  112.     ; Xend   = xy(0,1)         ;经度最大
  113.     ; Ystart = xy(1,0)         ;纬度最小
  114.     ; Yend   = xy(1,1)         ;纬度最大

  115.     ;找到dbz最大值的点
  116.     xy = wrf_user_ll_to_ij(a,(/startlon,endlon/),(/startlat,endlat/),res)  ;AB经纬度转换为坐标点
  117.     latlon = wrf_user_ij_to_ll(a,xy(0,:),xy(1,:),res)
  118.     print("----")
  119.     print("Requested  min/max xlon = " + startlon + "/" + endlon)
  120.     print("Actual     min/max xlon = " + xlon(xy(1,0)-1,xy(0,0)-1) + "/" + \
  121.                                          xlon(xy(1,1)-1,xy(0,1)-1))
  122.     print("Calculated min/max xlon = " + latlon(0,0) + "/" + latlon(0,1))

  123.     print("----")
  124.     print("Requested  min/max xlat = " + startlat + "/" + endlat)
  125.     print("Actual     min/max xlat = " + xlat(xy(1,0)-1,xy(0,0)-1) + "/" + \
  126.                                          xlat(xy(1,1)-1,xy(0,1)-1))
  127.     print("Calculated min/max xlat = " + latlon(1,0) + "/" + latlon(1,1))
  128.     print("----")
  129.     xy = xy-1  ;转换为NCL下标

  130.     print("经度格点0 = "+ xy(0,0)) ;经度
  131.     print("经度格点1 = "+ xy(0,1)) ;经度
  132.     print("----")
  133.     print("纬度格点0 = "+ xy(1,0)) ;纬度
  134.     print("纬度格点1 = "+ xy(1,1)) ;纬度
  135.     print("----")

  136.     ;--------------------------------------------------------------------------
  137.     ;本次经纬度转为网格点后对应的最小最大值
  138.     ;--------------------------------------------------------------------------
  139.     ; startlon = xy(0,0)         ;经度最小
  140.     ; startlat = xy(1,0)         ;纬度最大
  141.    
  142.     ; endlon   = xy(0,1)         ;经度最大
  143.     ; endlat   = xy(1,1)         ;纬度最小
  144.    
  145.     ;plot AB垂直剖面的位置

  146.     plane = new(4,float)                 ; a new array
  147.     ;plane  = (/ xy(0,0), xy(1,1)/)              ;A(x,y)点坐标
  148.     ;plane  = (/ xy(0,1), xy(1,0)/)              ;B(x,y)点坐标
  149.     plane  = (/ xy(0,0), xy(1,0), xy(0,1), xy(1,1) /)              ;A(x,y)点到B(x,y)坐标
  150.     alfa   = atan2((endlat-startlat),(endlon-startlon))         ;计算线段AB与x轴方向的夹角
  151.     angle  = 90-alfa*180/PI
  152.     uv     = u*cos(alfa)+v*sin(alfa)                                 ;由于atan2函数正北为0度,因此线段AB与x轴方向的夹角angle和atan2函数出来的alfa之间要换算过
  153.     print("angle of Line AB with X in plot is " +angle+" Degree")

  154.    

  155.     ;opts = False                                   ; start and end point not supplied
  156.     opts = True ; 给定线段AB的话,这个要设置为True才行!!!

  157.     ;--------------------------------------------------------------------------
  158.     ;由于atan2函数正北为0度,因此线段AB的夹角angle和atan2函数出来的alfa之间要换算过
  159.     ;--------------------------------------------------------------------------
  160.     if(ip .eq. 1) then
  161.       ;angle = 90.
  162.       X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) ;This function interpolates 2D model data onto a specified line.
  163.       X_desc = "longitude"
  164.   end if
  165.   if(ip .eq. 2) then
  166.       ;angle = 0.
  167.       X_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
  168.       X_desc = "latitude"
  169.   end if
  170.   if(ip .eq. 3) then
  171.       ;angle = 45.
  172.       X_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
  173.       X_desc = "longitude"
  174.   end if

  175.   ;把数据插值到各个高度上
  176.   dbz_plane = wrf_user_intrp3d(dbz,z,"v",plane,angle,opts)
  177.   u_plane = wrf_user_intrp3d(uv,z,"v",plane,angle,opts)
  178.   w_plane = wrf_user_intrp3d(w,z,"v",plane,angle,opts)
  179.   tc_plane  = wrf_user_intrp3d(tc,z,"v",plane,angle,opts)
  180.   ter_plane = wrf_user_intrp2d(ter,plane,angle,opts)
  181.   mdbz_plane = wrf_user_intrp2d(mdbz,plane,angle,opts)
  182.   p_plane  = wrf_user_intrp3d( p,z,"v",plane,angle,opts)
  183.   rh_plane  = wrf_user_intrp3d(rh,z,"v",plane,angle,opts)
  184.   
  185.   printVarSummary(tc_plane)
  186.   printVarSummary(rh_plane)
  187.   printVarSummary(p_plane)

  188.   print("Max terrain height in plot is " + max(ter_plane)+"m")


  189.   ;*计算假相当位温
  190.   pressure_levels = (/ 1000.,950.,900.,850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,\
  191.                      350.,300.,250.,200.,150.,100./)   ; pressure levels to plot
  192.   nlevels         = dimsizes(pressure_levels)     ; number of pressure levels

  193.   tk = tc_plane+273.166 ; W8 ]: G. [% v
  194.   printVarSummary(tk)
  195.   printVarSummary(pressure_levels)
  196.   ;按指定匹配维扩展到与指定变量同等大小
  197.   ;假定数组q的维数 为 nt x ny x nx x nz  ,而数组dz 为一维数组,长度为nz
  198.   ;按 q 的第3维dz扩展 到 q 的大小
  199.   ;Assume ntim1=100 and ntim2=1, and that T1 is a 4D array with dimensions ntim1 x klev x nlat x mlon and dp is an array with dimensions ntim2 x klev. Note that ntim2 is a degenerate dimension.
  200.   ;In NCL V6.3.0 and later, if you want to conform dp to the same size of T1:
  201.   ;dpC = conform_dims(dimsizes(T), dp, (/0,1/))  ; => (ntim1,klev,nlat,mlon)

  202.   P = p_plane

  203.   a1 = where(tk .gt. 263.0, 0.622*6.11*exp(17.26*(tk-273.16)/(tk-35.86)), \
  204.              0.622*6.11*exp(21.87*(tk-273.16)/(tk-7.66)))   ;
  205.   b1= where(tk .gt. 263.0, P-0.278*exp(17.26*(tk-273.16)/(tk-35.86)),\
  206.             P-0.278*exp(21.87*(tk-273.16)/(tk-7.66)))   ;
  207.   qs1=a1/b1
  208.   q1=qs1*rh_plane
  209.   e1=P*q1/100./(0.62197+q1/100.0)

  210.   tk1=55.0+2840.0/(3.5*log(tk)-log(e1)-4.805)

  211.   pot1=tk*(1000/P)^(0.2854*(1.0-0.28*q1/100.0))
  212.   ept1=pot1*exp(((3376./tk1)-2.54)*q1/100.0*(1.0+0.81*q1/100.0)) ;
  213.   ept1@description = "0se"
  214.   ept1@units = "K"  
  215.   copy_VarCoords(tc_plane,ept1) ; assign coordinates
  216.   printVarSummary(ept1)



  217.   ; Find the index where 6km is - only need to do this once
  218.   if ( FirstTime ) then
  219.       zz = wrf_user_intrp3d(z,z,"v",plane,angle,opts)
  220.       b = ind(zz(:,0) .gt. zmax*1000. )
  221.       zmax_pos = b(0) - 1
  222.       if ( abs(zz(zmax_pos,0)-zmax*1000.) .lt. abs(zz(zmax_pos+1,0)-zmax*1000.) ) then
  223.         zspan = b(0) - 1
  224.     else
  225.         zspan = b(0)
  226.     end if
  227.     delete(zz)
  228.     delete(b)
  229.     FirstTime = False
  230. end if

  231. ; X-axis lables
  232. dimsX = dimsizes(X_plane)
  233. xmin  = X_plane(0)
  234. xmax  = X_plane(dimsX(0)-1)
  235. xspan = dimsX(0)-1
  236. nx    = floattoint( (xmax-xmin)/2 + 1)
  237. ;nx    = 10
  238. print("xmin and xmax in plot is " + xmin+" and "+xmax)
  239. print("xspan and nx in plot is " + xspan+" and "+nx)
  240. ;---------------------------------------------------------------

  241. ; Options for XY Plots
  242. opts_xy                         = res
  243. opts_xy@tiXAxisString           = X_desc
  244. opts_xy@tiYAxisString           = "Height (km)"
  245. opts_xy@cnMissingValPerimOn     = True
  246. opts_xy@cnMissingValFillColor   = 0
  247. opts_xy@cnMissingValFillPattern = 11
  248. opts_xy@tmXTOn                  = False
  249. opts_xy@tmYROn                  = False
  250. opts_xy@tmXBMode                = "Explicit"


  251. ;;-------------------想要绘制的数据的经纬度范围-------------------------------------------------
  252. minlat = startlat
  253. maxlat = endlat
  254. minlon = startlon
  255. maxlon = endlon
  256. loc = wrf_user_ll_to_ij(a,(/minlon,maxlon/),(/minlat,maxlat/),res)  ;查找目标经纬度对应的格点
  257. loc = loc-1   ;转换为NCL下标
  258. ; Xstart = loc(0,0)         ;经度最小
  259. ; Xend   = loc(0,1)         ;经度最大
  260. ; Ystart = loc(1,0)         ;纬度最小
  261. ; Yend   = loc(1,1)         ;纬度最大

  262. ; X-axis lables
  263. xmin  = minlon
  264. xmax  = maxlon
  265. nx    = 10
  266. ;;-----------------------------------------------------------------------------------------------



  267. opts_xy@tmXBMode                = "Explicit"
  268. opts_xy@tmXBValues              = fspan(0,xspan,nx)                    ; Create tick marks
  269. opts_xy@tmXBLabels              = sprintf("%.2f",fspan(xmin,xmax,nx))  ; Create labels
  270. ;opts_xy@tmXBLabels              = (/"29.60N,119.20E","29.40N,119.40E","29.20N,119.60E","29.00N,119.80E","28.80N,120.00E"/)  ; Create labels
  271. opts_xy@tmXBLabelFontHeightF    = 0.010
  272. opts_xy@tmYLMode                = "Explicit"
  273. opts_xy@tmYLValues              = fspan(0,zspan,nz)                    ; Create tick marks
  274. opts_xy@tmYLLabels              = sprintf("%.1f",fspan(zmin,zmax,nz))  ; Create labels
  275. opts_xy@tiXAxisFontHeightF      = 0.020
  276. opts_xy@tiYAxisFontHeightF      = 0.020
  277. opts_xy@tmXBMajorLengthF        = 0.02
  278. opts_xy@tmYLMajorLengthF        = 0.02
  279. opts_xy@tmYLLabelFontHeightF    = 0.015


  280. ; Plotting options for dBz
  281. opts_dbz = opts_xy                        
  282. opts_dbz@cnFillOn = True  
  283. ;opts_dbz@ContourParameters = (/ 5., 65., 5./)
  284. opts_dbz@cnLevelSelectionMode= "ManualLevels"
  285. opts_dbz@cnMinLevelValF      = 0.
  286. opts_dbz@cnMaxLevelValF      = 65.
  287. opts_dbz@cnLevelSpacingF     = 5.0   ;色标间隔  


  288. ; Plotting options for mdBz
  289. opts_mdbz = opts_xy                        
  290. opts_mdbz@cnFillOn = True  
  291. ;opts_mdbz@ContourParameters = (/ 5., 65., 5./)
  292. opts_mdbz@cnLevelSelectionMode= "ManualLevels"
  293. opts_mdbz@cnMinLevelValF      = 0.
  294. opts_mdbz@cnMaxLevelValF      = 65.
  295. opts_mdbz@cnLevelSpacingF     = 5.0   ;色标间隔  


  296. ; Plotting options for terrain
  297. opts_ter = opts_xy
  298. opts_ter@gsnYRefLine = 0.0
  299. opts_ter@gsnAboveYRefLineColor = "black"
  300. opts_ter@gsnDraw = False
  301. opts_ter@gsnFrame = False
  302. opts_ter@trYMaxF = zmax*1000 ;地形数据为m,因此这里要放大倍数为km






  303. vcres = opts_xy
  304. ;vcres@vcGlyphStyle       = "WindBarb"
  305. vcres@vcGlyphStyle       = "CurlyVector";"LineArrow";   
  306. ;矢量设置包括: 
  307. ;“LineArrow”: 多段线和箭头,默认  
  308. ;“FillArrow”: 加填的多边形和箭头  
  309. ;“WindBarb”: 使用天气图上的标准风标图示  
  310. ;“CurlyVector”:卷曲矢量
  311. vcres@vcPositionMode     = "ArrowTail"
  312. vcres@vcRefAnnoOn        = True
  313. vcres@vcMinDistanceF     = 0.03
  314. vcres@vcRefMagnitudeF    = 8.   ;单位长度
  315. vcres@vcLineArrowThicknessF = 3.0
  316. vcres@vcMapDirection        = False
  317. vcres@vcFillArrowsOn     = True
  318. vcres@vcRefLengthF       = 0.015 ;矢量长度
  319. vcres@vcFillArrowHeadXF  = 0.2
  320. vcres@vcFillArrowHeadYF  = 0.2
  321. vcres@vcFillArrowHeadInteriorXF = 0.5


  322. ;;u ,v for streamline
  323. uvres                    = opts_xy                     ; plot mods desired
  324. uvres@tiMainString       = "Streamlines"            ; title
  325. uvres@stArrowLengthF     = 0.004                    ; size of the arrows.
  326. uvres@stMinArrowSpacingF = 0.004                    ; arrow spacing.
  327. uvres@stArrowStride      = 1                        ; arrows start every third
  328. uvres@stArrowLengthF     = 0.008                ; default is dynamic
  329. uvres@stLengthCheckCount = 15                   ; default is 35
  330. uvres@stLineStartStride  = 3.5                    ; default is 2 流线的密度
  331. uvres@stMinArrowSpacingF = 0.035                ; default is 0.0            
  332. uvres@stStepSizeF        = 0.001                ; default is dynamic
  333. uvres@stLineThicknessF   = 2.0           ; changes the line thickness



  334. ;;plot for tc温度
  335. opts_tc = opts_xy                        
  336. opts_tc@cnFillOn = False
  337. opts_tc@cnLineColor="Black"
  338. opts_tc@cnLineThicknessF = 2.0
  339. ;opts_tc@cnLineLabelFontColor="Black"
  340. opts_tc@cnInfoLabelOrthogonalPosF=0.00
  341. opts_tc@ContourParameters=(/ 2.5 /)
  342. opts_tc@gsnContourPosLineDashPattern=1
  343. opts_tc@cnLineDashPattern=1


  344. ;;plot for 0se相当位温
  345. opts_0se = opts_xy                        
  346. opts_0se@cnFillOn = False
  347. opts_0se@cnLineColor="Blue"
  348. opts_0se@cnLineThicknessF = 4.0
  349. opts_tc@cnLineLabelFontColor="Blue"
  350. opts_0se@cnInfoLabelOrthogonalPosF=0.00
  351. opts_0se@ContourParameters=(/ 4.0 /)
  352. opts_0se@gsnContourPosLineDashPattern=0
  353. opts_0se@cnLineDashPattern=0



  354. ; Get the contour info for the rh and temp
  355. contour_dbz = wrf_contour(a,wks,dbz_plane(0:zmax_pos,:),opts_dbz)  ;;loc(0,0):loc(0,1)为想要绘制的经度范围
  356. contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter)
  357. ;vector = gsn_streamline(wks,u_plane(0:zmax_pos,:),w_plane(0:zmax_pos,:)*10.,uvres)
  358. vector = gsn_vector(wks,u_plane(0:zmax_pos,:),w_plane(0:zmax_pos,:)*10.,vcres)
  359. contour_tc  = wrf_contour(a,wks,tc_plane(0:zmax_pos,:),opts_tc)  ;;loc(0,0):loc(0,1)为想要绘制的经度范围
  360. contour_0se  = wrf_contour(a,wks,ept1(0:zmax_pos,:),opts_0se)  ;;loc(0,0):loc(0,1)为想要绘制的经度范围
  361. ; contour_dbz = wrf_contour(a,wks,dbz_plane(0:zmax_pos,loc(0,0):loc(0,1)),opts_dbz)  ;;loc(0,0):loc(0,1)为想要绘制的经度范围
  362. ; contour_ter = gsn_csm_xy(wks,X_plane,ter_plane(loc(0,0):loc(0,1)),opts_ter)
  363. ; vector = gsn_streamline(wks,u_plane(0:zmax_pos,loc(0,0):loc(0,1)),w_plane(0:zmax_pos,loc(0,0):loc(0,1))*10.,uvres)

  364. ; MAKE PLOTS         

  365. if (FirstTimeMap) then

  366.   ;lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
  367.   ;lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts)

  368.   pltres = True
  369.   pltres@FramePlot = False
  370.   ;mpres = True

  371.   optsM = res
  372.   optsM@NoHeaderFooter = True
  373.   optsM@cnFillOn = True
  374.   optsM@lbTitleOn = True;False
  375.   optsM@cnLevelSelectionMode= "ManualLevels"
  376.   optsM@cnMinLevelValF      = 0.
  377.   optsM@cnMaxLevelValF      = 65.
  378.   optsM@cnLevelSpacingF     = 5.0   ;色标间隔
  379.   optsM@pmLabelBarOrthogonalPosF = -0.12 ;!!!!!!!!!!!!!色标与图的距离  !!!!!!!!!!!!!


  380.   ;设置mdbz绘图区域
  381.   minlat=28.0
  382.   minlon=118.5
  383.   maxlat=30.5
  384.   maxlon=122.0
  385.   loc = wrf_user_ll_to_ij(a,(/minlon,maxlon/),(/minlat,maxlat/),res)  ;查找目标经纬度对应的格点
  386.   contour  = wrf_contour(a,wks,mdbz(loc(1,0):loc(1,1),loc(0,0):loc(0,1)),optsM)

  387.   ; Plotting options for Wind Vectors                 
  388.   opts = res         
  389.   opts@FieldTitle = "Wind 10 m"       ; overwrite Field Title
  390.   opts@NumVectors = 47           ; density of wind barbs
  391.   opts@vcWindBarbScaleFactorF  = 2.5                ; 风杆符合国内习惯(4m/s一个长杆)
  392.   opts@vcRefMagnitudeF         = 8.                ; make vectors larger
  393.   opts@vcRefLengthF            = 0.040              ; ref vec length
  394.   opts@vcGlyphStyle            = "LineArrow"         ; select wind barbs
  395.   opts@vcMinDistanceF          = 0.025              ; thin out windbarbs
  396.   opts@vcLineArrowThicknessF =1.5
  397.   opts@vcRefAnnoOn        = True  ;风速单位标注
  398.   ;;绘制水平风场
  399.   vector_uv = wrf_vector(a,wks,u10(loc(1,0):loc(1,1),loc(0,0):loc(0,1)),v10(loc(1,0):loc(1,1),loc(0,0):loc(0,1)),opts)
  400.   delete(opts)

  401.   ;设置绘图区域
  402.   mpres@mpLeftCornerLatF  = minlat
  403.   mpres@mpRightCornerLatF = maxlat
  404.   mpres@mpLeftCornerLonF  = minlon
  405.   mpres@mpRightCornerLonF = maxlon
  406.   ;plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
  407.   plot = wrf_map_overlays(a,wks,(/contour,vector_uv/),pltres,mpres)

  408.   ;绘制垂直剖面的直线AB
  409.   lnres = True
  410.   lnres@gsLineThicknessF = 6.0
  411.   lnres@gsLineColor = "Black"
  412.   ; do ii = 0,dimsX(0)-2
  413.   ;   gsn_polyline(wks,plot,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres)
  414.   ; end do
  415.   gsn_polyline(wks,plot,(/startlon,endlon/),(/startlat,endlat/),lnres)

  416.   frame(wks)
  417.   ;delete(lon_plane)
  418.   ;delete(lat_plane)
  419.   pltres@FramePlot = True
  420. end if

  421. ;plot = wrf_overlays(a,wks,(/contour_dbz,contour_tc,vector,contour_ter/),mpres)    ; plot x-section
  422. plot = wrf_overlays(a,wks,(/contour_dbz,vector,contour_ter,contour_0se/),mpres)    ; plot x-section

  423. ; Delete options and fields, so we don't have carry over
  424. delete(opts_xy)
  425. delete(opts_dbz)
  426. delete(vcres)
  427. delete(X_plane)
  428. delete(ter_plane)
  429. delete(mdbz_plane)


  430. end do  ; make next cross section

  431. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  432. FirstTimeMap = False
  433. end do        ; END OF TIME LOOP

  434. end
复制代码

密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2019-4-30 19:41:51 | 显示全部楼层
_﹎__╱/●勄 发表于 2019-1-9 09:43
楼主,请教一下,我想把剖面的纬度缩小一点也就是X轴范围变小,应该怎么办,用哪个命令?

; X-axis lables
xmin  = minlon
xmax  = maxlon
nx    = 10

opts_xy@tmXBMode                = "Explicit"
opts_xy@tmXBValues              = fspan(0,xspan,nx)                    ; Create tick marks
opts_xy@tmXBLabels              = sprintf("%.2f",fspan(xmin,xmax,nx))  ; Create labels
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-5-8 19:05:14 | 显示全部楼层
小其其格 发表于 2019-4-30 19:37
用以下代码绘制试试看,我的出图没问题:

好!感谢您!我试一下看看!这几天在复习考试没有及时回复
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-5-8 19:08:06 | 显示全部楼层
小其其格 发表于 2019-4-30 19:37
用以下代码绘制试试看,我的出图没问题:

好!感谢您!我试一下看看!这几天在复习考试没有及时回复~
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-5-9 10:33:01 | 显示全部楼层
本帖最后由 雨萌萌 于 2019-5-9 10:34 编辑
小其其格 发表于 2019-4-30 19:41
; X-axis lables
xmin  = minlon
xmax  = maxlon

我试了,画出来还是窜位的~哭
plt_UW_mdBz_08.ab3.ncl.000003.png
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-5-14 15:13:48 | 显示全部楼层
由于atan2函数正北为0度,因此线段AB的夹角angle和atan2函数出来的alfa之间要换算过

这句话不是很理解,如果用grads做剖面这个角度也是要换算吗,我看其他人没有换算过啊,好像是直接用的atan2的值
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-5-22 21:16:34 | 显示全部楼层
请问以下怎么设置截面的纬度
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-5-22 21:18:28 | 显示全部楼层
谢谢楼主,感谢感谢
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

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

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

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