Plot and Analyze ERA5 Wind Data

Inspired by: (https://github.com/MarieHofmann/Working-with-Climate-Data-in-Python)

Author: Zuhayr Shahid Ishmam

Print names of NetCDF files in the current working directory

import glob
for file_name in glob.iglob('*.nc', recursive=True):
  print(file_name)

ERA5_Wind_monthly_1_deg_Pres_lvls_2018_2020.nc
ERA5_Wind_monthly_Geopot_Vertvel_1_deg_Pres_lvls_2018_2020.nc

Import required Libraries

# Display the plots in the notebook:
%matplotlib inline

# Some defaults:
#plt.rcParams['figure.figsize'] = (12, 5)  # Default plot size
import numpy as np  # numerical library
np.set_printoptions(threshold=20)  # avoid to print very large arrays on screen
# The commands below are to ignore certain warnings.
import warnings
warnings.filterwarnings('ignore')

import matplotlib.pyplot as plt  # plotting library

import xarray as xr  # netCDF library
import cartopy  # Map projections libary, needed to display world maps
import cartopy.crs as ccrs  # Projections list
import logging

Data Download (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-pressure-levels-monthly-means?tab=overview)

import cdsapi

c = cdsapi.Client()

c.retrieve( 'reanalysis-era5-pressure-levels-monthly-means', { 'product_type': 'monthly_averaged_reanalysis', 'variable': [ 'u_component_of_wind', 'v_component_of_wind', ], 'pressure_level': [ '50', '100', '200', '300', '400', '500', '700', '850', '1000', ], 'year': [ '2018', '2019', '2020', ], 'month': [ '01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12', ], 'grid': [1.0, 1.0], 'time': '00:00', 'format': 'netcdf', }, 'ERA5_Wind_monthly_1_deg_Pres_lvls_2018_2020.nc')

Load the Data

ds = xr.open_dataset('ERA5_Wind_monthly_1_deg_Pres_lvls_2018_2020.nc')
ds

<xarray.Dataset>
Dimensions:    (longitude: 360, latitude: 181, level: 9, time: 36)
Coordinates:
  * longitude  (longitude) float32 0.0 1.0 2.0 3.0 ... 356.0 357.0 358.0 359.0
  * latitude   (latitude) float32 90.0 89.0 88.0 87.0 ... -88.0 -89.0 -90.0
  * level      (level) int32 50 100 200 300 400 500 700 850 1000
  * time       (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2020-12-01
Data variables:
    u          (time, level, latitude, longitude) float32 ...
    v          (time, level, latitude, longitude) float32 ...
Attributes:
    Conventions:  CF-1.6
    history:      2022-08-20 11:32:26 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...

Calculate the averaged u and v windfields at the Lowest level for the month of August

Here, the main dataset is global wind data (u and v component, at different pressure levels, 2018 - 2020 - 1 degee res)

u_Jan_1000 = ds.u.groupby('time.month').mean(dim = 'time').sel(level = 1000, month=8).load()
v_Jan_1000 = ds.v.groupby('time.month').mean(dim = 'time').sel(level = 1000, month=8).load()

Quiver Plots

fig = plt.figure(figsize=(15, 10))
ax = plt.axes(projection=ccrs.PlateCarree())  
pu, pv = u_Jan_1000[::15,::15], v_Jan_1000[::15,::15]  
qv = ax.quiver(pu.longitude, pu.latitude, pu, pv, transform=ccrs.PlateCarree())
ax.coastlines(color='black')
plt.title('Average WindSpeed and Direction - Month of August (2018 - 2020) - Pressure Level - 1000hPa',
         fontsize=14)
ax.gridlines(draw_labels = True);

Plot Average WindSpeed and Direction over the Indian Subcontinent - Month of May 2020 - Pressure Level - 1000hPa

south_asia_u_1000_may_2020 = ds.u.sel(latitude = slice(40,7), longitude = slice(70,100), level = 1000, time = '2020-05-01')
south_asia_v_1000_may_2020 = ds.v.sel(latitude = slice(40,7), longitude = slice(70,100), level = 1000, time = '2020-05-01')

windspeed_south_asia_1000_may_2020 = (south_asia_u_1000_may_2020 ** 2 + south_asia_v_1000_may_2020 ** 2) ** 0.5

fig = plt.figure(figsize=(18, 11))
ax = plt.axes(projection=ccrs.PlateCarree())  


pu, pv = south_asia_u_1000_may_2020[::2,::2], south_asia_v_1000_may_2020[::2,::2]  



windspeed_south_asia_1000_may_2020.plot.contourf(ax=ax, transform=ccrs.PlateCarree(), 
                                                 cmap='hot_r', levels = range(0,8), cbar_kwargs={'label':'Windspeed [m/s]'})




qv = ax.quiver(pu.longitude, pu.latitude, pu, pv, transform=ccrs.PlateCarree(), color = 'b')
#qv = ax.quiver(south_asia_u_1000_may_2020.longitude, south_asia_u_1000_may_2020.latitude, south_asia_u_1000_may_2020, 
               #south_asia_v_1000_may_2020, transform=ccrs.PlateCarree())

ax.coastlines(color='black')
ax.add_feature(cartopy.feature.BORDERS)
ax.gridlines(draw_labels = True);
plt.title('Average WindSpeed and Direction - Month of May 2020 - Pressure Level - 1000hPa', fontsize=18)

Text(0.5, 1.0, 'Average WindSpeed and Direction - Month of May 2020 - Pressure Level - 1000hPa')

Plot Average WindSpeed and Direction over South America - Month of May 2020 - Pressure Level - 1000hPa

south_am_u_1000_may_2020 = ds.u.sel(longitude = slice(270,330), latitude = slice(20,-60), level = 1000, time = '2020-05-01')
south_am_v_1000_may_2020 = ds.v.sel(longitude = slice(270,330), latitude = slice(20,-60), level = 1000, time = '2020-05-01')

windspeed_south_am_1000_may_2020 = (south_am_u_1000_may_2020 ** 2 + south_am_v_1000_may_2020 ** 2) ** 0.5

fig = plt.figure(figsize=(18, 11))
ax = plt.axes(projection=ccrs.PlateCarree())  


pu, pv = south_am_u_1000_may_2020[::5,::5], south_am_v_1000_may_2020[::5,::5]  



windspeed_south_am_1000_may_2020.plot.contourf(ax=ax, transform=ccrs.PlateCarree(), 
                                                 cmap='cool', levels = range(0,10), cbar_kwargs={'label':'Windspeed [m/s]'})




qv = ax.quiver(pu.longitude, pu.latitude, pu, pv, transform=ccrs.PlateCarree(), color = 'b')
#qv = ax.quiver(south_asia_u_1000_may_2020.longitude, south_asia_u_1000_may_2020.latitude, south_asia_u_1000_may_2020, 
               #south_asia_v_1000_may_2020, transform=ccrs.PlateCarree())

ax.coastlines(color='black')
ax.add_feature(cartopy.feature.BORDERS)
ax.gridlines(draw_labels = True);
plt.title('Average WindSpeed and Direction - Month of May 2020 - Pressure Level - 1000hPa', fontsize=18)

Text(0.5, 1.0, 'Average WindSpeed and Direction - Month of May 2020 - Pressure Level - 1000hPa')

Streamline Plot of Wind over the Indian Subcontinent - Month of May 2020 - Pressure Level - 1000hPa

fig = plt.figure('figsize', [18,13])

ax = plt.axes(projection=ccrs.PlateCarree())

strm = ax.streamplot(south_asia_u_1000_may_2020.longitude, south_asia_u_1000_may_2020.latitude, 
              south_asia_u_1000_may_2020.values, south_asia_v_1000_may_2020.values, transform=ccrs.PlateCarree(), 
              linewidth = 5, arrowsize = 3, density=10, color=windspeed_south_asia_1000_may_2020.values)

ax.coastlines(color='black')
ax.add_feature(cartopy.feature.BORDERS)
ax.gridlines(draw_labels = True);
bar = plt.colorbar(strm.lines) # add a colorbar defined by the lines of the streamplot (strm.lines)
bar.set_label('Windspeed [m/s]') #give a label to the colorbar

plt.title('Streamline of Wind - Month of May 2020 - Pressure Level - 1000hPa', fontsize=18)

Text(0.5, 1.0, 'Streamline of Wind - Month of May 2020 - Pressure Level - 1000hPa')

Monthly means of the Windspeed in Coastal Bangladesh

u_Cst_BD_monthly = ds.u.groupby('time.month').mean(dim = 'time').sel(level = 1000, 
                                                                    latitude = 22, longitude = 91) 

v_Cst_BD_monthly = ds.v.groupby('time.month').mean(dim = 'time').sel(level = 1000, 
                                                                    latitude = 22, longitude = 91)

ws_Cst_BD_monthly = (u_Cst_BD_monthly**2 + v_Cst_BD_monthly**2)**0.5


fig = plt.figure('figsize', [10,7])

months = ['Jan','Feb', 'Mar','Apr','May','Jun','Jul','Aug','Sep', 'Oct','Nov','Dec']

plt.bar(months, ws_Cst_BD_monthly )


logging.basicConfig(level='INFO')
mlogger = logging.getLogger('matplotlib')
mlogger.setLevel(logging.WARNING)

plt.title("Monthly Average Windspeed (2018-2020) in 22N, 91E Bangladesh", fontsize=18);
plt.ylabel('m s$^{-1}$');
plt.xlabel('');

Display different Pressure Levels of the main dataset

print(ds.level)

<xarray.DataArray 'level' (level: 9)>
array([  50,  100,  200,  300,  400,  500,  700,  850, 1000])
Coordinates:
  * level    (level) int32 50 100 200 300 400 500 700 850 1000
Attributes:
    units:      millibars
    long_name:  pressure_level

Download Geopotential and Vertical Velocity within the same time frame

import cdsapi

c = cdsapi.Client()

c.retrieve( 'reanalysis-era5-pressure-levels-monthly-means', { 'product_type': 'monthly_averaged_reanalysis', 'variable': [ 'geopotential', 'vertical_velocity', ], 'pressure_level': [ '50', '100', '200', '300', '400', '500', '700', '850', '1000', ], 'year': [ '2018', '2019', '2020', ], 'month': [ '01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12', ], 'grid': [1.0, 1.0], 'time': '00:00', 'format': 'netcdf', }, 'ERA5_Wind_monthly_Geopot_Vertvel_1_deg_Pres_lvls_2018_2020.nc')

Load Geopotential and Vertical Velocity Dataset

ds_gp_vv = xr.open_dataset('ERA5_Wind_monthly_Geopot_Vertvel_1_deg_Pres_lvls_2018_2020.nc')
ds_gp_vv

<xarray.Dataset>
Dimensions:    (longitude: 360, latitude: 181, level: 9, time: 36)
Coordinates:
  * longitude  (longitude) float32 0.0 1.0 2.0 3.0 ... 356.0 357.0 358.0 359.0
  * latitude   (latitude) float32 90.0 89.0 88.0 87.0 ... -88.0 -89.0 -90.0
  * level      (level) int32 50 100 200 300 400 500 700 850 1000
  * time       (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2020-12-01
Data variables:
    z          (time, level, latitude, longitude) float32 ...
    w          (time, level, latitude, longitude) float32 ...
Attributes:
    Conventions:  CF-1.6
    history:      2022-08-21 07:43:51 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...

print(ds_gp_vv.level)

<xarray.DataArray 'level' (level: 9)>
array([  50,  100,  200,  300,  400,  500,  700,  850, 1000])
Coordinates:
  * level    (level) int32 50 100 200 300 400 500 700 850 1000
Attributes:
    units:      millibars
    long_name:  pressure_level

Select South America @500 hPa, average over time and plot the Vertical Velocity, w

w_S_Am = ds_gp_vv.w.mean(dim='time').sel(level=500, longitude = slice(270,330), latitude = slice(20,-60))*-1

fig = plt.figure('figsize', [18,13])
ax = plt.axes(projection=ccrs.PlateCarree())
w_S_Am.plot(transform = ccrs.PlateCarree() )
ax.set_title('500 hPa Mean Vertical Winds',fontsize=18)
ax.coastlines(color='black'); 
ax.gridlines(draw_labels = True);
ax.add_feature(cartopy.feature.BORDERS)

<cartopy.mpl.feature_artist.FeatureArtist at 0x2d1c22d6d60>

Average the zonal winds over the African Horn and plot them as a function of the month

w_horn=ds_gp_vv.w.sel(longitude=slice(40, 52), latitude=slice(15, -5))
w_horn = w_horn.groupby('time.month').mean(dim = ['time', 'latitude','longitude'])
w_horn *= -1

plt.figure(figsize=(20,7))
w_horn.T.plot.contourf(levels=30, cbar_kwargs={'label':'m s$^{-1}$'})
plt.ylim([1000,50])
#plt.xlim([1,12])
#plt.yscale('log')
plt.title('Vertical Wind distibution over the North- and South-Eastern Horn (40° to 55°E and 5°S to 15°N)',fontsize=18);

Average the zonal winds over the Indian Subcontinent and plot them as a function of the month

w_ind_bd=ds_gp_vv.w.sel(latitude = slice(40,7), longitude = slice(70,100))
w_ind_bd = w_ind_bd.groupby('time.month').mean(dim = ['time', 'latitude','longitude'])
w_ind_bd *= -1

plt.figure(figsize=(20,7))
w_ind_bd.T.plot.contourf(levels=30, cbar_kwargs={'label':'m s$^{-1}$'})
plt.ylim([1000,50])
#plt.xlim([1,12])
#plt.yscale('log')
plt.title('Monthly Vertical Wind distibution over the Indian Subcontinent (Averaged - 2018-2020)', fontsize=18);

Average over longitude and time and plot the vertical wind as a function of latitude

w_long_time_ave=ds_gp_vv.w.mean(dim= ['longitude', 'time']).sel(latitude = slice(30,-30))
w_long_time_ave *= -1
plt.figure(figsize=(20,7))
w_long_time_ave.plot.contourf(levels=30, cbar_kwargs={'label':'m s$^{-1}$'})
plt.ylim([1000,50])
plt.ylabel('Pressure [hPa]')
plt.xlabel('Latitude [°N]')
#plt.yscale('log')
plt.title('Mean Vertical Wind Distribution from -30 to 30°N');

Zonal plot of the vertical structure of the atmosphere

w_monthly_mean = ds_gp_vv.w.groupby('time.month').mean('time')

#w_monthly_mean.sel(month = 7)

#w_monthly_mean.sel(month = 8).mean(dim='longitude')

plt.figure(figsize=(20,7))
(w_monthly_mean.sel(month = 8).mean(dim='longitude')*-1).plot.contourf();
plt.ylim([1000, 50]); 
plt.xlim([-90,90]);
plt.yscale('log')

plt.title('Average Vertical Wind Distribution August');