Page History

...

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. 

   Info
   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: 

   Code Block
   language py linenumbers true
   import matplotlib.pyplot as plt import netCDF4   ncfile = './Water level (op).nc'

4. To display the contents of the file, use the following build in function. This will print a list of all variables and sizes to the command window. Next, we build a dataset object and print all the variables in the netcdf file

   Code Block
   linenumbers true
   ncdisp(ncfile)ncfid = netCDF4.Dataset(path) 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'! Check this using ncdisp. To Let's change this to MATLAB compatible time, use the following codepython datetime objects

      Code Block
      linenumbers true
      %# Read data from nc file time = ncread(ncfile, np.array(ncfid.variables['time'])   %# convertConvert to MATLABdatetime datesobjects time = time./(1000*3600*24[datetime(year=1970, month=1, day=1) + datenum('01-01-1970', 'dd-mm-yyyy')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:

      Code Block
      linenumbers true
      %# Read data from nc file waterlevels = ncread(ncfile, 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 cell array immediately: 

      Code Block
      linenumbers true
      % Read data from nc file names = cellstr(ncread(path, 'feature_name')')

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

   Code Block
   linenumbers true
   % Change the index to cycle through the stations station_index = 1 % Plot the data figure(1), clf(1) plot(time, waterlevels(station_index, :), 'color', [.3 .3 .9], 'linewidth', 2) title(sprintf('Station: %s', names{station_index})) datetick('x', 'dd mmm yyyy' ,'keepticks', 'keeplimits') ylabel('Water level [m]')   % Some extra code to make the figure less default-looking set(gca, ... 'Box' , 'on' , ... 'TickDir' , 'out' , ... 'TickLength' , [.02 .02] , ... 'XMinorTick' , 'on' , ... 'YMinorTick' , 'on' , ... 'YGrid' , 'on' , ... 'XColor' , [.3 .3 .3], ... 'YColor' , [.3 .3 .3], ... 'LineWidth' , 1 );

...