- 积分
- 8270
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-12-27
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
分别利用ncl和grads对同一个斜线作垂直剖面
出来的图效果感觉ncl的画的不对，请各位大神指点。
----------------------------------
NCL出图效果
ncl1
ncl2
----------------------------------
grads出图效果
---------------------------------------------------------------------------------------------------------------------------------------
ncl绘图脚本
- ; Example script to produce plots for a WRF real-data run,
- ; with the ARW coordinate dynamics option.
- ; Plot data on a cross section
- ; This script will plot data at a set angle through a specified point
- ; This script adds lon/lat info along X-axis
- load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
- load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
- begin
- ;
- ; The WRF ARW input file.
- ; This needs to have a ".nc" appended, so just do it.
- a = addfile("./wrfout_d03_2016-05-05_08_00_00.nc","r")
- ; We generate plots, but what kind do we prefer?
- ; type = "x11"
- ; type = "pdf"
- ; type = "ps"
- ; type = "ncgm"
- type = "png"
- wks = gsn_open_wks(type,"plt_UW_mdBz_08.ncl")
- ;gsn_define_colormap(wks,"WhViBlGrYeOrReWh") ; Overwrite the standard color map
- ;gsn_define_colormap(wks, "Radar2") ;加载自己定义的色标
- gsn_define_colormap(wks, "Radar.mdbz") ;加载自己定义的色标
- ; Set some basic resources
- res = True
- res@MainTitle = "REAL-TIME WRF"
- res@Footer = False
- pltres = True
- mpres = True
- mpres@mpGeophysicalLineColor = "Black"
- mpres@mpNationalLineColor = "Black"
- mpres@mpUSStateLineColor = "Black"
- mpres@mpGridLineColor = "Black"
- mpres@mpLimbLineColor = "Black"
- mpres@mpPerimLineColor = "Black"
- mpres@mpDataSetName = "Earth..4"
- mpres@mpDataBaseVersion = "MediumRes"
- mpres@mpOutlineOn = True
- mpres@mpOutlineSpecifiers = (/"China:states","Taiwan","Zhejiang"/)
- mpres@mpGeophysicalLineThicknessF = 2
- mpres@mpNationalLineThicknessF = 2
- mpres@mpNationalLineColor = "black"
- mpres@mpGeophysicalLineColor = "black"
- mpres@mpProvincialLineColor = "black"
- mpres@mpProvincialLineThicknessF = 2
- mpres@mpFillOn = True
- mpres@mpFillAreaSpecifiers = (/"water", "land" /)
- mpres@mpSpecifiedFillColors = (/"skyblue1","white"/)
- mpres@mpMaskAreaSpecifiers = (/"China:states","Taiwan","Zhejiang"/)
- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
- FirstTime = True
- FirstTimeMap = True
- times = wrf_user_getvar(a,"times",-1) ; get times in the file
- ntimes = dimsizes(times) ; number of times in the file
- mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file
- dNames = getfilevardims(a,"P") ; NCL 6.4.0 and earlier
- nd = dimsizes(mdims) ; 查看mdims为几围数组
- print(dNames + " : " + mdims) ; 查看各个变量维数
- xlat = wrf_user_getvar(a, "XLAT",0)
- xlon = wrf_user_getvar(a, "XLONG",0)
- ter = wrf_user_getvar(a, "HGT",0)
- printVarSummary(xlat)
- ;---------------------------------------------------------------
- do it = 0,ntimes-1 ; TIME LOOP
- print("Working on time: " + times(it) )
- res@TimeLabel = times(it) ; Set Valid time to use on plots
-
- z = wrf_user_getvar(a, "z",it) ; grid point height
- dbz = wrf_user_getvar(a,"dbz",it) ;dBz
- mdbz = wrf_user_getvar(a,"mdbz",it) ;dBz
- printVarSummary(dbz)
- u = wrf_user_getvar(a,"ua",it)
- v = wrf_user_getvar(a,"va",it)
- w = wrf_user_getvar(a,"wa",it)
- if ( FirstTime ) then ; get height info for labels
- zmin = 0.
- zmax = 10. ; We are only interested in the first 6km
- nz = floattoint(zmax + 1)
- end if
- ;---------------------------------------------------------------
- do ip = 1, 1,1 ; we are doing 3 plots；我们这里只绘制2张图
- ; all with the pivot point (plane) in the center of the domain
- ; at angles 0, 45 and 90
- ;
- ; |
- ; angle=0 is | ；垂直剖面直线AB的形状
- ; |
- ;
- ;线段AB的起点和终点经纬度设置
-
- ;08 UTC
- startlon = 119.2
- endlon = 120.0
- startlat = 28.8
- endlat = 29.6
- PI=3.1415926
- ;--------------------------------------------------------------------------
- ;经纬度转为网格点后对应的最小最大值
- ;--------------------------------------------------------------------------
- ; Xstart = xy(0,0) ;经度最小
- ; Xend = xy(0,1) ;经度最大
- ; Ystart = xy(1,0) ;纬度最小
- ; Yend = xy(1,1) ;纬度最大
-
- xy = wrf_user_ll_to_ij(a,(/startlon,endlon/),(/startlat,endlat/),res) ;AB经纬度转换为坐标点
- xy = xy-1 ;转换为NCL下标
- ;plot AB垂直剖面的位置
- ; xlat2d=xlat(xy(1,0):xy(1,1),xy(0,0):xy(0,1))
- ; xlon2d=xlon(xy(1,0):xy(1,1),xy(0,0):xy(0,1))
- ; dimsX=dimsizes(xlat2d)
- ; dimsY=dimsizes(xlon2d)
- ; print(dimsX)
- ; pointX = xy(0,0)+dimsX(1)/2.
- ; pointY = xy(1,0)+dimsX(0)/2.
-
- plane = new(4,float) ; a new array
- plane = (/ xy(0,0),xy(1,1),xy(0,1),xy(1,0)/) ;start x;y & end x;y point
- alfa = atan2((startlat-endlat),(endlon-startlon)) ;计算线段AB与x轴方向的夹角
- angle = 90-alfa*180/PI
- uv = u*cos(alfa)+v*sin(alfa) ;由于atan2函数正北为0度，因此线段AB与x轴方向的夹角angle和atan2函数出来的alfa之间要换算过
- print("angle of Line AB with X in plot is " +angle+" Degree")
-
-
- opts = False ; start and end point not supplied
- ;--------------------------------------------------------------------------
- ;由于atan2函数正北为0度，因此线段AB的夹角angle和atan2函数出来的alfa之间要换算过
- ;--------------------------------------------------------------------------
- if(ip .eq. 1) then
- ;angle = 90.
- X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) ;This function interpolates 2D model data onto a specified line.
- X_desc = "longitude"
- end if
- if(ip .eq. 2) then
- ;angle = 0.
- X_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
- X_desc = "latitude"
- end if
- if(ip .eq. 3) then
- ;angle = 45.
- X_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
- X_desc = "longitude"
- end if
- ;把数据插值到各个高度上
- dbz_plane = wrf_user_intrp3d(dbz,z,"v",plane,angle,opts)
- u_plane = wrf_user_intrp3d(uv,z,"v",plane,angle,opts)
- w_plane = wrf_user_intrp3d(w,z,"v",plane,angle,opts)
- ter_plane = wrf_user_intrp2d(ter,plane,angle,opts)
- mdbz_plane = wrf_user_intrp2d(mdbz,plane,angle,opts)
- print("Max terrain height in plot is " + max(ter_plane)+"m")
-
- ; Find the index where 6km is - only need to do this once
- if ( FirstTime ) then
- zz = wrf_user_intrp3d(z,z,"v",plane,angle,opts)
- b = ind(zz(:,0) .gt. zmax*1000. )
- zmax_pos = b(0) - 1
- if ( abs(zz(zmax_pos,0)-zmax*1000.) .lt. abs(zz(zmax_pos+1,0)-zmax*1000.) ) then
- zspan = b(0) - 1
- else
- zspan = b(0)
- end if
- delete(zz)
- delete(b)
- FirstTime = False
- end if
- ; X-axis lables
- dimsX = dimsizes(X_plane)
- xmin = X_plane(0)
- xmax = X_plane(dimsX(0)-1)
- xspan = dimsX(0)-1
- nx = floattoint( (xmax-xmin)/2 + 1)
- ;nx = 10
- ;print("xspan and nx in plot is " + xmin+" and "+xmax)
- print("xspan and nx in plot is " + xspan+" and "+nx)
- ;---------------------------------------------------------------
-
- ; Options for XY Plots
- opts_xy = res
- opts_xy@tiXAxisString = X_desc
- opts_xy@tiYAxisString = "Height (km)"
- opts_xy@cnMissingValPerimOn = True
- opts_xy@cnMissingValFillColor = 0
- opts_xy@cnMissingValFillPattern = 11
- opts_xy@tmXTOn = False
- opts_xy@tmYROn = False
- opts_xy@tmXBMode = "Explicit"
-
- ;;-------------------想要绘制的数据的经纬度范围-------------------------------------------------
- minlat = startlat
- maxlat = endlat
- minlon = startlon
- maxlon = endlon
- loc = wrf_user_ll_to_ij(a,(/minlon,maxlon/),(/minlat,maxlat/),res) ;查找目标经纬度对应的格点
- loc = loc-1 ;转换为NCL下标
- ; Xstart = loc(0,0) ;经度最小
- ; Xend = loc(0,1) ;经度最大
- ; Ystart = loc(1,0) ;纬度最小
- ; Yend = loc(1,1) ;纬度最大
- xspan = abs(loc(0,1)-loc(0,0))
- xmin = minlon
- xmax = maxlon
- nx = 10
- ;;-----------------------------------------------------------------------------------------------
- opts_xy@tmXBMode = "Explicit"
- opts_xy@tmXBValues = fspan(0,xspan,10) ; Create tick marks
- opts_xy@tmXBLabels = sprintf("%.2f",fspan(xmin,xmax,nx)) ; Create labels
- ;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
- opts_xy@tmXBLabelFontHeightF = 0.010
- opts_xy@tmYLMode = "Explicit"
- opts_xy@tmYLValues = fspan(0,zspan,nz) ; Create tick marks
- opts_xy@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels
- opts_xy@tiXAxisFontHeightF = 0.020
- opts_xy@tiYAxisFontHeightF = 0.020
- opts_xy@tmXBMajorLengthF = 0.02
- opts_xy@tmYLMajorLengthF = 0.02
- opts_xy@tmYLLabelFontHeightF = 0.015
- ; Plotting options for dBz
- opts_dbz = opts_xy
- opts_dbz@cnFillOn = True
- ;opts_dbz@ContourParameters = (/ 5., 65., 5./)
- opts_dbz@cnLevelSelectionMode= "ManualLevels"
- opts_dbz@cnMinLevelValF = 15.
- opts_dbz@cnMaxLevelValF = 65.
- opts_dbz@cnLevelSpacingF = 1.0 ;色标间隔
- ; Plotting options for mdBz
- opts_mdbz = opts_xy
- opts_mdbz@cnFillOn = True
- ;opts_mdbz@ContourParameters = (/ 5., 65., 5./)
- opts_mdbz@cnLevelSelectionMode= "ManualLevels"
- opts_mdbz@cnMinLevelValF = 15.
- opts_mdbz@cnMaxLevelValF = 65.
- opts_mdbz@cnLevelSpacingF = 1.0 ;色标间隔
- ; Plotting options for terrain
- opts_ter = opts_xy
- opts_ter@gsnYRefLine = 0.0
- opts_ter@gsnAboveYRefLineColor = "black"
- opts_ter@gsnDraw = False
- opts_ter@gsnFrame = False
- opts_ter@trYMaxF = zmax*1000 ;地形数据为m，因此这里要放大倍数为km
- vcres = opts_xy
- ;vcres@vcGlyphStyle = "WindBarb"
- vcres@vcGlyphStyle = "LineArrow";"CurlyVector"
- vcres@vcPositionMode = "ArrowTail"
- vcres@vcRefAnnoOn = True
- vcres@vcMinDistanceF = 0.03
- vcres@vcRefMagnitudeF = 8. ;单位长度
- vcres@vcLineArrowThicknessF = 3.0
- vcres@vcMapDirection = False
- vcres@vcFillArrowsOn = True
- vcres@vcRefLengthF = 0.02
- vcres@vcFillArrowHeadXF = 0.2
- vcres@vcFillArrowHeadYF = 0.2
- vcres@vcFillArrowHeadInteriorXF = 0.15
-
- ;;u ,v for streamline
- uvres = opts_xy ; plot mods desired
- uvres@tiMainString = "Streamlines" ; title
- uvres@stArrowLengthF = 0.004 ; size of the arrows.
- uvres@stMinArrowSpacingF = 0.004 ; arrow spacing.
- uvres@stArrowStride = 1 ; arrows start every third
- uvres@stArrowLengthF = 0.008 ; default is dynamic
- uvres@stLengthCheckCount = 15 ; default is 35
- uvres@stLineStartStride = 3.5 ; default is 2 流线的密度
- uvres@stMinArrowSpacingF = 0.035 ; default is 0.0
- uvres@stStepSizeF = 0.001 ; default is dynamic
- uvres@stLineThicknessF = 2.0 ; changes the line thickness
-
-
- ; Get the contour info for the rh and temp
- contour_dbz = wrf_contour(a,wks,dbz_plane(0:zmax_pos,xy(0,0):xy(0,1)),opts_dbz) ;;loc(0,0):loc(0,1)为想要绘制的经度范围
- contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter)
- vector = gsn_streamline(wks,u_plane(0:zmax_pos,xy(0,0):xy(0,1)),w_plane(0:zmax_pos,xy(0,0):xy(0,1))*10.,uvres)
-
- ; MAKE PLOTS
- if (FirstTimeMap) then
- lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
- lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
-
- pltres = True
- pltres@FramePlot = False
- ;mpres = True
- optsM = res
- optsM@NoHeaderFooter = True
- optsM@cnFillOn = True
- optsM@lbTitleOn = False
- optsM@cnLevelSelectionMode= "ManualLevels"
- optsM@cnMinLevelValF = 15.
- optsM@cnMaxLevelValF = 65.
- optsM@cnLevelSpacingF = 1.0 ;色标间隔
- ;设置mdbz绘图区域
- minlat = 27.
- maxlat = 31.5
- minlon = 118.0
- maxlon = 122.5
- loc = wrf_user_ll_to_ij(a,(/minlon,maxlon/),(/minlat,maxlat/),res) ;查找目标经纬度对应的格点
- contour = wrf_contour(a,wks,mdbz(loc(1,0):loc(1,1),loc(0,0):loc(0,1)),optsM)
- ;设置绘图区域
- mpres@mpLeftCornerLatF = minlat
- mpres@mpRightCornerLatF = maxlat
- mpres@mpLeftCornerLonF = minlon
- mpres@mpRightCornerLonF = maxlon
- plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
- ;绘制垂直剖面的直线AB
- lnres = True
- lnres@gsLineThicknessF = 3.0
- lnres@gsLineColor = "Black"
- do ii = 0,dimsX(0)-2
- gsn_polyline(wks,plot,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres)
- end do
- frame(wks)
- delete(lon_plane)
- delete(lat_plane)
- pltres@FramePlot = True
- end if
-
- plot = wrf_overlays(a,wks,(/contour_dbz,vector,contour_ter/),mpres) ; plot x-section
- ; Delete options and fields, so we don't have carry over
- delete(opts_xy)
- delete(opts_dbz)
- delete(vcres)
- delete(X_plane)
- delete(ter_plane)
- delete(mdbz_plane)
-
- end do ; make next cross section
- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
- FirstTimeMap = False
- end do ; END OF TIME LOOP
- end
|
|
窥视卡
雷达卡
发表于 2018-10-26 16:43:44
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主