[ncl-talk] WRF streamlines in ncl
Zilore Mumba
zmumba at gmail.com
Thu Oct 6 23:29:44 MDT 2022
Thanks Rick, I will try that example.
On Thu, Oct 6, 2022 at 10:37 PM Rick Brownrigg <brownrig at ucar.edu> wrote:
> Hi Zilore,
>
> There's a WRF examples page that's a bit obscurely referenced from the
> main examples pages:
>
> http://ncl.ucar.edu/Applications/wrfgsn.shtml#ex8
>
> Example #8 in particular shows how to set up the coordinates. I am not an
> expert at this, but hopefully there's enough there to get you going.
>
> Rick
>
>
> On Thu, Oct 6, 2022 at 2:04 PM Zilore Mumba via ncl-talk <
> ncl-talk at mailman.ucar.edu> wrote:
>
>>
>> I am wondering if there is anybody drawing streamlines with WRF output.
>> My search has yielded nothing and my attempt has been disastrous.
>> I am having two problems:
>> 1. how do I attach lon/lat units to the u, v components
>> (check_for_y_lat_coord: Warning: Data either does not contain
>> a valid latitude coordinate array or doesn't contain one at all...
>> errors).
>> See code pasted below.
>> 2. The streamlines are drawn, but in the wrong place over a global map.
>> If there is anyone who is able to draw streamlines, I will appreciate to
>> learn from them.
>>
>> my doe so far, below
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>>
>> begin
>> ;
>> ; Read in data.
>> ;;;;;;;;;;;;;;
>> ;
>> f =
>> addfile("/home/zmumba/DA/OUTPUT/2022092800/noda/wrfout_d01_2022-09-28_00:00:
>> 00.nc","r")
>>
>> ; What times and how many time steps are in the data set?
>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>> times = wrf_user_getvar(f,"times",-1) ; get all times in the file
>> ntimes = dimsizes(times) ; number of times in the file
>>
>> ; The specific pressure levels that we want the data interpolated to.
>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>> pressure_levels = (/ 850., 700., 500., 300./) ; pressure levels to
>> plot
>> nlevels = dimsizes(pressure_levels) ; number of pressure
>> levels
>>
>>
>> tc = wrf_user_getvar(f,"tc",-1) ; T in C
>> u = wrf_user_getvar(f,"ua",-1) ; u averaged to mass points
>> v = wrf_user_getvar(f,"va",-1) ; v averaged to mass points
>> p = wrf_user_getvar(f, "pressure",-1) ; pressure is our vertical
>> coordinate
>> z = wrf_user_getvar(f, "z",-1) ; grid point height
>> rh = wrf_user_getvar(f,"rh",-1) ; relative humidity
>>
>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>>
>> ;lat2d = wrf_user_getvar(f,"XLAT",0)
>> ;lon2d = wrf_user_getvar(f,"XLONG",0)
>>
>> lat2d = f->XLAT(0,:,:)
>> lon2d = f->XLONG(0,:,:)
>>
>> lat = lat2d(:,0) ; create classic 1D coordinate arrays
>> lon = lon2d(0,:)
>>
>> ;u&lat at units = "degrees_north"
>> ;u&lon at units = "degrees_east"
>> ;v&lat at units = "degrees_north"
>> ;v&lon at units = "degrees_east"
>>
>> ;lat at units= "degrees_north"
>> ;lon at units= "degrees_east"
>> ;lat!0 = "lat"
>> ;lon!0 = "lon"
>> ;lat&lat = lat
>> ;lon&lon = lon
>>
>>
>> ;pressure = pressure_levels(level)
>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>> pressure = pressure_levels(0)
>>
>> tc_plane = wrf_user_intrp3d(tc,p,"h",pressure,0.,False)
>> z_plane = wrf_user_intrp3d( z,p,"h",pressure,0.,False)
>> rh_plane = wrf_user_intrp3d(rh,p,"h",pressure,0.,False)
>> u_plane = wrf_user_intrp3d( u,p,"h",pressure,0.,False)
>> v_plane = wrf_user_intrp3d( v,p,"h",pressure,0.,False)
>>
>> spd = (u_plane*u_plane + v_plane*v_plane)^(0.5) ; m/sec
>> spd at description = "Wind Speed"
>> spd at units = "m/s"
>> u_plane = u_plane*1.94386 ; kts
>> v_plane = v_plane*1.94386 ; kts
>> u_plane at units = "kts"
>> v_plane at units = "kts"
>>
>> wks = gsn_open_wks("x11","gsn_csm_streamline_map")
>>
>> res = True
>> res at gsnAddCyclic = False
>> res at gsnMaximize = True
>> res at gsnDraw = False ; don't draw yet
>> res at gsnFrame = False ; don't advance yet
>>
>> res at tiMainString = "WRF Streamlines" ; draw a title
>> res at gsnLeftString = ""
>> res at gsnRightString = ""
>> res at gsnCenterString = " "
>>
>> res at mpProjection = "Mercator"
>>
>> ;---Zoom in on map and plot again
>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>> ;res at cnLineThicknessF = 2 ; make lines thicker, Default:
>> 1
>> ;res at mpMinLatF = min(lat) ; zoom in on lat/lon area
>> ;res at mpMaxLatF = max(lat)
>> ;res at mpMinLonF = min(lon)
>> ;res at mpMaxLonF = max(lon)
>>
>> res at mpDataBaseVersion = "MediumRes" ; better and more map
>> outlines
>> res at mpDataSetName = "Earth..4"
>> res at mpOutlineBoundarySets = "AllBoundaries"
>> res at mpOutlineOn = True
>> res at mpNationalLineThicknessF = 2.5 ; turn on country
>> boundaries
>>
>> res at lbOrientation = "Vertical"
>> res at tiMainOffsetYF = -0.03 ; Move the title down
>>
>> plot = gsn_csm_streamline_map(wks,u_plane,v_plane,res)
>>
>> ;Define shapefiles
>> ;*****************
>> shp_filename =
>> ("$HOME/DA/SCRIPTS/06_Utility_Files/Shapefiles/Lesotho/lso_admbnda_adm1_FAO_MLGCA_2019.shp")
>>
>> ;Set shapefile resources
>> ;***********************
>> shpres = True
>> shpres at gsLineThicknessF = 1 ; increase line thickness
>> shpres at gsLineColor = "Black" ; line
>> colorgsLineThicknessF
>> shpres at gsLineDashPattern = 1
>>
>> shp= gsn_add_shapefile_polylines(wks,plot,shp_filename,shpres)
>>
>> draw(plot)
>> frame(wks)
>> end
>> _______________________________________________
>> ncl-talk mailing list
>> ncl-talk at mailman.ucar.edu
>> List instructions, subscriber options, unsubscribe:
>> https://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20221007/c991f108/attachment.htm>
More information about the ncl-talk
mailing list