[ncl-talk] WRF streamlines in ncl
Dennis Shea
shea at ucar.edu
Thu Oct 6 16:06:18 MDT 2022
1] NCL calculates streamlines ONLY for global rectilinear grids using
spherical harmonics.
* https://www.ncl.ucar.edu/Applications/stream.shtml*
<https://www.ncl.ucar.edu/Applications/stream.shtml>
NCL does not provide any functions for creating streamlines for WRF
winds, which, I believe, are on a curvilinear staggered grid.
I am not aware of WRF streamlines being plotted in *ANY* language
(Matlab, Python, IDL, GrADS, ferret, ...).
[2] Perhaps, a WRF post-processing group could help.
[3] Generally, speaking, the following is not correct:
lat2d = f->XLAT(0,:,:)
lon2d = f->XLONG(0,:,:)
lat = lat2d(:,0) ; create classic 1D coordinate arrays
lon = lon2d(0,:)
[4] **IF IT WAS CORTRECT**, you would have to do
lat at units= "degrees_north"
lon at units= "degrees_east"
lat!0 = "lat"
lon!0 = "lon"
lat&lat = lat
lon&lon = lon
[5] Hypothetically, if the grid was {u,v}_plane was 2-dimensional
u_plane!0 = "lat"
u_plane!1 = "lon"
u_plane&lat = lat
u_plane&lon = lon
printVarSummary(uplane)
same for v_plane
[5]
plot = gsn_csm_streamline_map(wks,u_plane,v_plane,res)
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/20221006/648989d8/attachment.htm>
More information about the ncl-talk
mailing list