Hello Earth Science Community,
I am a beginner Python programmer trying to plot a single time step of 10m surface wind data I grabbed from the ECMWF ERA5 reanalysis single level dataset in Python 3.8.1 on my mac (Mac OS Mojave 10.14.6). I am having some trouble with the density of the wind vectors on my basemap plot, as the resulting plot from my code produces a near black screen. Below is the code I am working on. Just as a note, the plot also displays mean sea level pressure from one time step and I have been able to plot this successfully.
import pygrib
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import shiftgrid
#---------------------------------------------------------------------------------------------
#Helpful links where I found much of this code
#https://confluence.ecmwf.int/display/CKB/How+to+plot+GRIB+files+with+Python+and+matplotlib
#https://earthscience.stackexchange.com/questions/7012/plotting-wind-barbs-in-python
#-----------------------------------------------------------------------------------------
plt.clf()
fig = plt.figure(figsize = (10,10))
#---------------------------------------------------------------------------------------------
#Set lat/lon coordinates for basemap projection (default projection = cylindrical equadistant)
lllon = -170
lllat = 10
urlon = -110
urlat = 61
#set lat/lon coordinates as well as plotting interval for drawing lat/lon as axis of plot
latmin=0
latmax=70
lonmin=-180
lonmax=-100
latinterval=20
loninterval=20
#---------------------------------------------------------------------------------------------
#create basemap and add some geo content to create a good starting point for plotting data
bmap = Basemap(llcrnrlon=lllon, llcrnrlat=lllat, urcrnrlon=urlon, urcrnrlat=urlat, resolution='f',epsg=3311)
bmap.drawcoastlines(linewidth=2,zorder=1)
bmap.drawstates(linewidth=2, zorder=1)
bmap.drawcountries(linewidth=2,zorder=1)
bmap.drawlsmask(land_color='white',ocean_color='white',lakes='false',resolution='f')
bmap.drawparallels(np.arange(latmin, latmax, latinterval), labels= [1,0,0,0],color='k',textcolor='k',linewidth = 2, fontsize=14)
bmap.drawmeridians(np.arange(lonmin, lonmax, loninterval), labels=[0,0,0,1],color='k',textcolor='k',linewidth = 2, fontsize=14)
#---------------------------------------------------------------------------------------------
#Read ERA5 GRIB file and Grab MSLP values
file = 'download.grib'
grbs = pygrib.open(file)
grb_mslp = grbs.select()[2] #MSLP for 2012-03-26 @00UTC
data = grb_mslp.values #MSLP values in Pa
#---------------------------------------------------------------------------------------------
#Set up plotting parameters lons & lats. Shift lons so projection is -180 to 180 (not 0-360)
#Create evenly spaced lons using first to last grid point. int(grb[Ni]) = 1440 (360/0.25)
lons = np.linspace(float(grb_mslp['longitudeOfFirstGridPointInDegrees']), \
float(grb_mslp['longitudeOfLastGridPointInDegrees']), int(grb_mslp['Ni']) )
#Create evenly spaced lats using first to last grid point. int(grb[Nj]) = 721 (180/.25 + 1)
lats = np.linspace(float(grb['latitudeOfFirstGridPointInDegrees']), \
float(grb['latitudeOfLastGridPointInDegrees']), int(grb_mslp['Nj']) )
#Grid shifting of lons. Not exactly sure how function works.
data, lons = shiftgrid(180., data, lons, start=False)
grid_lon, grid_lat = np.meshgrid(lons, lats) #regularly spaced 2D grid. Still not sure how this works.
#---------------------------------------------------------------------------------------------
#Plotting MSLP on basemap
x,y = bmap(grid_lon,grid_lat) #Pass ERA5 lat/lon to basemap
cs = bmap.contour(x,y,data/100,10,colors='r') #Plot MSLP contours. Divide by 100 = mb
plt.clabel(cs,inline=1,inline_spacing=5,fontsize=10,fmt='%1.1f') #Add labels to contours
#---------------------------------------------------------------------------------------------
#Repeat process for other parameters until I figure out a better method
grb_uwind = grbs.select()[0] #U 10m wind component for 2012-03-26 @00UTC
grb_vwind = grbs.select()[1] #V 10m wind component for 2012-03-26 @00UTC
data_uwind = grb_uwind.values #U 10m wind values in m/s
data_vwind = grb_vwind.values #V 10m wind values in m/s
#---------------------------------------------------------------------------------------------
#Plotting wind on basemap
x,y = bmap(grid_lon,grid_lat) #Pass ERA5 lat/lon to basemap
barbs = plt.quiver(x,y,data_uwind,data_vwind) #Plot 10m wind in m/s
The only solution I took a shot at was trying to plot every 50th data point using the notation:
barbs = plt.quiver(x[::50],y[::50],data_uwind[::50],data_vwind[::50]) #Plot every 50th 10m wind point in m/s
However this produced a map with no almost no results. Would anyone have an idea of how I may be approaching this problem incorrectly? Would it be best to try and guess and check various values to plot instead of every 50 points?
I sincerely appreciate any recommendations on how to proceed forward. Thanks for taking the time to read this.
Google Drive Link for Grib file: https://drive.google.com/file/d/1XEuatdp5piyTIiADgZhxFHQ1_WghMwPJ/view?usp=sharing