This tutorial requires Python 2.7 with the NetCDF4 and Matplotlib (1.5.1 or higher) packages. Python is a free and open source language. If you do not have Python installed, we recommend installing the Anaconda distribution.

This tutorial assumed basic knowledge of Python.

Note: this tutorial does not use the build-in IronPython distribution of DeltaShell.

Step-by-step guide

SOBEK outputs NetCDF files (*.nc). This tutorial shows you how to access these files in Python. We will investigate the contents of the netCDF file and plot a timeseries. 

 

  1. You can find SOBEK output in the <projectname>.dsproj_data folder. Here you will find several files with the *.nc extension. These are NetCDF files. They correspond to the output selected in DeltaShell. Note: if this folder contains files with the extension *.nc.changes you have not saved your project after running the model. Save the project in DeltaShell before continuing. 

    NetCDF (Network Common Data Form) is an open standard for storing scientific data. Delft3D Flexible Mesh (and SOBEK) use the CF-1.0 convention.

  2. For this tutorial we will use the observation point output for water level. If you do not have a SOBEK model readily available, download the following output file: Water level (op).nc

  3. Open your Python editor of choice. First we import the necessary modules and define a variable pointing to the *.nc file: 

    import matplotlib.pyplot as plt
import netCDF4
import numpy as np
from datetime import datetime, timedelta

ncfile = './Water level (op).nc'

  4. Next, we build a dataset object and print all the variables in the netcdf file

    ncfid = netCDF4.Dataset(ncfile)
for variable in ncfid.variables:
	print variable

  5. In the next steps, we will extract the times, values and observation point names from the NetCDF file

    1. The variable for the timestamps is called 'time', but note that it is in 'milliseconds since 1970-01-01'! Let's change this to python datetime objects

      # Read data from nc file
time = np.array(ncfid.variables['time'])
 
# Convert to datetime objects
time = [datetime(year=1970, month=1, day=1) + timedelta(seconds=t/1000.) for t in time]

    2. The values (which are water levels, in this case) are easily extracted. We convert it to a numpy array to easily transpose:

      # Read data from nc file
waterlevels = np.array(ncfid.variables['value']).T

    3. Finally, the names of the observation stations are stored under the attribute 'feature_name'. Reading this data will return a transposed character array. It is convenient to transform this to a proper list immediatly: 

      # Read data from nc file
names = [''.join(i).strip() for i in ncfid.variables['feature_name']]

  6. With the data read from file, we are ready to plot the result:

    # Note in Python the first index is 0 (in MATLAB it is 1)
station_index = 0
fig, ax = plt.subplots()
ax.plot(time, waterlevels[station_index])
plt.title(names[station_index])
 
# Tweak the appearance of the plot
ax.get_yaxis().set_tick_params(direction='out')
ax.get_xaxis().set_tick_params(direction='out')
plt.grid(b=True, which='both', color='0.65',linestyle='-')
for spine in ['bottom', 'top', 'right', 'left']:
	ax.spines[spine].set_color((0.3, 0.3, 0.3))

plt.show()




    That's it! 

Disclaimer