Note
Click here to download the full example code
SyntaxErrorΒΆ
Example script with invalid Python syntax
"""
Plotting Radial WS with Leading Times
====================================
This script plots the 850 hPa radial wind speed with leading times. The `radialAvg.ncl <https://github.com/E3SM-Project/ACME-ECP/blob/master/components/homme/dcmip_tests/dcmip2016_test2_tropical_cyclone/preqx/radialAvg.ncl>`_ needs to be staged to the same directory as the example ncl script below.
"""
loadscript( "radialAvg.ncl")
begin
; Define plot name
pngname="ufs_GFSv16beta_radial_ws_time_plot"
wks=gsn_open_wks("png", pngname)
; Read GFSv16beta vortext tracker results
tcfile="GFSv16beta/fort.69"
delim=","
tclines=asciiread(tcfile, -1, "string")
leadtimestr=tointeger(str_get_field(tclines, 6, delim))
tclatstr=str_get_field(tclines, 7, delim)
tclonstr=str_get_field(tclines, 8, delim)
tcRMWstr=str_get_field(tclines, 20, delim)
tcRMW=tofloat(tcRMWstr)
tcdimsize=dimsizes(leadtimestr)-1
critstr=str_get_field(tclines, 12, delim)
crit=toint(critstr)
count=0
newtcRMW=new((/29/),float)
do i=0,tcdimsize,1
if(crit(i).eq.34) then
newtcRMW(count)=tcRMW(i)
count=count+1
end if
end do
do i=0,tcdimsize,1
if(leadtimestr(i).eq.leadtime) then
tclat=tofloat(str_get_cols(tclatstr(i), 0, 3))*0.1
tclon=tofloat(str_get_cols(tclonstr(i), 0, 4))*0.1
print(tclat)
print(tclon)
end if
end do
psminlat= tclat
psminlon= tclon*(-1)+360
; Use wgrib2 to convert all the GFSPRS* outputs to netcdf format, and read in all the nc files
ncfili=systemfunc("ls GFSv16beta/GFSPRS.GrbF*.nc")
ncfiles=addfiles(ncfili,"r")
UGRD850=ncfiles[:]->UGRD_850mb
VGRD850=ncfiles[:]->VGRD_850mb
WS850=(wind_speed(UGRD850,VGRD850))*1.944
; Make a array for leading time after landfall from f66 to f120
time=(/66,72,78,84,90,96,102,108,114,120/)
dsizes=dimsizes(UGRD850)
;Define a new array for 850 hPa wind speed, /Pressure, latitude, longitude/
verTMP=new((/dsizes(0),dsizes(1),dsizes(2)/),float)
verTMP!0 ="Pressure"
verTMP&Pressure=time ;Trick to replace pressure with leading time data
verTMP&Pressure@units="hPa"
verTMP!1="latitude"
verTMP&latitude=UGRD850&latitude
verTMP!2="longitude"
verTMP&longitude=UGRD850&longitude
verTMP(0,:,:)=(/WS850(0,:,:)/)
verTMP(1,:,:)=(/WS850(1,:,:)/)
verTMP(2,:,:)=(/WS850(2,:,:)/)
verTMP(3,:,:)=(/WS850(3,:,:)/)
verTMP(4,:,:)=(/WS850(4,:,:)/)
verTMP(5,:,:)=(/WS850(5,:,:)/)
verTMP(6,:,:)=(/WS850(6,:,:)/)
verTMP(7,:,:)=(/WS850(7,:,:)/)
verTMP(8,:,:)=(/WS850(8,:,:)/)
verTMP(9,:,:)=(/WS850(9,:,:)/)
; Using the radialAvg3D function from the radialAvg.ncl
outerRad=700.
mergeInnerBins=True
radiaverWS850=radialAvg3D(verTMP(:,:,:),lat,lon,verTMP&Pressure,psminlat,psminlon,outerRad,mergeInnerBins)
radiaverWS850f=tofloat(radiaverWS850)
copy_VarCoords(radiaverWS850, radiaverWS850f)
; Plot the contour field of wind speed at 850hPa
resx=True
resx@gsnDraw = False
resx@gsnFrame=False
resx@cnFillOn = True ; turn on color fill
resx@cnLinesOn = False ; turn lines on/off ; True is default
resx@cnLineLabelsOn = False ; turn line labels on/off ; True is default
resx@cnFillPalette="WhiteBlueGreenYellowRed";"temp_19lev"
resx@cnLevelSelectionMode="ManualLevels"
resx@tmXTOn=False
resx@tmYROn=False
resx@lbOrientation="Vertical"
resx@tiYAxisString ="Forecast Hour"
resx@tiXAxisString="Radius (km)"
radiaverWS850f@units="knots"
radiaverWS850f@long_name="GFSv16beta 850hPa Wind Speed"
resx@cnLevelSelectionMode="ManualLevels"
resx@cnMinLevelValF= 10
resx@cnMaxLevelValF= 60
resx@cnLevelSpacingF= 2
resx@trYMinF=66
resx@trYMaxF=96
resx@tmYLMode="Explicit"
resx@tmYLValues=(/66,72,78,84,90,96/)
resx@tmYLLabels=(/66,72,78,84,90,96/)
plot=gsn_csm_contour(wks, radiaverWS850f(0:5,:), resx)
; Overlay the whiteline of radius of the maximum wind (RMW) to the wind speed contour plot
res=True
res@gsnDraw = False
res@gsnFrame=False
res@xyLineColors = (/"white"/)
res@xyLineThicknesses = (/5.0/)
plotxy=gsn_csm_xy(wks, newtcRMW(10:15), time(0:5),res)
overlay(plot, plotxy)
draw(plot)
frame(wks)
end
Total running time of the script: ( 0 minutes 0.000 seconds)