Plotting NetCDF with Python

This example will examine how to plot time series wind measurements stored as NetCDF datasets, using Python3 (for info on installing Python3 and packages, see our previous blog).

If you would like to use the same NetCDF files, they can be retrieved from ECMWF using their web API. Once you have created an account with ECMWF, you can find the Python retrieval scripts for U and V wind components on my Github here.

Python has a wealth of package options for processing different types of data. For plotting NetCDF files in this example, the following packages are required:

netCDF4

matplotlib

numpy

xarray

Once the necessary packages and dependancies are installed, import them:

# importing of NetCDF files
import netCDF4

# plotting and data visualisations
from matplotlib import pyplot as plt

# n-dimensional arrays
import numpy as np

# high level functions for working with NetCDF files
import xarray as xr

# add colour labels to plot
import matplotlib.patches as mpatches

 

The xarray package provides some very convenient tools for working with NetCDF files. This example uses the xarray open_dataset function, to open two .nc files (for u and v wind components):

# load NetCDF file into variable - using xarray open dataset function
# if using jupyter, make sure NetCDF files are in same directory
ds1 = xr.open_dataset('ERA5_analysis_10U_165_2008.nc')
ds2 = xr.open_dataset('ERA5_analysis_10V_166_2008.nc')

 

Typing the variable name and executing gives you some quick info on the dataset:

ds2

 

You can see the dataset ‘ds2’ contains hourly time-series, grid referenced data. With the variable v10 (describing the v component taken at 10 metre height). ERA5 data is taken at 0.25 degree resolution:

<xarray.Dataset>
Dimensions:    (latitude: 185, longitude: 271, time: 8784)
Coordinates:
  * longitude  (longitude) float32 -22.0 -21.75 -21.5 -21.25 -21.0 -20.75 ...
  * latitude   (latitude) float32 72.5 72.25 72.0 71.75 71.5 71.25 71.0 ...
  * time       (time) datetime64[ns] 2008-01-01 2008-01-01T01:00:00 ...
Data variables:
    v10        (time, latitude, longitude) float32 -0.5332937 -0.42049408 ...
Attributes:
    Conventions:  CF-1.6
    history:      2018-06-22 08:30:01 GMT by grib_to_netcdf-2.8.0: grib_to_ne...

 

From these two datasets you can calculate the total wind speed. It’s important to perform this calculation before doing any time or location extractions, to ensure integrity of the data:

# calculate windspeed
ds12 = ds1
ds12['uv10'] = (ds1['u10']**2 + ds2['v10']**2)**(1/2)

 

As there is one month of hourly data, the time series would be too long to efficiently plot across your x-axis. Luckily, xarray supports ‘group by’ operations. So you can quickly take an day’s hourly average, using the xarray groupby and mean function:

# get daily average
daily_data12 = ds12.groupby('time.hour').mean('time')

 

Next, select your desired grid reference, again using xarray’s sel function, calling the ‘nearest‘ method:

# select a grid location (lon/lat)
dsloc12 = daily_data12.sel(longitude=14,latitude=40,method='nearest')

 

Finally, you can plot the data using matplotlib:

# plot time series
dsloc12['uv10'].plot.line('o-',color='blue',figsize=(15,10))

# change y axis to show zero reference
plt.ylim((-0.5,7.5))

# plot zero reference line
plt.axhline(0, color='grey')

# add colour reference
blue_patch = mpatches.Patch(color='blue', label='10m analysis')
plt.legend(handles=[blue_patch,])

# add titles
plt.ylabel('wind speed m/s')
plt.title('Day Average, 2008 (40°N, 14°E)')

 

Depending on your grid locations, you should get a plot that looks something like this:

 

Repeating the above steps you can easily build multiple plots to compare different data sets:

For more examples on this, see matplotlib documentation.

The code for this example is on my GitHub

by Luke Sanger (WEMC, 2018)

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.