How to read NetCDF file and write to CSV using Python

I would load the data into Pandas, which facilitates the analysis and plotting of time series data, as well as writing to CSV.

So here's a real working example which pulls a time series of wave heights from a specified lon,lat location out of a global forecast model dataset.

Note: here we access an OPeNDAP dataset so we can just extract the data we need from a remote server without downloading files. But netCDF4 works exactly the same for a remove OPeNDAP dataset or a local NetCDF file, which is a very useful feature!

import netCDF4
import pandas as pd
import matplotlib.pyplot as plt

# NetCDF4-Python can read a remote OPeNDAP dataset or a local NetCDF file:
url='http://thredds.ucar.edu/thredds/dodsC/grib/NCEP/WW3/Global/Best'
nc = netCDF4.Dataset(url)
nc.variables.keys()

lat = nc.variables['lat'][:]
lon = nc.variables['lon'][:]
time_var = nc.variables['time']
dtime = netCDF4.num2date(time_var[:],time_var.units)

# determine what longitude convention is being used [-180,180], [0,360]
print lon.min(),lon.max()

# specify some location to extract time series
lati = 41.4; loni = -67.8 +360.0  # Georges Bank

# find closest index to specified value
def near(array,value):
    idx=(abs(array-value)).argmin()
    return idx

# Find nearest point to desired location (could also interpolate, but more work)
ix = near(lon, loni)
iy = near(lat, lati)

# Extract desired times.      
# 1. Select -+some days around the current time:
start = dt.datetime.utcnow()- dt.timedelta(days=3)
stop = dt.datetime.utcnow()+ dt.timedelta(days=3)
#       OR
# 2. Specify the exact time period you want:
#start = dt.datetime(2013,6,2,0,0,0)
#stop = dt.datetime(2013,6,3,0,0,0)

istart = netCDF4.date2index(start,time_var,select='nearest')
istop = netCDF4.date2index(stop,time_var,select='nearest')
print istart,istop

# Get all time records of variable [vname] at indices [iy,ix]
vname = 'Significant_height_of_wind_waves_surface'
#vname = 'surf_el'
var = nc.variables[vname]
hs = var[istart:istop,iy,ix]
tim = dtime[istart:istop]

# Create Pandas time series object
ts = pd.Series(hs,index=tim,name=vname)

# Use Pandas time series plot method
ts.plot(figsize(12,4),
   title='Location: Lon=%.2f, Lat=%.2f' % ( lon[ix], lat[iy]),legend=True)
plt.ylabel(var.units);

#write to a CSV file
ts.to_csv('time_series_from_netcdf.csv')

which both creates this plot to verify that you've got the data you wanted:

and also writes the desired CSV file time_series_from_netcdf.csv to disk.

You can also view, download and/or run this example on Wakari.

Rich Signell's answer was incredibly helpful! Just as a note, it's important to also import datetime, and when extracting times, it's necessary to use the following code:

import datetime
import netCDF4
import pandas as pd
import matplotlib.pyplot as plt

...

# 2. Specify the exact time period you want:
start = datetime.datetime(2005,1,1,0,0,0)
stop = datetime.datetime(2010,12,31,0,0,0)

I then looped over all the regions that I needed for my dataset.

Not sure what you're still having trouble with, this looks good. I do see:

# convert date, how to store date only strip away time?
 print "Converting Dates"
 units = nc.variables['time'].units
 dates = num2date (times[:], units=units, calendar='365_day')

you now have the dates as python datetime objects

 #print [dates.strftime('%Y%m%d%H') for date in dates]

and this is what you need if you want them as strings -- but if you only want the day, remove the %H:

date_strings = [dates.strftime('%Y%m%d') for date in dates]

if you want the year, month day as numbers, datetime objects have attributes for that:

dt.year, dt.month, dt.day

As for your sfc variable -- is a 3-d array, so to get a particular value, you can do:

sfc[time_index, lat_index, lon_index]

being 3-d there are more than one way to write it to a csv file, but I'm guessing you might want something like:

for time_index, time in enumerate(time): # pull out the data for that time data = sfc[time_index, :, :] # write the date to the file (maybe) # .... Now loop through the "rows" for row in data: outputwriter.writerow( [str(val) for val in row] )

Or something like that....

The problem with the attribute error is because content needs to be a list, and you initialize it with lat, which is just a number. You need to stuff that into a list.

Regarding the 3D variables, lats = nc.variables['latitude'][:] is sufficient to read all the data.

Update: Iterate over lon/lat together

Here's your code with the mod for the list and iteration:

# the csv file is closed when you leave the block
with open('output.csv', 'wb') as csvFile:
    outputwriter = csv.writer(csvFile, delimiter=',')
    for time_index, time in enumerate(times): # pull the dates out for the header
        t = num2date(time, units = units, calendar='365_day')
        header.append(t)
    outputwriter.writerow(header)

    for latlon_index, (lat,lon) in enumerate(zip(lats, lons)):
        content = [lat, lon] # Put lat and lon into list
        print latlon_index
        for time_index, time in enumerate(times): # for a date
            # pull out the data 
            data = sfc[time_index,lat_index,lon_index]
            content.append(data)
            outputwriter.writerow(content)``

I haven't actually tried to run this, so there may be other problems lurking.