How to Properly Rotate WRF Winds to Earth-Relative Coordinates
Using Python, GEMPAK, and NCL
Updated November 18, 2022 to include Cartopy work-arounds.

by David Ovens, Department of Atmospheric Sciences, University of Washington

Background

WRF is most often used as a regional model, using a Lambert Conformal Conic, Polar Stereographic, or Mercator projection. In all projections, WRF outputs the U and V components of the wind relative to the model grid, not relative to earth coordinates (where U would be west-east and V would be north-south). In the case of a Mercator projection, the model grid does align with earth coordinates and no rotation of the winds is needed. However, for the Lambert Conformal Conic and Polar Stereographic projections, a rotation is necessary whenever one wants to reproject the output or to calculate wind direction.

In my work, I have used a number of tools to rotate WRF winds to earth coordinates. The purpose of this page is to summarize and illustrate the problems, correct some formulas, and give some examples to help future users of WRF find a solution that works for them. I will illustrate solutions using NCL (the NCARGraphics Command Language), Python (using Basemap from matplotlib), GEMPAK, and read_wrf_nc (a WRF utility).

Python Illustration

In this figure below I have attempted to illustrate the problem by modifying a WRF output file's U10 and V10 fields. The V10 field is set to 0 and the U10 field is simply 7 times the value of XLAND. The nco command used to generate these fields is:
ncap2 -O -s 'V10=V10*0.' -s 'U10=XLAND*7.' -s 'V=V*0.' -s 'U=U*0.+10.' wrfout_d01_2016-03-02_15:00:00 wrfout_d1.uvrewrite.nc  
The blue wind barbs are correctly rotated -- they line up with the model's grid points. They are produced by first finding the correct earth-relative components of the wind using the COSALPHA and SINALPHA fields from the WRF output:
  Uearth = U*cosalpha - V*sinalpha
  Vearth = V*cosalpha + U*sinalpha

  NOTE: this is the correct version of these equations. Some sources on the Web incorrectly reverse the signs on the sinalpha terms.
  
Then, those earth-relative components are used along with the rotate_vector() function of matplotlib's basemap to get the new coordinates on the map's x/y grid.
  Urot, Vrot = map.rotate_vector(Uearth,Vearth,lons,lats)
  
map.barbs() can then plot the winds correctly based on whatever projection you have chosen:
  map.barbs(x, y, Urot, Vrot, color='blue')
  
The green wind barbs are the result of simply plotting U10 and V10. Note how they are parallel to the x axis of the plot.
The red wind barbs are the result of plotting the Uearth and Vearth winds without rotating them to the map's projection.
Python with basemapPython with Cartopy

Here is the basemap python code
that was used to generate this plot from this WRF output file:
from netCDF4 import Dataset as NetCDFFile
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
import sys
# allow specification of the wrfout file on the command line.
if (len(sys.argv) > 1):
    ncfile = sys.argv[1]
else:
    ncfile = 'wrfout_d1.uvrewrite.nc'
print 'Using wrfout file: "','\b'+ncfile,'\b"' # use \b to delete space
nc     = NetCDFFile(ncfile, 'r')
#
# CEN_LAT and CEN_LON are specific to the nest and are not
# actually what you want.  If available, you should use
# MOAD_CEN_LAT and STAND_LON
cenlat = nc.getncattr('CEN_LAT')
cenlon = nc.getncattr('CEN_LON')
cenlon = nc.getncattr('STAND_LON')
cenlat = nc.getncattr('MOAD_CEN_LAT')
lat1   = nc.getncattr('TRUELAT1')
lat2   = nc.getncattr('TRUELAT2')
#
# get the actual longitudes, latitudes, and corners
lons = nc.variables['XLONG'][0]
lats = nc.variables['XLAT'][0]
lllon = lons[0,0]
lllat = lats[0,0]
urlon = lons[-1,-1]
urlat = lats[-1,-1]
#
# get the grid-relative wind components and covert from m/s to knots
ur = nc.variables['U10'][0] * 1.94386
vr = nc.variables['V10'][0] * 1.94386
#
# get the WRF local cosine and sines of the map rotation
cosalpha = nc.variables['COSALPHA'][0]
sinalpha = nc.variables['SINALPHA'][0]
#
# get the temperature field to plot a background color contour
t2m = nc.variables['T2'][0]
#
# Make our map which is on a completely different projection.
# To illustrate the problems, go to the upper right corner of the domain,
# or lower right corner in the southern hemisphere, where the distortion
# is the greatest
fig = plt.figure(1, figsize=(12.92, 9.85))
ax = plt.subplot(111)

if urlat >= 0:
    map = Basemap(projection='cyl',llcrnrlat=urlat-5,urcrnrlat=urlat+3,
                llcrnrlon=urlon-10,urcrnrlon=urlon,
                resolution='h')
else:
    map = Basemap(projection='cyl',llcrnrlat=lllat-5,urcrnrlat=lllat+3,
                llcrnrlon=lllon,urcrnrlon=lllon+10,
                resolution='h')
#
x, y = map(lons[:,:], lats[:,:])
cf = map.contourf(x, y, t2m, alpha = 0.4) 
#plt.colorbar(cf, orientation='horizontal', pad=0, aspect=50)
#
# overlay wind barbs doing nothing to show the problem
map.barbs(x, y, ur, vr, color='green',label='Unrotated - wrong')
# rotate winds to earth-relative using the correct formulas
ue = ur * cosalpha - vr * sinalpha
ve = vr * cosalpha + ur * sinalpha
# overlay wind barbs using the earth-relative winds showing we're still wrong
map.barbs(x, y, ue, ve, color='red',label='Rotated to latlon - insufficient')
# now rotate our earth-relative winds to the x/y space of this map projection
urot, vrot = map.rotate_vector(ue,ve,lons,lats)
map.barbs(x, y, urot, vrot, color='blue',label='Rotated to latlon and then to map - CORRECT')
# finish up map drawing
parallels = np.arange(-85.,85,1.)
meridians = np.arange(-180.,180.,1.)
map.drawcoastlines(color = '0.15')
map.drawparallels(parallels,labels=[False,True,True,False])
map.drawmeridians(meridians,labels=[True,False,False,True])
# add a legend.
# NOTE: for wind barb/vector plots, each component of the wind got its own 
# legend entry, so this line
#   plt.legend()
# resulted in duplication.  The next 3 lines handle sifting the legend.
#ax = plt.gca()
handles, labels = ax.get_legend_handles_labels()
legend = plt.legend([handles[0],handles[2],handles[4]], [labels[0],labels[2],labels[4]],loc='upper right')
#
# plot the points of the grid to illustrate where the vectors should be pointing
map.plot(x,y,'bo')
plt.title('How to properly rotate WRF winds to Earth-relative coordinates\nand then reproject them to a new map projection. Using basemap\n')
#plt.show()
plt.savefig('basemapwinds.jpg')
  

Here is the cartopy python code
that was used to generate this plot from this WRF output file:
######################################################################
# Import the needed modules.
#
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER


######################################################################
# The following code reads the example data using the xarray open_dataset
# function and prints the coordinate values that are associated with the
# various variables contained within the file.

ds = xr.open_dataset('wrfout_d1.uvrewrite.nc')
ds.coords

# Grab lat/lon values
lats = ds.XLAT.data[0]
lons = ds.XLONG.data[0]
lllon = lons[0][0]
lllat = lats[0][0]
urlon = lons[-1][-1]
urlat = lats[-1][-1]

# Select and grab data
t2m  = ds['T2'].data[0]
uwnd = ds['U10']
vwnd = ds['V10']
cosalpha = ds['COSALPHA']
sinalpha = ds['SINALPHA']

# properly convert to earth-relative U and V and convert m/s to kt
ue = 1.94386 *(uwnd.data[0]*cosalpha.data[0] - vwnd.data[0]*sinalpha.data[0]) 
ve = 1.94386 *(vwnd.data[0]*cosalpha.data[0] + uwnd.data[0]*sinalpha.data[0]) 

# make grid-relative U and V in knots
ukt = 1.94386 * uwnd.data[0]
vkt = 1.94386 * vwnd.data[0]

# Set up the projection of the data; if lat/lon then PlateCarree is what you want
plotcrs = ccrs.PlateCarree()

# this correction comes from
# https://github.com/SciTools/cartopy/issues/1179
rho = np.pi/180.
u_src_crs = ue / np.cos(lats*rho)
v_src_crs = ve
magnitude = np.sqrt(ue**2 + ve**2)
magn_src_crs = np.sqrt(u_src_crs**2 + v_src_crs**2)
ux = u_src_crs * magnitude / magn_src_crs
vx = v_src_crs * magnitude / magn_src_crs

# Set up the projection that this WRF model actually uses
#datacrs = ccrs.LambertConformal(central_longitude=-121.0,
#                               central_latitude=45.66481,
#                               standard_parallels=(60, 30),
#                               globe=ccrs.Globe(ellipse='WGS84',semimajor_axis=6370000.,semiminor_axis=6370000.))
# The following
#   ux,vx = datacrs.transform_vectors(plotcrs,lons,lats,ue,ve)
# should have worked, but it the U component
# you retrieve needs to be divided by cos(lat) to get the 
# direction correct, but then you mess up the magnitude.
# That's why we had to calculate ux and vx on our own above.
##

# Start the figure and create plot axes with proper projection
fig = plt.figure(1, figsize=(12.92, 9.85))
ax = plt.subplot(111, projection=plotcrs)
if urlat >= 0:
    ax.set_extent([urlon-10, urlon, urlat-5, urlat+3], plotcrs)
else:
    ax.set_extent([lllon, lllon+10, lllat-5, lllat+3], plotcrs)


# Add geopolitical boundaries for map reference
states = cfeature.NaturalEarthFeature(category="cultural", scale="50m",
                             facecolor="none",
                             name="admin_1_states_provinces_shp")
ax.add_feature(states, linewidth=0.5, edgecolor="black")
ax.coastlines('50m', linewidth=0.8)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.add_feature(cfeature.STATES.with_scale('50m'))

# Plot 500-hPa wind barbs in knots, regrid to reduce number of barbs
# Plot 500-hPa wind barbs in knots, regrid to reduce number of barbs
levs = np.arange(234,288,6)
cf = ax.contourf(lons, lats, t2m, levs, alpha = 0.4, transform=plotcrs)
#plt.colorbar(cf, orientation='horizontal', pad=0, aspect=50)
ax.barbs(lons, lats, ukt, vkt, pivot='tip',
         color='green', transform=plotcrs,label='Unrotated - wrong')
ax.barbs(lons, lats, ue, ve, pivot='tip',
         color='red', transform=plotcrs,label='Rotated to latlon - insufficient')
ax.barbs(lons, lats, ux, vx, pivot='tip',
         color='blue', transform=plotcrs,label='Rotated to latlon and then scaled - CORRECT')

# Make some nice titles for the plot (one right, one left)
ax = plt.gca()
handles, labels = ax.get_legend_handles_labels()
legend = plt.legend([handles[0],handles[1],handles[2]], [labels[0],labels[1],labels[2]])
# plot the points of the grid to illustrate where the vectors should be pointing
ax.plot(lons,lats,'bo', transform=plotcrs)
ax.gridlines()
plt.title('How to properly rotate WRF winds to Earth-relative coordinates\nand then reproject them to a new map projection.  Using Cartopy\n')

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=2, color='gray', alpha=0.5, linestyle='--')
gl.top_labels = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.ylabel_style = {'size': 10, 'color': 'black'}
gl.xlabel_style = {'size': 10, 'color': 'black'}


plt.savefig('cartopywinds.jpg')
  

GEMPAK Illustration

In GEMPAK, pass the U and V components into GEMPAK as UREL and VREL and then you can use GEMPAK's built-in VECR() function to tell programs that the winds are grid relative. The figures below plot winds with
    gvect = vecr(urel,vrel)
  or
    gdpfun = vecr(urel,vrel)     
  
NOTE: I was unable to get the proper winds setting PROJ = CED. However, setting PROJ to STR, LCC, MER all work fine.


Here are the settings used to produce this last figure from gdplot2 using this simple10.gem file:

 GDFILE   = simple10.gem
 GDATTIM  = first
 GLEVEL   = 10
 GVCORD   = hagl
 PANEL    = 0
 SKIP     = 0/5
 SCALE    = 0
 GDPFUN   = tmpc@2 !vecr(urel,vrel)
 TYPE     = f ! BM
 CONTUR   = 3
 CINT     = 2
 LINE     = 3/-3/3/2    ! 17/-3/2/2
 FINT     = 0
 FLINE    = 24-8
 HILO     =  
 HLSYM    =  
 CLRBAR   = 32
 WIND     = bm4/0.7
 REFVEC   = 10
 TITLE    = 32//GEMPAK winds vecr(urel,vrel)
 TEXT     = 1/22/hw
 CLEAR    = YES
 GAREA    = 20;-170;65;-85
 IJSKIP   =  
 PROJ     = mer
 MAP      = 32
 MSCALE   = 32;32
 LATLON   = 30/2//5;5
 DEVICE   = XW
 STNPLT   =  
 SATFIL   =  
 RADFIL   =  
 IMCBAR   =  
 LUTFIL   =  
 STREAM   =  
 POSN     = 0
 COLORS   = 1
 MARKER   = 4/17/0.5
 GRDLBL   = 0
 FILTER   = no
  

NCL Illustration

In NCL (the NCARGraphics Command Language) WRF wind rotation is made extremely easy by the use of the uvmet and uvmet10 fields in wrf_user_getvar:
    wrf_user_getvar(filehandle,"uvmet",time) 
    wrf_user_getvar(filehandle,"uvmet10",time) 
    
The figure below uses this function.
Here is the code wrf_wind_reproj.ncl that produced that figure. Note, a typo in this script was corrected Nov 17, 2022.
;----------------------------------------------------------------------
; wrf_wind_reproj.ncl
;----------------------------------------------------------------------
; Concepts illustrated:
;   - using uvmet to rotate WRF grid-relative winds to earth-relative
;     and confirming this by reprojecting WRF output.
;   - Using gsn_csm_contour_map to plot WRF-ARW data
;   - Drawing a WRF lat/lon grid using gsn_coordinates
;----------------------------------------------------------------------
; This example is similar to wrf_gsn_3.ncl, except it zooms in further
; on the map and draws the WRF grid as points.
; Contributed March 2016 by David Ovens at
; http://www.atmos.washington.edu/~ovens/wrfwinds.html
;----------------------------------------------------------------------; 
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
;---Open WRF output file.
  if (.not. isvar("f")) then      ; is f on command line?
      f = "./wrfout_d1.uvrewrite.nc";
  end if
  a = addfile(f,"r")   ; Open a file

;---Read terrain height and lat/lon off file.
  it        = 0     ; first time step
  lats      = wrf_user_getvar(a,"XLAT",it)   ; latitude/longitude
  lons      = wrf_user_getvar(a,"XLONG",it)  ; latitude/longitude

  if(isatt(a,"STAND_LON")) then
     cenlon = a@STAND_LON
  else
     cenlon = a@CEN_LON
  end if
  print("cenlon = "+cenlon)

  u  = wrf_user_getvar(a,"U10",it)        ; u10 already on mass points
  v  = wrf_user_getvar(a,"V10",it)        ; v10 already on mass points
  t2m  = wrf_user_getvar(a,"T2",it)       ; temp at 2 meters in kelvin
  dims = dimsizes(t2m)                     ;
  print("dims = "+dims)
  t2m  = t2m - 273.15                     ; now in celsius
  nd = dimsizes(dims);
  dimU = dims(nd-1) ;
  dimV = dims(nd-2) ;
  print("dimU = "+dimU+" dimV = "+dimV)

  lllat = lats(0,0) ;
  lllon = lons(0,0) ;
  urlat = lats(dimV-1,dimU-1) ;
  urlon = lons(dimV-1,dimU-1) ;
  print("ll latlon ="+lllat+lllon+" ur latlon ="+urlat+urlon)

; NOTE: a lot of people assume that cen_lon is correct, but truly
; you want to use STAND_LON.  This is because you can set your
; longitude that goes exactly vertical in your WRF domain to be off
; the center.  STAND_LON is what you want.
;
; no need to do this:
;  uvm = wrf_uvmet(u,v,hgt@lat2d,hgt@lon2d,-121.,0.716)
; NCL can do the transformation internally:
  uvm = wrf_user_getvar(a,"uvmet10",it) ; get it 
; uc and vc calculated the long way:
;  cosalpha  = wrf_user_getvar(a,"COSALPHA",it)        ; 
;  sinalpha  = wrf_user_getvar(a,"SINALPHA",it)        ; 
;  uc = 1.94386 * (u*cosalpha - v*sinalpha)
;  vc = 1.94386 * (v*cosalpha + u*sinalpha)
; are identical to u,v calculated the easy way:
  u = uvm(0,:,:)*1.94386     ; kts
  v = uvm(1,:,:)*1.94386     ; kts
  u@units = "kts"
  v@units = "kts"
  t2m@lat2d = lats           ;  required for plotting
  t2m@lon2d = lons;          ;  required for plotting
  t2m@description = "2m Temperature"
  t2m@units = "C"
  u@lat2d = lats           ;  required for plotting
  u@lon2d = lons;          ;  required for plotting
  v@lat2d = lats           ;  required for plotting
  v@lon2d = lons;          ;  required for plotting
;  uc@lat2d = lats           ;  required for plotting
;  uc@lon2d = lons;          ;  required for plotting
;  vc@lat2d = lats           ;  required for plotting
;  vc@lon2d = lons;          ;  required for plotting

  wks = gsn_open_wks("png","wrf_gsn")

;---Set some basic plot options
  res               = True

  res@gsnMaximize   = True   ; maximize plot in frame
  res@gsnDraw       = False  ; don't draw
  res@gsnFrame      = False  ; don't advance frame
  res@gsnMaximize   = True   ; use full page
  res@cnInfoLabelOn = False  ; turn off cn info label

  res@mpProjection  = "CylindricalEquidistant"    ; The default

; Usually, when data is placed onto a map, it is TRANSFORMED to the
; specified projection. Since this model is already on a native lambert
; conformal grid, we want to turn OFF the tranformation.
;
  res@tfDoNDCOverlay           = False

;---Zoom in on plot
  res@mpMinLatF     = urlat-5
  res@mpMaxLatF     = urlat+3.
  res@mpMinLonF     = urlon-10
  res@mpMaxLonF     = urlon

  res@gsnAddCyclic             = False            ; regional data 
  res@vcRefMagnitudeF          = 10.0             ; define vector ref mag
  res@vcRefLengthF             = 0.04            ; define length of vec ref
  res@vcGlyphStyle             = "WindBarb"       ; WindBarb or CurlyVector
  res@vcPositionMode           = "ArrowHead"      ; stick tail at value
  res@vcMinDistanceF           = 0.001
  res@gsnScalarContour        = True
  res@cnFillOn                = True
  res@cnFillDrawOrder         = "Draw"
  res@cnFillPalette           = "gui_default"     ; set color map
  res@cnLinesOn               = False             ; turn off contour lines
  res@cnLineLabelsOn          = False             ; turn off contour lines

;---Additional resources desired
  res@pmTickMarkDisplayMode = "Always"   ; nicer tickmarks

  res@mpDataBaseVersion     = "MediumRes"       ; better and more map outlines
  res@mpDataSetName         = "Earth..4"
  res@mpOutlineBoundarySets = "AllBoundaries"
  res@mpOutlineOn           = True
  res@lbLabelBarOn          = False
  res@lbTitleOn             = False

  res@gsnCenterString          = "NCL:  Rotation handled with: ~C~ ~Z75~ uvm = wrf_user_getvar(a,uvmet10,it)"
  res@gsnLeftString          = " ";
  res@gsnRightString          = " ";
  plot = gsn_csm_vector_scalar_map(wks,u,v,t2m,res)
;  plot = gsn_csm_vector_scalar_map(wks,uc,vc,t2m,res)

  mkres                 = True
  mkres@gsMarkerSizeF  = 0.008
  mkres@gsMarkerColor  = "blue"
  gsn_coordinates(wks,plot,t2m,mkres)

end