[ncl-talk] WRF streamlines in ncl
Dennis Shea
shea at ucar.edu
Fri Oct 7 09:15:16 MDT 2022
Oops. Apologies. I had a 'brain freeze'
For some unknown reason, I 'conflated' stream function and streamline.
stream function: See Example wind_3.ncl @ *
http://ncl.ucar.edu/Applications/wind.shtml*
<http://ncl.ucar.edu/Applications/wind.shtml>
On Fri, Oct 7, 2022 at 2:01 AM Zilore Mumba via ncl-talk <
ncl-talk at mailman.ucar.edu> wrote:
> Rick,
> Thank you once again. wrf_gsn_8.ncl does exactly what I was trying to do.
> I do not know how I missed it. I spent the whole of Thursday 6th Oct.
> trying to make my code work.
> best regards.
>
> On Fri, Oct 7, 2022 at 7:29 AM Zilore Mumba <zmumba at gmail.com> wrote:
>
>> 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
>>>>
>>> _______________________________________________
> 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/9c8ff5e5/attachment.htm>
More information about the ncl-talk
mailing list