<p>Hello Earth Science Community,</p>

<p>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.</p>



<pre><code>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 &amp; 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
</code></pre>

<p>The only solution I took a shot at was trying to plot every 50th data point using the notation:</p>



<pre><code>barbs = plt.quiver(x[::50],y[::50],data_uwind[::50],data_vwind[::50]) #Plot every 50th 10m wind point in m/s
</code></pre>

<p>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? </p>

<p>I sincerely appreciate any recommendations on how to proceed forward. Thanks for taking the time to read this.</p>

<p>Google Drive Link for Grib file: <a href="https://drive.google.com/file/d/1XEuatdp5piyTIiADgZhxFHQ1_WghMwPJ/view?usp=sharing" rel="noreferrer">https://drive.google.com/file/d/1XEuatdp5piyTIiADgZhxFHQ1_WghMwPJ/view?usp=sharing</a></p>

CC BY-SA 4.0
2020-01-18T03:57:09.343Z
Thanks for the comment! I just added a Google Drive link with the grib file. Let me know if you have any issues with it.
mariandob
CC BY-SA 4.0
2020-01-19T00:14:23.037Z
I would sincerely appreciate it. I started using basemap because a friend provided me with code to get my plots going, but it would probably serve me well to begin to learn how to make use of the cartopy package.
mariandob
CC BY-SA 4.0
2020-01-19T07:26:15.983Z
Sure thing. If I have questions i'll be sure to reach out in the comments section here. Thank you for your assistance.
mariandob
<meteorology><wind><era><python><visualization>
Plotting ERA5 U/V Wind Data Using Python
mariandob
# Plotting ERA5 U/V Wind Data Using Python

(Asked by: username_0 on 2020-01-18)

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>

**Tags:** meteorology, wind, era, python, visualization

> **Comments:**
> **username_0**: Thanks for the comment! I just added a Google Drive link with the grib file. Let me know if you have any issues with it.  
> **username_0**: I would sincerely appreciate it. I started using basemap because a friend provided me with code to get my plots going, but it would probably serve me well to begin to learn how to make use of the cartopy package.  
> **username_0**: Sure thing. If I have questions i'll be sure to reach out in the comments section here. Thank you for your assistance.  

---

## 0 answers

