Python | Tutorial: Introduction to reading netCDF files

Introduction

NetCDF (Network Common Data Format) is an extremely useful way to store, transfer, subset, and otherwise work with gridded data sets.  This format is commonly encountered in the atmospheric sciences, climate, and oceanography communities.  It is much easier to work with than the "grib" data format, as it allows for metadata attributes to be attached directly to the datafile.  Finally, there are numerous stand alone programs (e.g., ncview, nco) and bindings (Fortran, Python, IDL, Matlab, GrADS:see here for my post on netCDF with GrADS) that make it easy to work with NetCDF data regardless of platform, or programming language.  Because NetCDF is so common, and because most of the blog posts on this site require some knowledge on how to work with NetCDF data, I'm going to provide a short tutorial on how to read and work with NetCDF files in Python.  In this tutorial, you will learn how to:

• Open at NetCDF file using "NetCDF4"
• Query the file attributes, dimensions, and variables
• Convert NetCDF integer time to Python "datetime" objects
• Extract monthly temperature anomaly
• Interface with the NOMADS OpenDAP-alt data repository

Important Note: There are two supported packages (that I'm aware of) that can open and translate NetCDF data (NetCDF4 and xarray).  For this tutorial, I will focus on the NetCDF4 package, perhaps I will post an xarray tutorial later on, but if you're already reading and working with NetCDF data using xarray, and you have no intention on changing, then maybe this tutorial isn't for you.

Opening, and Querying a NetCDF file

First off, if you are going to work with the NetCDF data, you will need to have the NetCDF4 package installed and imported.

Once you have that package installed, we'll need to grab some data.  For this tutorial, I'll be using monthly mean air temperature data from the NCEP reanalysis.  To get this data, go to this website and download the "air.mon.mean.nc" file and place it in the same directory as the script you'll be working with.  Once that is done, you're ready to proceed with the tutorial.

The first step is to import the relevant packages:

import numpy as np
import netCDF4 as NC
from matplotlib import pyplot as plt 
from cartopy import crs as ccrs

import datetime as DT

I included the Cartopy import to plot the data on a map, if you don't have this package and you don't want to install it, that's fine you can simply remove that line from the import section, and skip the mapping portion of the tutorial.

Once the packages are imported, set the path to the file, and use the "Dataset" function to open the data file:

## Since file is in current folder, the path is just the file name
ncdata_path='air.mon.mean.nc'
ncdata=NC.Dataset(ncdata_path,'r')
print(ncdata)

This will result in the following output to the screen:

root group (NETCDF4_CLASSIC data model, file format HDF5):
    description: Data from NCEP initialized reanalysis (4x/day).  These are the 0.9950 sigma level values
    platform: Model
    Conventions: COARDS
    NCO: 20121012
    history: Thu May  4 20:11:16 2000: ncrcat -d time,0,623 /Datasets/ncep.reanalysis.derived/surface/air.mon.mean.nc air.mon.mean.nc
Thu May  4 18:11:50 2000: ncrcat -d time,0,622 /Datasets/ncep.reanalysis.derived/surface/air.mon.mean.nc ./surface/air.mon.mean.nc
Mon Jul  5 23:47:18 1999: ncrcat ./air.mon.mean.nc /Datasets/ncep.reanalysis.derived/surface/air.mon.mean.nc /dm/dmwork/nmc.rean.ingest/combinedMMs/surface/air.mon.mean.nc
/home/hoop/crdc/cpreanjuke2farm/cpreanjuke2farm Mon Oct 23 21:04:20 1995 from air.sfc.gauss.85.nc
created 95/03/13 by Hoop (netCDF2.3)
Converted to chunked, deflated non-packed NetCDF4 2014/09
    title: monthly mean air.sig995 from the NCEP Reanalysis
    References: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.derived.html
    dataset_title: NCEP-NCAR Reanalysis 1
    dimensions(sizes): lat(73), lon(144), time(863)
    variables(dimensions): float32 lat(lat), float32 lon(lon), float64 time(time), float32 air(time,lat,lon)
    groups: 

Note that a lot of the information above is probably not relevant to you, but you can see some pieces of information in there that might be useful.  One thing that you notice, is that there printed output has colon (:) separated text, these correspond to the global metadata attributes, and structures attached to the file, and they can be accessed individually using the following syntax:
For example, accessing the file "description"

print(ncdata.description)

Returns:

Data from NCEP initialized reanalysis (4x/day).  These are the 0.9950 sigma level values

Which you'll notice matches the text after the colon separating the word "description"

Though, it can be kind of burdensom to have to print the entire netCDF file, just to see what's in it before you can start extracting data.  One way to findout what attributes the file has without dumping the entire file out to screen is to use the "ncattrs" function which is a shorthand for "NetCDF (nc) + attributes (attrs) - ncattrs" to find out what attributes the file has:

print(ncdata.ncattrs())

Returns:

['description', 'platform', 'Conventions', 'NCO', 'history',
'title', 'References', 'dataset_title']

You can use the "ncgetattr" function to return an attribute value for a given attribute.  Following the example above:

print(ncdata.getncattr('description'))

Returns:

Data from NCEP initialized reanalysis (4x/day).  These are the 0.9950 sigma level values

You may be asking the question: why would I use the "getncattr" function to get an attribute instead of just using the shorthand "ncdata.description"?  The advantage is that the getncattr function allows you to pass string information, where as the shorthand method does not, which makes it more versitile incase you don't know what the attributes are.

Loading NetCDF variables and Dimensions

Dimensions:

The dimensions structure includes the dimensions used in the file.  In geographic datasets, these are typically x,y,z,time.  x and y are often latitude and longitude, but not always (for example, WRF data has the x and y dimensions as "west_east" and "south_north").  z is usually a vertical coordinate (sigma, height, pressure, etc), and time is usually a floating point value referenced to a specific starting time.  The file dimensions can be printed simply:

print(ncdata.dimensions)

Returns:

OrderedDict([('lat', <class 'netCDF4._netCDF4.Dimension'>: name = 'lat', size = 73), 

('lon', <class 'netCDF4._netCDF4.Dimension'>: name = 'lon', size = 144), 

('time', <class 'netCDF4._netCDF4.Dimension'> (unlimited): name = 'time', size = 863)])

What you get is a dictionary of dimensions with names and sizes: in this instance, latitude, longitude, and time.  Now, the key about the dimensions, is that EACH variable will have dimensions that will be some combination of the dimensions within the dimensions dictionary.

Variables:

Simiarly, the variables, can be accessed:

for v in ncdata.variables:
    print(v)

Returns:
lat
lon
time
air

Similar to the dimensions object, the variables object is a dictionary, where the variables can be accessed using keys.  To demonstrate, I'm going to access the "time" variable and convert the time values into Python datetime.
Printing the time variable:

print(ncdata.variables['time'])

Returns:
<class 'netCDF4._netCDF4.Variable'>
float64 time(time)
    long_name: Time
    delta_t: 0000-01-00 00:00:00
    avg_period: 0000-01-00 00:00:00
    prev_avg_period: 0000-00-01 00:00:00
    standard_name: time
    axis: T
    units: hours since 1800-01-01 00:00:0.0
    actual_range: [1297320. 1927008.]
unlimited dimensions: time
current shape = (863,)
filling on, default _FillValue of 9.969209968386869e+36 used

This statement tells us that the time data is stored as a 64bit (double) float, and is hours since 1800-01-01 00 UTC. We also see that it's only dimension is the "time" dimension.  So, from this data, we can use the datetime strptime and timedelta functions to convert the floating point values to a list of datetime objects.  Knowing the string format of the "units" attribute, the starting "strp" time object is set following:

starttime=DT.datetime.strptime(ncdata.variables['time'].units,
        'hours since %Y-%m-%d %H:%M:0.0')

Now, pulling the time data into a Numpy array, and using list iteration, we can build a list of datetime objects:

time_values=ncdata.variables['time'][:]
strptimes=[starttime+DT.timedelta(hours=i) for i in time_values]

Okay, now that we have a list of times, let's use it to get a time index for a specific input time.  For example: March 1997:

choice='1997-03-01'
choice_strp=DT.datetime.strptime(choice,'%Y-%m-%d')
tidx=np.argmin(np.abs(np.array(strptimes)-choice_strp))

print("Index:%i, time:"%tidx, strptimes[tidx])

Returns:

Index:590, time: 1997-03-01 00:00:00

Now, finally, lets plot the data!  Looking back at the variables, data, we see a lat,lon and air variable as well.  Each variable can be printed to find out it's metadata but, for brevity, I'll just tell you that, lon and lat are each 1-dimensional arrays that vary with their corresponding dimensions, and air is air temperature varying in time,lat,lon.

Simply accessing the data for each dimension:

t2=ncdata.variables['air'][tidx,:]
lon=ncdata.variables['lon'][:]
lat=ncdata.variables['lat'][:]

and plotting it:

plt.figure(figsize=(10,7))
plt.pcolormesh(lon,lat,t2,cmap='jet')
plt.show()

Simple air temperature plot

Or, if you want to include the Cartopy mapping:

plt.figure(figsize=(10,7))
ax=plt.subplot(111,projection=ccrs.PlateCarree())
Z=ax.pcolormesh(lon,lat,t2,cmap='jet')
plt.colorbar(Z)
ax.coastlines()
plt.show()

With a map and colorbar.

Finally, since 1997 was a famous En Nino year, we can plot the annual temperature anomaly by passing a tuple of time indices to the time index of the data array and plot the annual mean temperature anomaly.

month_idx=[i for i in range(len(strptimes)) if strptimes[i].strftime('%Y')=='1997']
nino_mean=np.mean(ncdata.variables['air'][tuple(month_idx),:],axis=0)
total_mean=np.mean(ncdata.variables['air'][:],axis=0)

plt.figure(figsize=(10,7))
ax=plt.subplot(111,projection=ccrs.PlateCarree())
Z=ax.pcolormesh(lon,lat,nino_mean-total_mean,cmap='coolwarm',vmin=-2.5,vmax=2.5)
plt.title("Annual Mean Temperature Anomaly: 1997")
plt.colorbar(Z)
ax.coastlines()
plt.show()

Basically, we build a tuple of indices where that matches where the time array indicates 1997 and then average the temperature data for all of those indicies to get the 1997 mean temperature.  We then compute a total mean, which takes the mean along the entire time axis, and difference the two values.  This results in the following image:

1997 ENSO event shown in the 0.995 sigma NCEP reanalysis temperature field

Interfacing with the NOMADS server

One great aspect of using the NetCDF library is that you can access near realtime model data from the NOAA NOMADS data server.  This section will be a short add on, since there isn't a whole lot new to add, but it will show you how to access near real time model data.

First, a simple line of code is required to a) get todays date, and b) define a filepath to the data.  In this example, I'll use the 0.25 degree GFS data, explore the NOMADS server for more options.  Assuming you have imported the packages above:

today=DT.datetime.now().strftime('%Y%m%d')
input_file='http://nomads.ncep.noaa.gov:80/dods/gfs_0p25/gfs%s/gfs_0p25_00z'%today
gfsdata=NC.Dataset(input_file,'r')

Now, you are welcome (encouraged?) to query the file using the same techniques as above, but for brevity, I'll just write down the code to store the lat/lon values into arrays (which similarly, are 1-D arrays cooresponding to x and y).  Note that time is stored in a very similar way as in the NCEP reanalysis file and can be handled in an identical way, but won't be demonstrated here.

lons=gfsdata.variables['lon'][:]
lats=gfsdata.variables['lat'][:] 

Then the data can be extracted and plotted (in this example, 2 meter temperature at time index 6):

t2=gfsdata.variables['tmp2m'][6,:]

plt.figure(figsize=(11,9))
ax=plt.subplot(111,projection=ccrs.PlateCarree())
plt.pcolormesh(lons,lats,t2,cmap='jet')
plt.colorbar()
plt.title(gfsdata.variables['tmp2m'].long_name)
ax.coastlines()
plt.show()

Note that the above example gets the figure title directly from the NetCDF metadata.

Finally, it's always a good practice to close your NetCDF files when you are done using them:

gfsdata.close()

Note that when you close a NetCDF file, any data that you have loaded into a Numpy array is still accessible.  However!  Any variables that are just pointing to the NetCDF object are no longer accesible.  For example:

t2_pointer=gfsdata.variables['tmp2m']
t2_values=gfsdata.variables['tmp2m'][6,:]

gfsdata.close()
 
t2_values ## Still here!
t2_pointer ## GONE!

Wrap up

That's it for the introduction to NetCDF tutorial.  There is a lot left to be discovered and learned when it comes to working with NetCDF data using Python, however I tried to provide a basic overview on how to open, query, and load data from netCDF files with a eye towards more common applications. There is a lot of great documentation out there if you want to learn more, a good place to start is the Unidata Documentation.

So get out there and start working with NetCDF data!

Read More

Tutorial: The difference between self describing and non-self describing NETCDF files, and how to write them

In writing posts for this blog, I try my best to avoid what could be considered "mission creep," or writing coding tutorials that at best only tangentially related to GrADS.  This tutorial is verging upon what could be considered mission creep and has as much to do with teaching "good data" practices as it does with teaching GrADS skills.  Furthermore, there are no GrADS scripting lessons to be learned here, and the only code I provide is written for Python. 

All that said, in this tutorial, I will be showing you how to write NETCDF files that can be read using GrADS using both the self describing 'sdfopen' and non-self describing 'xdfopen' commands. I do this because, you may one day find yourself with some gridded data set that you want to plot using GrADS, however this data set may be in a text file format, or perhaps maybe you need to perform higher level calculations on the data that cannot be performed in GrADS before you display it.  While you could write the data out as a binary data file with a GrADS control file, maybe you are more comfortable with the NETCDF format (I know I certainly am).  Plus there are other benefits to using NETCDF instead of binary; one example is cross compatibility. NETCDF files can be read in with pretty much any data analysis software, so not only can you read your file in GrADS, but you can read it in with Matlab, IDL, R, Python, etc., as well.

What makes a NETCDF file self describing?

It's all about the metadata!  Metadata is often referred to "data about the data."  For example, metadata includes things like:

• The data and time a file was created
• What data is in the file
• Units of variables within the file
• Information about the coordinates (e.g., lat, lon vs. x,y)
• If the data file is for model output it might include information such as 

• Resolution, parameterizations, time step

Metadata is super important, and learning how to include it in any data set you are producing/sharing is a critical skill in a scientific career.  We are lucky in meteorology because there have been solid efforts to standardize metadata practices in the past couple of decades.  Now, NETCDF makes it very easy to provide good metadata, so there is no excuse not to do it.  So if nothing else, I hope to get that point across to you with this tutorial.  To really drive this point home, here is a meme.

 

To make a file self-describing (that is such that GrADS can read the file in with NO accompanying control file required using the 'sfdopen' command) the metadata including in the NETCDF file must adhere to the COARDS data standards: More information on COARDS here.  Personally, I didn't find this website all that useful when figuring out how to make self-describing NETCDF files, but it does provide some nice background.

How to write a self describing NETCDF file (SDFOPEN)

Here I will provide the key metadata required to write a self-describing NETCDF file.  The example, I provide is sub-setting the first 24 hours of GFS temperature data from the 00z December 30th 2015 run.  If you choose to try and reproduce these results, you will need to update the date.  I will us python to read in the data from the GrADS Data Server (GDS) and to write the local NETCDF file.  You could use any program of your choice to do this, so long as it reads and writes NETCDF files.  The python script is included in this tutorial.

So the most important thing to understand when writing your datafile is that your coordinate dimensions and accompanying dimension variables (e.g., lat,lon) adhere to the COARDS guidelines, otherwise GrADS will not open your file.  For example, you might get an error like:

     gadsdf: SDF file has no discernable Y coordinate

The metadata for the non-dimension variables (e.g., temperature, u,v) is more up to the person writing the file.

Assuming you know how to read in the NETCDF file and data from the GDS, in python:

   datafile='http://nomads.ncep.noaa.gov:9090/dods/gfs_0p25/gfs20151230/gfs_0p25_00z'
   cdata=netCDF4.Dataset(datafile,'r')
 
To start writing the file, we will first look at the dimensions and dimension variables.
You must include the following dimensions and variables in your NETCDF file:

• lat
• lon
• time
• lev

 In python you will make these dimensions and variables as such:

   newdata=netCDF4.Dataset(newfile,'w')
   newdata.createDimension('lat', len(lat[lmin:lmax]))
   newdata.createDimension('lon', len(lon[lomin:lomax]))
   newdata.createDimension('lev', len(lev[:]))
   newdata.createDimension('time', 8)

Where lmin,lmax,lomin,lomax are indices that mark your subset boundaries.
Then in addition to creating the dimensions, you must make the variables as well (For example latitude):

   latvar=newdata.createVariable('lat','f8',('lat'))

Now you need to populate latvar with metadata.  To determine what metadata you need, I simply printed out the 'lat' variable from within the GDS datafile as this datafile is self-describing and adheres to the COARDS standards.  If you've seen my other tutorials I often access these files using the 'sdfopen' command.  In doing that, I get the following (again this will look different if you aren't using python):

<type 'netCDF4.Variable'>
float64 lat(lat)
    grads_dim: y
    grads_mapping: linear
    grads_size: 721
    units: degrees_north
    long_name: latitude
    minimum: -90.0
    maximum: 90.0
    resolution: 0.25
unlimited dimensions:
current shape = (721,)

These are the metadata that you need to include for the latitude variable in your NETCDF file.  The dimensions will change, because you are sub setting the data, but things like 'grads_dim' and 'grads_mapping' need to be the same.  In python adding metadata to a NETCDF variable is easy:

   latvar.grads_dim='y'
   latvar.grads_mapping='linear'
   latvar.grads_size=str(len(lat[lmin:lmax]))
   latvar.units='degrees_north'
   latvar.long_name='latitude'
   latvar.minimum=str(lat[lmin])
   latvar.maximum=str(lat[lmax])
   latvar.resolution='0.25'
   latvar[:]=lat[lmin:lmax]

Again, remember with the sub setting, the lmin and lmax indices are required to match the output file.
Now if I print out MY latitude variable within the output NETCDF file it looks similar to the GDS file with only minor differences:

<type 'netCDF4.Variable'>
float64 lat(lat)
    grads_dim: y
    grads_mapping: linear
    grads_size: 200
    units: degrees_north
    long_name: latitude
    minimum: 10.0
    maximum: 60.0
    resolution: 0.25
unlimited dimensions:
current shape = (200,)

You will need to do the same thing for the other dimension variables (lon,lev,time) with a few minor differences.  I don't explicitly show the steps here to save space and redundancy.  You can find the steps written down in the accompanying .py file.

Once you have your dimension variables, you need to add the temperature variables.  In the example, I have you add both the 2 meter temperature and the temperature as a function of pressure level.  This is to show you how to add variables with differing dimensions.

As I stated earlier, there are fewer restrictions for these variables than the dimension ones.  But it is still a good practice to add metadata.  When printing out the GDS 'tmp2m' var, you get this:

<type 'netCDF4.Variable'>
float32 tmp2m(time, lat, lon)
    _FillValue: 9.999e+20
    missing_value: 9.999e+20
    long_name: ** 2 m above ground temperature [k]
unlimited dimensions:
current shape = (81, 721, 1440)

So for our variable, we will try to match this as close as possible:
  gtmpvar=newdata.createVariable('GTMP','f4',('time','lat','lon'),fill_value=9.999E20)
  gtmpvar.long_name='2 Meter Temperature'
  gtmpvar.units='K'
  gtmpvar[:]=var1[:8,lmin:lmax,lomin:lomax] ## Add the actual data from the GDS file

And that's it for the 2 meter temperature.  Now, you can do the exact same thing for the vertically varying temperature with only minor changes:
   tmpvar=newdata.createVariable('TEMP','f4',('time','lev','lat','lon'),fill_value=9.999E20)
   tmpvar.long_name='Pressure Interpolated Temperature'
   tmpvar.units='K'
   tmpvar[:]=var2[:8,:,lmin:lmax,lomin:lomax]

The big difference here, is that in the dimensions section of the createVariable function includes a dimension for 'lev'.

When you are done, close out your file and you can then open it in GrADS using the 'sdfopen' command.  Then if you display the temperature data, you should get something that looks like this:

How to read and write a NON self describing NETCDF file (xdfopen)

If you want to read a file that is not self-describing into GrADS, you need to include a control file.
Essentially, you need to treat the NETCDF file as if it were a binary data file.

Now for this tutorial, we will run through the exact same process, highlighting the differences in how we write metadata to the NETCDF file.

Starting with the dimension variables, we will no longer include the metadata that we did for the self describing file:

    newdata=netCDF4.Dataset(xdffile,'w')
    newdata.createDimension('lat', len(lat[lmin:lmax]))
    latvar=newdata.createVariable('lat','f8',('lat'))
    latvar[:]=lat[lmin:lmax]

Similarly, I'll skip the rest of the dimension variables, but as you see, all I did was define the variable and add the data to it.  For simplicity, and since it doesn't matter, I'll keep the temperature variables as they were in the self describing example.

Now the NETCDF file is created.  If however you try to open this file using the 'sdfopen' command, you will get an error.  So we need to create a control file.

The control file can be easily written using python (as in this example) or whatever program you are using.  I'm going to copy the code to write the control file below, and explain it after.

    ### Now write the control file ###
    xdfctl='C:/OpenGrADS/xdffile.ctl'
    xdflocal='xdffile.nc'
    xdf = open(xdfctl, 'w')
    xdf.write('DSET ^' + xdflocal +'\n')
    xdf.write('TITLE XDF example File 00z30dec2015 \n')
    xdf.write('UNDEF 99999.0\n')
    xdf.write('XDEF lon '+str(len(lon[lomin:lomax]))+' LINEAR ' + str(lon[lomin]) +' 0.25\n')
    xdf.write('YDEF lat '+str(len(lat[lmin:lmax]))+' LINEAR ' + str(lat[lmin]) +' 0.25\n')
    xdf.write('ZDEF lev '+str(len(lev[:]))+' LEVELS ')
    for i in lev[:]:
        xdf.write(str(int(i))+' ')
    xdf.write('\n')
    xdf.write('TDEF time 8 LINEAR 00z30dec2015 3hr\n')
    xdf.write('VARS 2\n')
    xdf.write('GTMP=>GTMP 0 99 2-meter Temperature [k]\n')
    xdf.write('TEMP=>GTMP '+str(len(lev[:]))+' 99 Temperature [K]\n')
    xdf.write('ENDVARS')
    xdf.close()

The most important take away to learn about writing a control for a NETCDF file is that you match the variable/dimension names in the control file to the variable/dimension names in the NETCDF file (case sensitive).

For example: 'GTMP=>GTMP' will only work if there is a variable in the NETCDF file called 'GTMP'.  Similarly, following the XDEF definition, you need to point to the longitude variable (lon) in the NETCDF file.  This is different than writing out control files for binary datasets which don't require anything after XDEF.

Other than that, you see some other things, the levels are defined in a special way since the vertical coordinate is not incremented linearly.  Also the lat/lon sizes are not hard coded, to allow you to make bigger or smaller spatial subsets without a lot of extra work.

Open the control file in GrADS using the 'xdfopen' command:

'xdfopen xdffile.ctl'

A few notes on using the control file:

Note that the coordinate data is not actually read
from the NETCDF file, it is inferred from the control file.  So it is important that the resolution and the starting points are right, otherwise you will get an incorrect result, even if your longitude data is correct in the NETCDF file.  For example, comparing the correct (0.25) resolution to an incorrect resolution (0.45) gives the two different maps shown on the right.  So get the resolution right!

One final note: This kind of control file only seems to work with nice rectangular datasets.  That is lat/lon coordinates that vary independently from one another.  If your lat/lon grid does not vary independently (i.e., longitude is a function of y and visa-versa), then using this method will warp tour map projection.  In these situations it is recommended that you write your data out to a binary file with the PDEF option.  I have yet to get the PDEF option to work with a NETCDF file.






So that is it for the tutorial, hopefully you have a better understanding regarding the difference between self describing and non self describing NETCDF files, and how to write and open them in GrADS.  No GrADS scripts to download here, but a full python script (version 2.7) is included that will reproduce the files and plots generated in this tutorial.  You will need the netCDF4 module installed for this script to work however.  Numpy and sys are standard I believe.  Also you may need to update the date.

Download Python Script here

Read More