Script used to generate the netCDF files of the Savi et al experimental data.

cube_Savi_data.py

"""
Script used to generate the netCDF files of the Savi et al experimental data.

Frames extracted from T_NC2_Complete21fps.wmv video file using the command:

     ffmpeg -i T_NC2_Complete21fps.wmv -r 21 T_NC2_frames/%04d.png

after a subdirectory T_NC2_frames was created to hold the output png files.
"""
import os
import numpy as np
import matplotlib.pyplot as plt
import tifffile as tiff
from netCDF4 import Dataset
import time

# two paths to be set to the locations of the DEM and image data
basepath = 'T_NC2/'
vidpath = 'T_NC2_frames/'

# first handle the scan data
flist = os.listdir(basepath)
# shorten list to tif files only
short_list = []
for i in flist:
    if i[-3:] == 'tif':
        short_list.append(i)
# sort
short_list.sort()
# drop duplicate scan (#11)
del short_list[9]
# read data into a 3-D t-x-y array
data = np.zeros((len(short_list), 3400, 1700))
# loop and place data into array
for i in range(len(short_list)):
    data[i, ...] = tiff.imread(basepath + short_list[i]).T

# define vectors of time and space dimensions
# scans taken every 30 minutes
t = np.arange(0, 480+30, 30)
# space dimensions - scan resolution in mm
x = np.arange(data.shape[1])
y = np.arange(data.shape[2])

# create and fill the netcdf file
scan_file = Dataset('T_NC2_scans.nc', 'w', format='NETCDF4')

# description information
scan_file.description = 'Scan data from experiment T_NC2'
scan_file.history = 'Created ' + time.ctime(time.time())
scan_file.source = 'Savi et al 2020; https://doi.org/10.5194/esurf-8-303-2020'

# create time and spatial netCDF dimensions
scan_file.createDimension('time', data.shape[0])
scan_file.createDimension('x', data.shape[1])
scan_file.createDimension('y', data.shape[2])

# create time and space variables
v_time = scan_file.createVariable(
    'time', 'f4', ('time',))
v_time.units = 'minute'
v_x = scan_file.createVariable(
    'x', 'f4', ('x'))
v_x.units = 'millimeter'
v_y = scan_file.createVariable(
    'y', 'f4', ('y'))
v_y.units = 'millimeter'

# fill the variables with the coordinate information
v_time[:] = t
v_x[:] = x
v_y[:] = y

# set up variables for output data grids
v_eta = scan_file.createVariable(
    'eta', 'f4', ('time', 'x', 'y'))
v_eta.units = 'millimeter'
v_eta[:] = data

# set up metadata group and populate variables
scan_file.createGroup('meta')
# store input parameters for this run
v_MC_Qw = scan_file.createVariable(  # a scalar, main channel water discharge
    'meta/MC_Qw', 'f4', ())  # no dims for scalar
v_MC_Qw.units = 'mL/s'
v_MC_Qw[:] = 63
v_MC_Qs = scan_file.createVariable(  # a scalar, main channel sed discharge
    'meta/MC_Qs', 'f4', ())  # no dims for scalar
v_MC_Qs.units = 'mL/s'
v_MC_Qs[:] = 1.3
v_T_Qw = scan_file.createVariable(  # a scalar, tributary water discharge
    'meta/T_Qw', 'f4', ())  # no dims for scalar
v_T_Qw.units = 'mL/s'
v_T_Qw[:] = 41.5
v_T_Qs = scan_file.createVariable(  # a scalar, tributary sediment discharge
    'meta/T_Qs', 'f4', ())  # no dims for scalar
v_T_Qs.units = 'mL/s'
v_T_Qs[:] = 2.2

# close the netcdf file
scan_file.close()


# now handle the overhead photos (extracted from the video)
# load one image to get shape of data
img = plt.imread(vidpath + '0015.png')
# only going to take every 3rd image so get list of subset values
subset = np.arange(1, 1466, 3)
# initialize data arrays for the three color bands (RGB)
red = np.zeros((len(subset), img.shape[1], img.shape[0]))
green = np.zeros_like(red)
blue = np.zeros_like(red)
# load and place the data in the correct arrays
for idx, val in enumerate(subset):
    _img = plt.imread(vidpath + str(val).zfill(4) + '.png')
    red[idx, ...] = _img[:, :, 0].T
    green[idx, ...] = _img[:, :, 1].T
    blue[idx, ...] = _img[:, :, 2].T

# define vectors of time and space dimensions
# images roughly once a minute
t = np.arange(0, len(subset))
# backsolve number of millimeters per pixel based on scan information
dy = data.shape[1] / red.shape[1]
dx = data.shape[2] / red.shape[2]
# take average and round as we assume pixels are probably square
px = np.round((dy + dx) / 2, 2)
# space dimensions
y = np.arange(0, red.shape[1]*px, px)
x = np.arange(0, red.shape[2]*px, px)

# create and fill the netcdf file
img_file = Dataset('T_NC2_images.nc', 'w', format='NETCDF4')

# description information
img_file.description = 'Subset of overhead images from experiment T_NC2'
img_file.history = 'Created ' + time.ctime(time.time())
img_file.source = 'Savi et al 2020; https://doi.org/10.5194/esurf-8-303-2020'

# create time and spatial netCDF dimensions
img_file.createDimension('time', red.shape[0])
img_file.createDimension('x', red.shape[1])
img_file.createDimension('y', red.shape[2])

# create time and space variables
v_time = img_file.createVariable(
    'time', 'f4', ('time',))
v_time.units = 'minute'
v_x = img_file.createVariable(
    'x', 'f4', ('x'))
v_x.units = 'millimeter'
v_y = img_file.createVariable(
    'y', 'f4', ('y'))
v_y.units = 'millimeter'

# fill the variables with the coordinate information
v_time[:] = t
v_x[:] = y
v_y[:] = x

# set up variables for output data grids
v_red = img_file.createVariable(
    'red', 'f4', ('time', 'x', 'y'))
v_red.units = 'float'
v_red[:] = red

v_green = img_file.createVariable(
    'green', 'f4', ('time', 'x', 'y'))
v_green.units = 'float'
v_green[:] = green

v_blue = img_file.createVariable(
    'blue', 'f4', ('time', 'x', 'y'))
v_blue.units = 'float'
v_blue[:] = blue

# set up metadata group and populate variables
img_file.createGroup('meta')
# store input parameters for this run
v_MC_Qw = img_file.createVariable(  # a scalar, main channel water discharge
    'meta/MC_Qw', 'f4', ())  # no dims for scalar
v_MC_Qw.units = 'mL/s'
v_MC_Qw[:] = 63
v_MC_Qs = img_file.createVariable(  # a scalar, main channel sed discharge
    'meta/MC_Qs', 'f4', ())  # no dims for scalar
v_MC_Qs.units = 'mL/s'
v_MC_Qs[:] = 1.3
v_T_Qw = img_file.createVariable(  # a scalar, tributary water discharge
    'meta/T_Qw', 'f4', ())  # no dims for scalar
v_T_Qw.units = 'mL/s'
v_T_Qw[:] = 41.5
v_T_Qs = img_file.createVariable(  # a scalar, tributary sediment discharge
    'meta/T_Qs', 'f4', ())  # no dims for scalar
v_T_Qs.units = 'mL/s'
v_T_Qs[:] = 2.2

# close the netcdf file
img_file.close()