[ncl-talk] Wind barb orientation in python plots

[Zilore Mumba]

I am wondering if python does not recognise that wind barbs are directed to
the right in the southern hemisphere and to the left in the northern
hemisphere. This is the same problem in ncl which seems to be carried over
to python.
If it is possible to change this I would appreciate assistance. I reproduce
below a script to plot 850 winds, if it is possible to indicate where to
make any change.

from netCDF4 import Dataset
import cartopy.io.shapereader as shpreader
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature

from wrf import (getvar, interplevel, to_np, latlon_coords, get_cartopy,
 cartopy_xlim, cartopy_ylim)

#Open the NetCDF file
#--------------------
ncfile = Dataset("wrfout_d01_2021-08-19_06:00:00")

#Extract the pressure, geopotential height, and wind variables
#-------------------------------------------------------------
p = getvar(ncfile, "pressure")
z = getvar(ncfile, "z", units="dm")
ua = getvar(ncfile, "ua", units="kt")
va = getvar(ncfile, "va", units="kt")
wspd = getvar(ncfile, "wspd_wdir", units="kts")[0,:]

#Interpolate geopotential height, u, and v winds to 850 hPa
#-----------------------------------------------------------
ht_850 = interplevel(z, p, 850)
u_850 = interplevel(ua, p, 850)
v_850 = interplevel(va, p, 850)
wspd_850 = interplevel(wspd, p, 850)

#Get the lat/lon coordinates
#---------------------------
lats, lons = latlon_coords(ht_850)

#Get the map projection information
#----------------------------------
cart_proj = get_cartopy(ht_850)

#Create the figure
#-----------------
fig = plt.figure(figsize=(12,9))
ax = plt.axes(projection=cart_proj)

#Download and add the states and coastlines
#-----------------------------------------
#states = 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)

fname = './ZMB_ADM1.shp'
adm1_shapes = list(shpreader.Reader(fname).geometries())

#Add the 850 hPa geopotential height contours
#--------------------------------------------
levels = np.arange(150., 155., 1.)
contours = plt.contour(to_np(lons), to_np(lats), to_np(ht_850),
   levels=levels, colors="blue",
   transform=crs.PlateCarree())
plt.clabel(contours, inline=1, fontsize=10, fmt="%i")

#Add the wind speed contours
#---------------------------
levels = [0, 5, 10, 15, 20, 25, 30]
wspd_contours = plt.contourf(to_np(lons), to_np(lats), to_np(wspd_850),
 levels=levels,
 cmap=get_cmap("rainbow"),
 transform=crs.PlateCarree())
plt.colorbar(wspd_contours, ax=ax, orientation="vertical", pad=.05)

#divider = make_axes_locatable(ax)
#ax_cb = divider.new_horizontal(size="5%", pad=0.1, axes_class=plt.Axes)

#plt.add_axes(ax_cb)
#plt.colorbar(wspd_contours, cax=ax_cb, pad=.05)

#Add the 850 hPa wind barbs, only plotting every 5th data point.
#-----------------------------------------------------------------
plt.barbs(to_np(lons[::10,::10]), to_np(lats[::10,::10]),
  to_np(u_850[::10, ::10]), to_np(v_850[::10, ::10]),
  transform=crs.PlateCarree(), length=6)


ax.add_geometries(adm1_shapes, crs.PlateCarree(),
                  edgecolor='black', facecolor='gray', alpha=0.5)

#Set the map bounds
#------------------
ax.set_xlim(cartopy_xlim(ht_850))
ax.set_ylim(cartopy_ylim(ht_850))

ax.gridlines()

plt.title("850 MB Height (dm), Wind Speed (kt), Barbs (kt)")

plt.show()
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20210826/252b7658/attachment.html>