|
Here is an example of how to read and write data with Unidata NetCDF (Network Common Data Form) files using the NetCDF4 Python module. In this example, I use a NetCDF file of 2012 air temperature on the 0.995 sigma level ('./air.sig995.2012.nc') from the NCEP/NCAR Reanalysis I (Kalnay et al. 1996) [NCEP/NCAR Reanalysis data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their Web site at http://www.esrl./psd/]. With the data, we make a contour fill map in a Mollweide equal area projection using the Matplotlib toolkit Basemap. After this, we will write the air temperature profile for Darwin, Australia (12.45°S, 130.83°E) into a NetCDF file entitled './darwin_2012.nc.' The temperature profile has a fixed location and is only dependent on the time dimension. We will create a simple line plot to visualize this data. [Please note: the NetCDF file generated in this example is not CF (Climate and Forecast) metadata convention compliant. Visit for more details.] Lastly, we will compute the global air temperature departure from its value at Darwin, Australia for all of 2012. We will create a corresponding NetCDF file entitled './air.departure.sig995.2012.nc' with the air temperature departure as a function of time, latitude, and longitude. In addition, we will create a contour fill plot of the departure. NCEP/NCAR Reanalysis I project products and documentation can be found at http://www.esrl./psd/data/gridded/data.ncep.reanalysis.html. For more examples on using NetCDF and Python, visit the Unidata NetCDF example programs page at http://www.unidata./software/netcdf/examples/programs/. Please feel free to contact me with any feedback, questions, comments, or concerns. My contact information can be found on my about page. Please note that Python - NetCDF reading and writing example with plotting by Chris Slocum is licensed under a Creative Commons Attribution-NonCommercial-NoDerivs 3.0 Unported License. Example output ¶FiguresAfter the data are read using Python, the air temperature is plotted using a Mollweide projection.  2012 temperature profile for Darwin, Australia.  The global air temperature departure from its value at Darwin, Australia.  Standard OutputNetCDF dimension and variable attributes for the './air.sig995.2012.nc' which is the 2012 air temperature output for the 0.9950 sigma level of the NCEP/NCAR Reanalysis. You can use the function ncdump() from the source code below with any NetCDF file to output similar file attribute information.NetCDF Global Attributes:
Conventions: u'COARDS'
title: u'mean daily NMC reanalysis (2012)'
history: u'created 2011/12 by Hoop (netCDF2.3)'
description: u'Data is from NMC initialized reanalysis\n(4x/day). These are the 0.9950 sigma level values.'
platform: u'Model'
references: u'http://www.esrl./psd/data/gridded/data.ncep.reanalysis.html'
NetCDF dimension information:
Name: lat
size: 73
type: dtype('float32')
units: u'degrees_north'
actual_range: array([ 90., -90.], dtype=float32)
long_name: u'Latitude'
standard_name: u'latitude'
axis: u'Y'
Name: lon
size: 144
type: dtype('float32')
units: u'degrees_east'
long_name: u'Longitude'
actual_range: array([ 0. , 357.5], dtype=float32)
standard_name: u'longitude'
axis: u'X'
Name: time
size: 366
type: dtype('float64')
units: u'hours since 1-1-1 00:00:0.0'
long_name: u'Time'
actual_range: array([ 17628096., 17636856.])
delta_t: u'0000-00-01 00:00:00'
standard_name: u'time'
axis: u'T'
avg_period: u'0000-00-01 00:00:00'
NetCDF variable information:
Name: air
dimensions: (u'time', u'lat', u'lon')
size: 3847392
type: dtype('int16')
long_name: u'mean Daily Air temperature at sigma level 995'
unpacked_valid_range: array([ 185.16000366, 331.16000366], dtype=float32)
actual_range: array([ 193.80001831, 317.52001953], dtype=float32)
units: u'degK'
add_offset: 512.81
scale_factor: 0.0099999998
missing_value: 32766
precision: 2
least_significant_digit: 1
GRIB_id: 11
GRIB_name: u'TMP'
var_desc: u'Air temperature'
dataset: u'NCEP Reanalysis Daily Averages'
level_desc: u'Surface'
statistic: u'Mean\nM'
parent_stat: u'Individual Obs\nI'
valid_range: array([-32765, -18165], dtype=int16)
Source Code¶'''
NAME
NetCDF with Python
PURPOSE
To demonstrate how to read and write data with NetCDF files using
a NetCDF file from the NCEP/NCAR Reanalysis.
Plotting using Matplotlib and Basemap is also shown.
PROGRAMMER(S)
Chris Slocum
REVISION HISTORY
20140320 -- Initial version created and posted online
20140722 -- Added basic error handling to ncdump
Thanks to K.-Michael Aye for highlighting the issue
REFERENCES
netcdf4-python -- http://code.google.com/p/netcdf4-python/
NCEP/NCAR Reanalysis -- Kalnay et al. 1996
http://dx./10.1175/1520-0477(1996)077<0437:TNYRP>2.0.CO;2
'''
import datetime as dt # Python standard library datetime module
import numpy as np
from netCDF4 import Dataset # http://code.google.com/p/netcdf4-python/
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, addcyclic, shiftgrid
def ncdump(nc_fid, verb=True):
'''
ncdump outputs dimensions, variables and their attribute information.
The information is similar to that of NCAR's ncdump utility.
ncdump requires a valid instance of Dataset.
Parameters
----------
nc_fid : netCDF4.Dataset
A netCDF4 dateset object
verb : Boolean
whether or not nc_attrs, nc_dims, and nc_vars are printed
Returns
-------
nc_attrs : list
A Python list of the NetCDF file global attributes
nc_dims : list
A Python list of the NetCDF file dimensions
nc_vars : list
A Python list of the NetCDF file variables
'''
def print_ncattr(key):
"""
Prints the NetCDF file attributes for a given key
Parameters
----------
key : unicode
a valid netCDF4.Dataset.variables key
"""
try:
print "\t\ttype:", repr(nc_fid.variables[key].dtype)
for ncattr in nc_fid.variables[key].ncattrs():
print '\t\t%s:' % ncattr, repr(nc_fid.variables[key].getncattr(ncattr))
except KeyError:
print "\t\tWARNING: %s does not contain variable attributes" % key
# NetCDF global attributes
nc_attrs = nc_fid.ncattrs()
if verb:
print "NetCDF Global Attributes:"
for nc_attr in nc_attrs:
print '\t%s:' % nc_attr, repr(nc_fid.getncattr(nc_attr))
nc_dims = [dim for dim in nc_fid.dimensions] # list of nc dimensions
# Dimension shape information.
if verb:
print "NetCDF dimension information:"
for dim in nc_dims:
print "\tName:", dim
print "\t\tsize:", len(nc_fid.dimensions[dim])
print_ncattr(dim)
# Variable information.
nc_vars = [var for var in nc_fid.variables] # list of nc variables
if verb:
print "NetCDF variable information:"
for var in nc_vars:
if var not in nc_dims:
print '\tName:', var
print "\t\tdimensions:", nc_fid.variables[var].dimensions
print "\t\tsize:", nc_fid.variables[var].size
print_ncattr(var)
return nc_attrs, nc_dims, nc_vars
nc_f = './air.sig995.2012.nc' # Your filename
nc_fid = Dataset(nc_f, 'r') # Dataset is the class behavior to open the file
# and create an instance of the ncCDF4 class
nc_attrs, nc_dims, nc_vars = ncdump(nc_fid)
# Extract data from NetCDF file
lats = nc_fid.variables['lat'][:] # extract/copy the data
lons = nc_fid.variables['lon'][:]
time = nc_fid.variables['time'][:]
air = nc_fid.variables['air'][:] # shape is time, lat, lon as shown above
time_idx = 237 # some random day in 2012
# Python and the renalaysis are slightly off in time so this fixes that problem
offset = dt.timedelta(hours=48)
# List of all times in the file as datetime objects
dt_time = [dt.date(1, 1, 1) + dt.timedelta(hours=t) - offset for t in time]
cur_time = dt_time[time_idx]
# Plot of global temperature on our random day
fig = plt.figure()
fig.subplots_adjust(left=0., right=1., bottom=0., top=0.9)
# Setup the map. See http:///basemap/users/mapsetup.html
# for other projections.
m = Basemap(projection='moll', llcrnrlat=-90, urcrnrlat=90, llcrnrlon=0, urcrnrlon=360, resolution='c', lon_0=0)
m.drawcoastlines()
m.drawmapboundary()
# Make the plot continuous
air_cyclic, lons_cyclic = addcyclic(air[time_idx, :, :], lons)
# Shift the grid so lons go from -180 to 180 instead of 0 to 360.
air_cyclic, lons_cyclic = shiftgrid(180., air_cyclic, lons_cyclic, start=False)
# Create 2D lat/lon arrays for Basemap
lon2d, lat2d = np.meshgrid(lons_cyclic, lats)
# Transforms lat/lon into plotting coordinates for projection
x, y = m(lon2d, lat2d)
# Plot of air temperature with 11 contour intervals
cs = m.contourf(x, y, air_cyclic, 11, cmap=plt.cm.Spectral_r)
cbar = plt.colorbar(cs, orientation='horizontal', shrink=0.5)
cbar.set_label("%s (%s)" % (nc_fid.variables['air'].var_desc, nc_fid.variables['air'].units))
plt.title("%s on %s" % (nc_fid.variables['air'].var_desc, cur_time))
# Writing NetCDF files
# For this example, we will create two NetCDF4 files. One with the global air
# temperature departure from its value at Darwin, Australia. The other with
# the temperature profile for the entire year at Darwin.
darwin = {'name': 'Darwin, Australia', 'lat': -12.45, 'lon': 130.83}
# Find the nearest latitude and longitude for Darwin
lat_idx = np.abs(lats - darwin['lat']).argmin()
lon_idx = np.abs(lons - darwin['lon']).argmin()
# Simple example: temperature profile for the entire year at Darwin.
# Open a new NetCDF file to write the data to. For format, you can choose from
# 'NETCDF3_CLASSIC', 'NETCDF3_64BIT', 'NETCDF4_CLASSIC', and 'NETCDF4'
w_nc_fid = Dataset('darwin_2012.nc', 'w', format='NETCDF4')
w_nc_fid.description = "NCEP/NCAR Reanalysis %s from its value at %s. %s" % (nc_fid.variables['air'].var_desc.lower(), darwin['name'], nc_fid.description)
# Using our previous dimension info, we can create the new time dimension
# Even though we know the size, we are going to set the size to unknown
w_nc_fid.createDimension('time', None)
w_nc_dim = w_nc_fid.createVariable('time', nc_fid.variables['time'].dtype, ('time',))
# You can do this step yourself but someone else did the work for us.
for ncattr in nc_fid.variables['time'].ncattrs():
w_nc_dim.setncattr(ncattr, nc_fid.variables['time'].getncattr(ncattr))
# Assign the dimension data to the new NetCDF file.
w_nc_fid.variables['time'][:] = time
w_nc_var = w_nc_fid.createVariable('air', 'f8', ('time'))
w_nc_var.setncatts({'long_name': u"mean Daily Air temperature", 'units': u"degK", 'level_desc': u'Surface', 'var_desc': u"Air temperature", 'statistic': u'Mean\nM'})
w_nc_fid.variables['air'][:] = air[time_idx, lat_idx, lon_idx]
w_nc_fid.close() # close the new file
# A plot of the temperature profile for Darwin in 2012
fig = plt.figure()
plt.plot(dt_time, air[:, lat_idx, lon_idx], c='r')
plt.plot(dt_time[time_idx], air[time_idx, lat_idx, lon_idx], c='b', marker='o')
plt.text(dt_time[time_idx], air[time_idx, lat_idx, lon_idx], cur_time, ha='right')
fig.autofmt_xdate()
plt.ylabel("%s (%s)" % (nc_fid.variables['air'].var_desc, nc_fid.variables['air'].units))
plt.xlabel("Time")
plt.title("%s from\n%s for %s" % (nc_fid.variables['air'].var_desc, darwin['name'], cur_time.year))
# Complex example: global temperature departure from its value at Darwin
departure = air[:, :, :] - air[:, lat_idx, lon_idx].reshape((time.shape[0], 1, 1))
# Open a new NetCDF file to write the data to. For format, you can choose from
# 'NETCDF3_CLASSIC', 'NETCDF3_64BIT', 'NETCDF4_CLASSIC', and 'NETCDF4'
w_nc_fid = Dataset('air.departure.sig995.2012.nc', 'w', format='NETCDF4')
w_nc_fid.description = "The departure of the NCEP/NCAR Reanalysis " + "%s from its value at %s. %s" % (nc_fid.variables['air'].var_desc.lower(), darwin['name'], nc_fid.description)
# Using our previous dimension information, we can create the new dimensions
data = {}
for dim in nc_dims:
w_nc_fid.createDimension(dim, nc_fid.variables[dim].size)
data[dim] = w_nc_fid.createVariable(dim, nc_fid.variables[dim].dtype, (dim,))
# You can do this step yourself but someone else did the work for us.
for ncattr in nc_fid.variables[dim].ncattrs():
data[dim].setncattr(ncattr, nc_fid.variables[dim].getncattr(ncattr))
# Assign the dimension data to the new NetCDF file.
w_nc_fid.variables['time'][:] = time
w_nc_fid.variables['lat'][:] = lats
w_nc_fid.variables['lon'][:] = lons
# Ok, time to create our departure variable
w_nc_var = w_nc_fid.createVariable('air_dep', 'f8', ('time', 'lat', 'lon'))
w_nc_var.setncatts({'long_name': u"mean Daily Air temperature departure", 'units': u"degK", 'level_desc': u'Surface', 'var_desc': u"Air temperature departure", 'statistic': u'Mean\nM'})
w_nc_fid.variables['air_dep'][:] = departure
w_nc_fid.close() # close the new file
# Rounded maximum absolute value of the departure used for contouring
max_dep = np.round(np.abs(departure[time_idx, :, :]).max()+5., decimals=-1)
# Generate a figure of the departure for a single day
fig = plt.figure()
fig.subplots_adjust(left=0., right=1., bottom=0., top=0.9)
m = Basemap(projection='moll', llcrnrlat=-90, urcrnrlat=90, llcrnrlon=0, urcrnrlon=360, resolution='c', lon_0=0)
m.drawcoastlines()
m.drawmapboundary()
dep_cyclic, lons_cyclic = addcyclic(departure[time_idx, :, :], lons)
dep_cyclic, lons_cyclic = shiftgrid(180., dep_cyclic, lons_cyclic, start=False)
lon2d, lat2d = np.meshgrid(lons_cyclic, lats)
x, y = m(lon2d, lat2d)
levels = np.linspace(-max_dep, max_dep, 11)
cs = m.contourf(x, y, dep_cyclic, levels=levels, cmap=plt.cm.bwr)
x, y = m(darwin['lon'], darwin['lat'])
plt.plot(x, y, c='c', marker='o')
plt.text(x, y, 'Darwin,\nAustralia', color='r', weight='semibold')
cbar = plt.colorbar(cs, orientation='horizontal', shrink=0.5)
cbar.set_label("%s departure (%s)" % (nc_fid.variables['air'].var_desc, nc_fid.variables['air'].units))
plt.title("Departure of Global %s from\n%s for %s" % (nc_fid.variables['air'].var_desc, darwin['name'], cur_time))
plt.show()
# Close original NetCDF file.
nc_fid.close()
Downloads ¶Image files: NetCDF files: Source code: |
|
|
air.sig995.2012.nc
