8
$\begingroup$

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

| improve this question | | | | |
$\endgroup$
  • $\begingroup$ can you put your sample grib file somewhere so that it can be downloaded ? $\endgroup$ – gansub Jan 18 at 3:22
  • 1
    $\begingroup$ 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. $\endgroup$ – mariandob Jan 18 at 3:57
  • $\begingroup$ thanks. i will take a look at it later unless someone gets to your question faster $\endgroup$ – gansub Jan 18 at 4:07
  • $\begingroup$ hello i am able to generate wind barbs plot but my code is very different. i use matplotlib and cartopy and no basemap. $\endgroup$ – gansub Jan 18 at 15:00
  • 1
    $\begingroup$ do you want me to show my solution that uses matplotlib and cartopy ? Because basemap has been deprecated for a while now and everyone uses cartopy $\endgroup$ – gansub Jan 18 at 16:30
11
$\begingroup$

First off very few people are using Basemap from Matplotlib these days. From this link matplot basemap

Basemap is deprecated in favor of the Cartopy project. See notes in Cartopy, New Management, and EoL Announcement for more details.

So we are going to use cartopy in addition with matplotlib to plot the grib file that you have provided.

Before I actually show you the code that produces the plot it maybe useful to look at the data structures and their respective shapes.

We basically are looking at four data structures i.e.

latitudes, longitudes and u and v vectors.

file = 'download.grib'
grbs = pygrib.open(file)
for g in grbs:
   lats, lons = g.latlons() # lats - latitude, lons - longitude
   lats, lons = np.array(lats), np.array(lons)
print(lats.shape)
print(lons.shape)
print(data_uwind.shape)
print(data_vwind.shape)

So this prints as

(721,1440)

(721,1440)

(721,1440)

(721,1440)

As one can make out all the shapes are equal. Secondly I am generating the latitude and longitude arrays from the grib file.

When it comes to plotting with cartopy we need to choose a projection. In my case I am choosing Plate Carree. You are welcome to choose another projection depending on your requirements.

So we will now look at the full code using pygrib and matplotlib and cartopy

import pygrib
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import sys

file = 'download.grib'
grbs = pygrib.open(file)
for g in grbs:
    lats, lons = g.latlons() # lats - latitude, lons - longitude
    lats, lons = np.array(lats), np.array(lons)

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

print(lons.shape)
print(lats.shape)
print(data_uwind.shape)
print(data_vwind.shape)

ax1 = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))

ax1.coastlines()

ax1.stock_img()

ax1.set_xticks([0, 60, 120, 180, 240, 300, 359.99], crs=ccrs.PlateCarree())
ax1.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True,
                               number_format='.0f')
lat_formatter = LatitudeFormatter()
ax1.xaxis.set_major_formatter(lon_formatter)
ax1.yaxis.set_major_formatter(lat_formatter)
ax1.barbs(lons[::50,::50],lats[::50,::50],data_uwind[::50,::50],
  data_vwind[::50,::50],transform=ccrs.PlateCarree())

plt.savefig('esse.png')
plt.show()

and finally here is the plot. You can adjust the number of values you want to skip in the barbs function if you want a greater or lesser density by using the Python slicing syntax and skipping either lesser or greater number of elements.

enter image description here

UPDATE

In response to the questions in the comments

If one actually prints what is the type of g in that for loop it turns out

for g in grbs:
    print(type(g))
    lats, lons = g.latlons() # lats - latitude, lons - longitude
    lats, lons = np.array(lats), np.array(lons)

class 'pygrib.gribmessage'

class 'pygrib.gribmessage'

class 'pygrib.gribmessage'

class 'pygrib.gribmessage'

class 'pygrib.gribmessage'

So it is of type gribmessage. The fact it prints five times points to the fact that there are five GRIB Messages in that GRIB file. Each GRIB Message contains both data as well as metadata.

You could if you wish just exit the for loop after just one iteration since the same longitude/latitude is contained in each GRIB Message.

| improve this answer | | | | |
$\endgroup$
  • $\begingroup$ thank you very much for the detailed response. I wanted to ask you a few questions about your solution. Looking at your first section of code where you gather information about the data structures of the grib file, what is the purpose of the for loop? Specifically, what is the first line after the for loop initialization doing? I attached a sketch of my thoughts here: link. I also wanted to ask if the need to shift the longitude data from 0-360 to -180-180 is not applicable here? Thanks again! $\endgroup$ – mariandob Jan 20 at 18:14
  • $\begingroup$ @mariandob i have added extra information. Regarding your second question i do not think so. $\endgroup$ – gansub Jan 21 at 11:44
  • 1
    $\begingroup$ Thank you for all of the helpful information. You have provided me with a great place to learn from. I went ahead and accepted your answer. The only thing I wanted to mention for anyone reading this answer is that when printing the type of "g" from the grib file linked above, using the "for" loop @gansub wrote, 5 "pygrib.messages" should be printed as opposed to 4. $\endgroup$ – mariandob Jan 21 at 12:50
  • 1
    $\begingroup$ @mariandob aah. sometimes i fail to count properly :) $\endgroup$ – gansub Jan 21 at 12:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.