[ncl-talk] WRF streamlines in ncl

[Zilore Mumba]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20221006/4e4c92b5/attachment.htm>