<div dir="ltr"><br>I am wondering if there is anybody drawing streamlines with WRF output. My search has yielded nothing and my attempt has been disastrous.<br>I am having two problems:<br>1. how do I attach lon/lat units to the u, v components (check_for_y_lat_coord: Warning: Data either does not contain<br>a valid latitude coordinate array or doesn't contain one at all... errors).<br>See code pasted below.<br>2. The streamlines are drawn, but in the wrong place over a global map.<br><div>If there is anyone who is able to draw streamlines, I will appreciate to learn from them.</div><div><br></div><div>my doe so far, below</div><div>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   <br>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   <br><br>begin<br>;<br>; Read in data.<br>;;;;;;;;;;;;;;<br>;<br>  f = addfile("/home/zmumba/DA/OUTPUT/2022092800/noda/wrfout_d01_2022-09-28_00:00:<a href="http://00.nc">00.nc</a>","r")<br><br>   ; What times and how many time steps are in the data set?<br>   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<br>   times = wrf_user_getvar(f,"times",-1)  ; get all times in the file<br>   ntimes = dimsizes(times)         ; number of times in the file<br><br>   ; The specific pressure levels that we want the data interpolated to.<br>   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<br>   pressure_levels = (/ 850., 700., 500., 300./)   ; pressure levels to plot<br>   nlevels         = dimsizes(pressure_levels)     ; number of pressure levels<br><br><br>    tc = wrf_user_getvar(f,"tc",-1)        ; T in C<br>    u  = wrf_user_getvar(f,"ua",-1)        ; u averaged to mass points<br>    v  = wrf_user_getvar(f,"va",-1)        ; v averaged to mass points<br>    p  = wrf_user_getvar(f, "pressure",-1) ; pressure is our vertical coordinate<br>    z  = wrf_user_getvar(f, "z",-1)        ; grid point height<br>    rh = wrf_user_getvar(f,"rh",-1)        ; relative humidity<br>    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<br><br>     ;lat2d    = wrf_user_getvar(f,"XLAT",0)<br>     ;lon2d    = wrf_user_getvar(f,"XLONG",0)<br><br>     lat2d    = f->XLAT(0,:,:)<br>     lon2d    = f->XLONG(0,:,:)<br><br>     lat      = lat2d(:,0)           ; create classic 1D coordinate arrays<br>     lon      = lon2d(0,:)<br><br>      ;u&lat@units = "degrees_north"<br>      ;u&lon@units = "degrees_east" <br>      ;v&lat@units = "degrees_north"<br>      ;v&lon@units = "degrees_east" <br><br>     ;lat@units= "degrees_north"<br>     ;lon@units= "degrees_east"<br>     ;lat!0    = "lat"<br>     ;lon!0    = "lon"<br>     ;lat&lat  =  lat<br>     ;lon&lon  =  lon<br><br><br>      ;pressure = pressure_levels(level)<br>    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<br>      pressure = pressure_levels(0)<br><br>      tc_plane = wrf_user_intrp3d(tc,p,"h",pressure,0.,False)<br>      z_plane  = wrf_user_intrp3d( z,p,"h",pressure,0.,False)<br>      rh_plane = wrf_user_intrp3d(rh,p,"h",pressure,0.,False)<br>      u_plane  = wrf_user_intrp3d( u,p,"h",pressure,0.,False)<br>      v_plane  = wrf_user_intrp3d( v,p,"h",pressure,0.,False)<br><br>      spd     = (u_plane*u_plane + v_plane*v_plane)^(0.5) ; m/sec<br>      spd@description = "Wind Speed"<br>      spd@units = "m/s"<br>      u_plane = u_plane*1.94386     ; kts<br>      v_plane = v_plane*1.94386     ; kts<br>      u_plane@units = "kts"<br>      v_plane@units = "kts"<br><br>     wks  = gsn_open_wks("x11","gsn_csm_streamline_map")<br><br>     res                  = True<br>     res@gsnAddCyclic     = False<br>     res@gsnMaximize      = True<br>     res@gsnDraw          = False         ; don't draw yet<br>     res@gsnFrame         = False         ; don't advance yet<br><br>     res@tiMainString    = "WRF Streamlines"           ; draw a title<br>     res@gsnLeftString   = ""<br>     res@gsnRightString  = ""<br>     res@gsnCenterString = " "<br><br>     res@mpProjection  = "Mercator"<br><br>    ;---Zoom in on map and plot again<br>       ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<br>     ;res@cnLineThicknessF = 2            ; make lines thicker, Default: 1<br>     ;res@mpMinLatF      = min(lat)   ; zoom in on lat/lon area<br>     ;res@mpMaxLatF      = max(lat)<br>     ;res@mpMinLonF      = min(lon)<br>     ;res@mpMaxLonF      = max(lon)<br><br>     res@mpDataBaseVersion     = "MediumRes"       ; better and more map outlines<br>     res@mpDataSetName         = "Earth..4"<br>     res@mpOutlineBoundarySets = "AllBoundaries"<br>     res@mpOutlineOn           = True<br>     res@mpNationalLineThicknessF    = 2.5          ; turn on country boundaries<br><br>     res@lbOrientation         = "Vertical"<br>     res@tiMainOffsetYF        = -0.03           ; Move the title down<br><br>     plot = gsn_csm_streamline_map(wks,u_plane,v_plane,res)<br><br>    ;Define shapefiles<br>    ;*****************<br>     shp_filename = ("$HOME/DA/SCRIPTS/06_Utility_Files/Shapefiles/Lesotho/lso_admbnda_adm1_FAO_MLGCA_2019.shp")<br><br>    ;Set shapefile resources<br>    ;***********************<br>     shpres                    =  True<br>     shpres@gsLineThicknessF   =  1             ; increase line thickness<br>     shpres@gsLineColor        =  "Black"       ; line colorgsLineThicknessF<br>     shpres@gsLineDashPattern  =  1<br><br>     shp= gsn_add_shapefile_polylines(wks,plot,shp_filename,shpres)<br><br>     draw(plot)<br>     frame(wks)<br>end</div></div>