Example processing scripts

A set of example python scripts are included to demonstrate how the various modules within the OCWM tools can be used to analyse model data, visualise model data, validate model data against in situ measurements, and convert ADCP data to the standardised netCDF file format and from netCDF to OCMW format. These python scripts all require access to the OCMW tools, so must be run from an environment where the tools have been installed. The example scripts were used for the analysis of an Open Telemac-Mascaret model of the Sound of Islay (soi), but can be used for any model by altering the data paths, input and output file names, etc. The analysis and validation processes assume that a timeseries of model data have been extracted for a specific spatial location using the otm_extract_by_location.py application, and the ADCP data have been converted to the internal OCMW data format using the adcp2ocmw function in the processADCPData module.

The example processing and analysis scripts provided are summarised in Table 4.

Table 4 : Example processing and analysis scripts.

Example

Purpose

Flow Characterisation

Apply flow characterisation to data extracted at a specific height above bed

Tidal Energy

Calculate and compare tidal energy estimates for a specified TEC based on extracted velocity profilies for a specific location

Tide Charts

Generate tide charts from surface elevation data extracted at a sepcific spatial location

Tidal Analysis

Apply tidal analysis to a timeseries of surface elevation data for a subset of model nodes

Map 2D model fields

Mapping 2D model fields extracted directly from the output Selafin file

Model Validation

Generate validation statitics from a comparison of extracted model data with in situ ADCP data

Convert ADCP to NetCDF

Method for converting ADCP data to standardised netCDF daily files

All of the python scripts take the standard form:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
'Description of script'
"""

# ----------------------------------------------------------------------------
#   IMPORTS
# ----------------------------------------------------------------------------
'Dependent package imports sorted by type'
# Standard Python Dependencies
# Non-Standard Python Dependencies
# Local Module Dependencies
# Other Dependencies

# ----------------------------------------------------------------------------
#   GLOBAL VARIABLES
# ----------------------------------------------------------------------------
'Global variable definitions - if used'

# ----------------------------------------------------------------------------
#   CLASS DEFINITIONS
# ----------------------------------------------------------------------------
'Internal class definitions - if used'

# ----------------------------------------------------------------------------
#   FUNCTION DEFINITIONS
# ----------------------------------------------------------------------------
'Internal function definitions - if used'

# ----------------------------------------------------------------------------
#   MAIN PROCESS
# ----------------------------------------------------------------------------
'Main processing code'

#--------------------------------------------------------------------------
# INTERFACE
#--------------------------------------------------------------------------
'Code to allow command line input via a call to "__main__" - if used'

The shebang line at the head of the file allows the script to be run direclty from the command line. Scripts that generate graphical output include a python input request to allow the figures to be seen when the script is run from the command line.

Flow Characterisation

The flow characterisation process takes a timeseries of data for a specified height above the seabea from a set of velocity profiles extracted for a single spatial location. The flow characterisation can be tied to a particular turbine design by extracting the data using the hub_height_data function, or to a user-defined height above bed by using the hab_data function. When calculated for a turbine, the generated characterisation plot includes information on the turbine cut-in and rated speeds.

The velocity data are preprocessed by sorting into ebb and flood phases, then the peak flow directions determined for each phase, flow asymmetry is defined, the percentage of time in the ebb and flood phases, and if the characterisations is tied to a turbine the percentage of time below cut-in, between cut-in and rated, and above rated speeds for both tide phases. These data are collated into a python dictionary structure along with timestamps and metadata describing the source data.

If a turbine is specified, then the power-weighted-rotor-averaged (PWRA) velocity is calculated from the pre-processed data and a separate python dictionary generated.

The function hub_height_data returns three dictionaries: (1) the extracted timeseries, (2) the pre-processed timeseries, and (3) the PWRA timeseries. The function hab_data only returns two dictionaries: (1) the extracted timeseries, and (2) the pre-processed timeseries. The data dictionaries are automatically saved to the ANALYSIS directory under the model data root directory. If the ANALYSIS directory does not exist, it is created.

The pre-processed data are passed to the graphics functions flow_class_diagram, speed_bin_histogram, and directional_spread_plot to present the characterisation in pictorial form.

The script shown in Listing 4 is for a Sound of Islay model where the characterisation is carried out for two spearate locations and for a Nova 100D tidal turbine.

Listing 4 soi_flowChar.py
  1#!/usr/bin/env python3
  2# -*- coding: utf-8 -*-
  3"""
  4Exmple of flow characterisation for data extracted from OTM model output at a
  5specific geosptatial location and at the hub height of a particular class of 
  6tidal energy converter (TEC).
  7
  8Chris Old
  9IES, School of Engineering, University of Edinburgh
 10Sep 2023
 11
 12"""
 13# ----------------------------------------------------------------------------
 14#   IMPORTS
 15# ----------------------------------------------------------------------------
 16
 17# Standard Python Dependencies
 18# Non-Standard Python Dependencies
 19# Local Module Dependencies
 20from ocmw.core.tecSpecs import getTec
 21from ocmw.core.graphics import flow_class_diagram
 22from ocmw.core.graphics import speed_bin_histograms
 23from ocmw.core.graphics import directional_spread_plot
 24from ocmw.dataproc.processModelData import hub_height_data
 25# Other Dependencies
 26
 27
 28# ----------------------------------------------------------------------------
 29#   GLOBAL VARIABLES
 30# ----------------------------------------------------------------------------
 31
 32
 33# ----------------------------------------------------------------------------
 34#   CLASS DEFINITIONS
 35# ----------------------------------------------------------------------------
 36
 37
 38# ----------------------------------------------------------------------------
 39#   FUNCTION DEFINITIONS
 40# ----------------------------------------------------------------------------
 41
 42
 43# ----------------------------------------------------------------------------
 44#   MAIN PROCESS
 45# ----------------------------------------------------------------------------
 46
 47# Load TEC data
 48tecName = 'Nova100D'
 49tec = getTec(tecName)
 50
 51# General file information
 52modelpath = 'D:/Models/soi_hr_msl_04'
 53datasrc = 'SoI MODEL (High Res TS=4s)'
 54prefix = 'soi_hr_msl_04_model'
 55
 56#------------------------------------------------------
 57# South deployment location
 58#------------------------------------------------------
 59filename = prefix+'_merged_southern.mat'
 60dataset = 'South Deployment'
 61outfile = prefix+'_sth_loc.mat'
 62
 63# Extract hub heigth data
 64ts_mod, data_mod, pwra_sth = hub_height_data(modelpath, filename,
 65                                             datasrc, dataset, tec,
 66                                             outfile)
 67
 68# Generate graphics
 69fh_mod = flow_class_diagram(data_mod,tec=tec,velmax=3.5)
 70fh_mod.show()
 71fh_spdh = speed_bin_histograms(data_mod,tec=tec,binwidth=0.1,maxfrac=0.06,maxspd=3.5)
 72fh_spdh.show()
 73fh_dsrd = directional_spread_plot(data_mod)
 74fh_dsrd.show()
 75
 76#------------------------------------------------------
 77# North deployment location
 78#------------------------------------------------------
 79filename = prefix+'_merged_northern.mat'
 80dataset = 'North Deployment'
 81outfile = prefix+'_nth_loc.mat'
 82
 83# Extract hub heigth data
 84ts_mod, data_mod, pwra_nth = hub_height_data(modelpath, filename,
 85                                             datasrc, dataset, tec,
 86                                             outfile)
 87# Generate graphics
 88fh_mod = flow_class_diagram(data_mod,tec=tec,velmax=3.5)
 89fh_mod.show()
 90fh_spdh = speed_bin_histograms(data_mod,tec=tec,binwidth=0.1,maxfrac=0.06,maxspd=3.5)
 91fh_spdh.show()
 92fh_dsrd = directional_spread_plot(data_mod)
 93fh_dsrd.show()
 94
 95print('Process complete.\n')
 96input('Press ENTER to exit...')
 97
 98#--------------------------------------------------------------------------
 99# INTERFACE
100#--------------------------------------------------------------------------
101

Tidal Energy

The estimation of energy yield is a key component of site assessment prior to site development. The yield will depend on the type of tidal turbine deployed and the spatial location of the turbine. Turbine energy yield is calculated from the PWRA velocity for the given turbine, the rotor area, the seawater density, and the turbines power coefficient. The code currently assumes a constant power coefficient, i.e. the coefficient is independent of flow speed or turbine controls.

The script shown in Listing 5 is for a Sound of Islay model where the power generated by a Nova 100D turbine is estimated for two separate locations and the predicted yields compared. For the purpose of this example, it has been assumed that the turbine is always oriented along the direction of peak flow, i.e. the turbine yaws. The power has also been calculated using the depth-averaged velocities to show the discrepancies that would arise from using a 2D depth-averaged model for site assessment.

There are two internal functions (get_soi_hubht_pwra_data, get_soi_depavg_pwr_data) used to determine where the PWRA data have already been extracted. If the PWRA data exist then the file is loaded, otherwise the data are generated prior to the analysis.

Listing 5 soi_Nova100D_tidalEnergy.py
  1#!/usr/bin/env python3
  2# -*- coding: utf-8 -*-
  3"""
  4Example of comparing predicted Turbine power production from data extracted from 
  5OTM model output at a pair of specfic geospatial locations.
  6
  7Chris Old
  8IES, School of Engineering, University of Edinburgh
  9Oct 2023
 10
 11"""
 12# ----------------------------------------------------------------------------
 13#   IMPORTS
 14# ----------------------------------------------------------------------------
 15
 16# Standard Python Dependencies
 17import os
 18# Non-Standard Python Dependencies
 19import numpy as np
 20import matplotlib.pyplot as plt
 21# Local Module Dependencies
 22from ocmw.core.fileTools import getListOfFiles
 23from ocmw.core.timeFuncs import dateStrFromDateNum, dateNumFromDateStr
 24from ocmw.core.flowChar import streamwise_flow, tec_power
 25from ocmw.core.tecSpecs import getTec
 26from ocmw.dataproc.processModelData import hub_height_data, dav_data
 27from ocmw.dataproc.ocmw_extract import load_ocmw_file
 28# Other Dependencies
 29
 30
 31# ----------------------------------------------------------------------------
 32#   GLOBAL VARIABLES
 33# ----------------------------------------------------------------------------
 34
 35
 36# ----------------------------------------------------------------------------
 37#   CLASS DEFINITIONS
 38# ----------------------------------------------------------------------------
 39
 40
 41# ----------------------------------------------------------------------------
 42#   FUNCTION DEFINITIONS
 43# ----------------------------------------------------------------------------
 44def get_soi_hubht_pwra_data(modelpath,filename,datasrc,dataset,tec,outfile):
 45    analpath = modelpath+'/ANALYSIS'
 46    if not os.path.isdir(analpath):
 47        os.mkdir(analpath)
 48    ofprts = outfile.split('.')
 49
 50    hab = tec.bottomMounted
 51    zhub = tec.hubHeight
 52    
 53    if hab:
 54        rstr = '_hab_'
 55    else:
 56        rstr = '_dbs_'
 57    
 58    pwrafile = ofprts[0]+rstr+str(np.int32(np.ceil(zhub))).zfill(2)+'m_PWRA.'+ofprts[1]
 59    
 60    if os.path.isdir(analpath):
 61        pwra_data = load_ocmw_file(analpath,pwrafile)
 62    else:
 63        modelt_ts, model_data, pwra_data = hub_height_data(modelpath,filename,datasrc,dataset,tec,outfile)
 64    return pwra_data
 65
 66
 67def get_soi_depavg_pwr_data(modelpath,filename,datasrc,dataset,tec,outfile):
 68    analpath = modelpath+'/ANALYSIS'
 69    if not os.path.isdir(analpath):
 70        os.mkdir(analpath)
 71    ofprts = outfile.split('.')
 72
 73    rdia = tec.rotorDiam
 74    Cp = tec.Cp
 75    
 76    files = getListOfFiles(analpath,ofprts[0]+'_depth_avg_TS')
 77    if len(files) == 0:
 78        model_ts, model_data = dav_data(modelpath, filename, datasrc, dataset, outfile)
 79    else:
 80        datafile = ofprts[0]+'_depth_avg.'+ofprts[1]
 81        model_data = load_ocmw_file(analpath,datafile)
 82    vel_strm = streamwise_flow(model_data)
 83    pwr_data = tec_power(vel_strm,rdia,Cp=Cp)
 84    return model_data, pwr_data
 85
 86
 87# ----------------------------------------------------------------------------
 88#   MAIN PROCESS
 89# ----------------------------------------------------------------------------
 90
 91# Load TEC data
 92tec = getTec('Nova100D')
 93
 94# General file information
 95modelpath = 'D:/Models/soi_hr_msl_04'
 96datasrc = 'SoI MODEL (High Res TS=4s)'
 97outprefix = 'soi_hr_msl_04_model'
 98
 99# South deployment location
100filename = outprefix+'_merged_sth_loc.mat'
101dataset = 'South Deployment'
102outfile = outprefix+'_sth_loc.mat'
103
104pwra_tecSth = get_soi_hubht_pwra_data(modelpath, filename, datasrc, dataset, tec, outfile)
105
106# North deployment location
107filename = outprefix+'_merged_nth_loc.mat'
108dataset = 'North Deployment'
109outfile = outprefix+'_nth_loc.mat'
110
111pwra_tecNth = get_soi_hubht_pwra_data(modelpath, filename, datasrc, dataset, tec, outfile)
112
113
114# Compare Power based on PWRA
115
116t = pwra_tecSth['times']+dateNumFromDateStr(pwra_tecSth['epoch'].replace('/','-'))
117t0 = np.floor(t[0])
118epstr = dateStrFromDateNum(t0)
119
120fig = plt.figure(figsize=(12,8))
121axs=fig.subplots(2,1)
122dP = pwra_tecNth['tec_pwr']/pwra_tecSth['tec_pwr']
123axs[0].plot(t-t0,pwra_tecSth['tec_pwr']/1.0e3,label='South')
124axs[0].plot(t-t0,pwra_tecNth['tec_pwr']/1.0e3,label='North')
125axs[0].set_xlim([0.0,27.0])
126axs[0].set_ylim([0.0,100.0])
127axs[0].grid('on')
128axs[0].legend()
129axs[0].set_ylabel('TEC Power  [ kW ]',fontsize=8)
130axs[0].set_title('Comparison of Estimated Power Production by South and North Turbines')
131axs[0].set_position([0.075,0.56,0.88,0.37])
132
133dP = pwra_tecNth['tec_pwr']-pwra_tecSth['tec_pwr']
134axs[1].plot(t-t0,dP/1.0e3)
135axs[1].set_xlim([0.0,27.0])
136axs[1].set_ylim([-50.0,0.0])
137axs[1].grid('on')
138axs[1].set_xlabel('Timestamp  [days from '+epstr+']',fontsize=8)
139axs[1].set_ylabel('North-South Power Differece  [ kW ]',fontsize=8)
140axs[1].set_title('Difference in Estimated Power Production by South and North Turbines')
141axs[1].set_position([0.075,0.09,0.88,0.37])
142
143fig.show()
144
145plt.figure(figsize=(9,9))
146plt.plot(pwra_tecSth['tec_pwr']/1.0e3,pwra_tecNth['tec_pwr']/1.0e3,'.')
147ax = plt.gca()
148ax.set_aspect('equal')
149plt.grid('on')
150plt.xlabel('South Power  [ kW ]')
151plt.ylabel('North Power  [ kW ]')
152
153plt.show()
154
155# Energy yield difference
156energy_tec1 = np.nansum(pwra_tecSth['tec_pwr'])*300.0
157energy_tec2 = np.nansum(pwra_tecNth['tec_pwr'])*300.0
158
159yield_difference = 100.*energy_tec2/energy_tec1
160
161print('Energy Yield Difference: North yield is '+"{0:+.1f}".format(yield_difference)+'% of South')
162
163# Exclude speeds below TEC cut-in
164tec1indx = np.where(pwra_tecSth['magn'] < tec.cutInSpeed)[0]
165tec1_pwr = pwra_tecSth['tec_pwr'].copy()
166tec1_pwr[tec1indx] = np.nan
167energy_tec1_limited = np.nansum(tec1_pwr)*300.0
168
169tec2indx = np.where(pwra_tecNth['magn'] < tec.cutInSpeed)[0]
170tec2_pwr = pwra_tecNth['tec_pwr'].copy()
171tec2_pwr[tec2indx] = np.nan
172energy_tec2_limited = np.nansum(tec2_pwr)*300.0
173
174yield_difference_limited = 100.*energy_tec2_limited/energy_tec1_limited
175
176print('Energy Yield Difference (above cut-in speed): North yield is '+"{0:+.1f}".format(yield_difference_limited)+'% of South')
177
178
179
180#=============================================
181# Depth-Averaged Velocity Anaysis
182#=============================================
183
184# South deployment location
185filename = outprefix+'_merged_southern.mat'
186dataset = 'South Deployment'
187outfile = outprefix+'_sth_loc.mat'
188
189data_tec1, pwr_da_tec1 = get_soi_depavg_pwr_data(modelpath, filename, datasrc, dataset, tec, outfile)
190
191# North deployment location
192filename = outprefix+'_merged_northern.mat'
193dataset = 'North Deployment'
194outfile = outprefix+'_nth_loc.mat'
195
196data_tec2, pwr_da_tec2 = get_soi_depavg_pwr_data(modelpath, filename, datasrc, dataset, tec, outfile)
197
198
199
200# Compare Power based on depth average velocity power
201
202t = data_tec1['time']+dateNumFromDateStr(data_tec1['epoch'].replace('/','-'))
203t0 = np.floor(t[0])
204epstr = dateStrFromDateNum(t0)
205
206fig = plt.figure(figsize=(12,8))
207axs=fig.subplots(2,1)
208dP = pwr_da_tec2/pwr_da_tec1
209axs[0].plot(t-t0,pwr_da_tec1/1.0e3,label='South')
210axs[0].plot(t-t0,pwr_da_tec2/1.0e3,label='North')
211axs[0].set_xlim([0.0,27.0])
212axs[0].set_ylim([0.0,200.0])
213axs[0].grid('on')
214axs[0].legend()
215axs[0].set_ylabel('TEC Power  [ kW ]',fontsize=8)
216axs[0].set_title('Comparison of Estimated Power Production by South and North Turbines (Depth Averaged Velocity)')
217axs[0].set_position([0.075,0.56,0.88,0.37])
218
219dP = pwr_da_tec2-pwr_da_tec1
220axs[1].plot(t-t0,dP/1.0e3)
221axs[1].set_xlim([0.0,27.0])
222axs[1].set_ylim([-150.0,0.0])
223axs[1].grid('on')
224axs[1].set_xlabel('Timestamp  [days from '+epstr+']',fontsize=8)
225axs[1].set_ylabel('North-South Power Differece  [ kW ]',fontsize=8)
226axs[1].set_title('Difference in Estimated Power Production by South and North Turbines (Depth Averaged)')
227axs[1].set_position([0.075,0.09,0.88,0.37])
228
229fig.show()
230
231# Energy yield difference
232energy_da_tec1 = np.nansum(pwr_da_tec1)*300.0
233energy_da_tec2 = np.nansum(pwr_da_tec2)*300.0
234
235yield_difference_da = 100.*energy_da_tec2/energy_da_tec1
236print('Depth-Averaged velocity analysis')
237print('Energy Yield Difference: North yield is '+"{0:+.1f}".format(yield_difference_da)+'% of South')
238
239# Exclude speeds below TEC cut-in
240tec1indx = np.where(pwr_da_tec1 < tec.cutInSpeed)[0]
241tec1_da_pwr = pwr_da_tec1.copy()
242tec1_da_pwr[tec1indx] = np.nan
243energy_da_tec1_limited = np.nansum(tec1_da_pwr)*300.0
244
245tec2indx = np.where(pwr_da_tec2 < tec.cutInSpeed)[0]
246tec2_da_pwr = pwr_da_tec2.copy()
247tec2_da_pwr[tec2indx] = np.nan
248energy_da_tec2_limited = np.nansum(tec2_da_pwr)*300.0
249
250yield_difference_limited_da = 100.*energy_da_tec2_limited/energy_da_tec1_limited
251
252print('Energy Yield Difference (above cut-in speed): North yield is '+"{0:+.1f}".format(yield_difference_limited_da)+'% of South')
253
254fig = plt.figure(figsize=(12,4))
255plt.plot(t-t0,pwr_da_tec1/1.0e3,label='South')
256plt.plot(t-t0,pwr_da_tec2/1.0e3,label='North')
257ax = plt.gca()
258ax.set_xlim([0.0,27.0])
259ax.set_ylim([0.0,200.0])
260ax.tick_params(axis='both', labelsize=8)
261ax.set_yticklabels([])
262ax.grid('on')
263plt.legend(fontsize=8)
264plt.xlabel('Timestamp  [days from '+epstr+']',fontsize=8)
265plt.ylabel('TEC Power  [ kW ]',fontsize=8)
266plt.title('Comparison of Estimated Power Production by South and North Turbines (Depth Averaged Velocity)',fontsize=10)
267ax.set_position([0.05,0.14,0.92,0.77])
268
269fig.show()
270
271print('Process complete.\n')
272input('Press ENTER to exit...')
273
274#--------------------------------------------------------------------------
275# INTERFACE
276#--------------------------------------------------------------------------
277

Tide Charts

Prediction of tidal information is important for operations planning in tehmarine environment. Free-surface hydrodynamic model generate both the 2D or 3D velocity fields and a 2D surface elevation field. From these data tidal information can be generated for specific spatial locations. Tidal charts typically present data on a daily basis, this has been followed in both the model runs and the generation of tidal charts.

The script shown in Listing 6 generated tide plots and summaries of the high and low water times and ranges, and the times of slack water for two separate locations based on the surface elevation data and surface and depth averaged velocity data from a Sound of Islay model. The data are calculated and presented for a single day of model run predictions. There is no dedicated graphics function for plotting tide charts, the script includes the plotting function. The charts for both locations are presented in a single plot. The speeds have been converted to knots and the a pair of upper and lower operting speed bands are included to support operations planning.

Listing 6 soi_tideCharts.py
  1#!/usr/bin/env python3
  2# -*- coding: utf-8 -*-
  3"""
  4Example of generating tide chart information from a time series of data 
  5extracted from OTM model output at a specfic geospatial location.
  6
  7Chris Old
  8IES, School of Engineering, University of Edinburgh
  9Aug 2023
 10
 11"""
 12# ----------------------------------------------------------------------------
 13#   IMPORTS
 14# ----------------------------------------------------------------------------
 15
 16# Standard Python Dependencies
 17# Non-Standard Python Dependencies
 18import matplotlib.pyplot as plt
 19import numpy as np
 20import utm
 21# Local Module Dependencies
 22from ocmw.core.timeFuncs import dateNumFromDateStr, datetimeFromNum
 23from ocmw.core.timeFuncs import dateStrFromDateNum
 24from ocmw.core.tidal import get_slack_water_times, get_tidal_range, print_tide_times
 25from ocmw.dataproc.processModelData import dav_data, srf_data
 26# Other Dependencies
 27
 28
 29# ----------------------------------------------------------------------------
 30#   GLOBAL VARIABLES
 31# ----------------------------------------------------------------------------
 32
 33
 34# ----------------------------------------------------------------------------
 35#   CLASS DEFINITIONS
 36# ----------------------------------------------------------------------------
 37
 38
 39# ----------------------------------------------------------------------------
 40#   FUNCTION DEFINITIONS
 41# ----------------------------------------------------------------------------
 42
 43
 44# ----------------------------------------------------------------------------
 45#   MAIN PROCESS
 46# ----------------------------------------------------------------------------
 47
 48#-------------------------------------------
 49# Load data
 50#-------------------------------------------
 51
 52# General file information
 53modelpath = 'D:/Models/soi_hr_msl_04'
 54datasrc = 'SoI HR MSL Model ( TS=4s )'
 55prefix = 'soi_hr_msl_04_model'
 56
 57#-------------------------------------------
 58# Southern Location
 59#-------------------------------------------
 60filename =prefix+'_merged_southern.mat'
 61dataset = 'Sth_Loc'
 62outfile = prefix+'_south_slack.mat'
 63
 64ts_dav, south_dav = dav_data(modelpath,filename,datasrc,dataset,outfile)
 65south_vmag = south_dav['magn']
 66south_elev = south_dav['surface_elev']
 67
 68ts_srf, south_srf = srf_data(modelpath,filename,datasrc,dataset,outfile)
 69south_vsrf = south_srf['magn']
 70
 71if ts_dav['crs'] in ['EPSG:32630']:
 72    south_loc = np.flip(np.asarray(utm.to_latlon(ts_dav['dataLoc'][0],ts_dav['dataLoc'][1],30,'N')))
 73else:
 74    south_loc = ts_dav['dataLoc']
 75south_msl = np.mean(ts_dav['depth'])
 76
 77dnum = south_dav['times']+dateNumFromDateStr(south_dav['epoch'].replace('/','-'))
 78dt = []
 79for t in dnum:
 80    dt.append(datetimeFromNum(t))
 81
 82#-------------------------------------------
 83# Northern Location
 84#-------------------------------------------
 85filename =prefix+'_merged_northern.mat'
 86dataset = 'Nth_Loc'
 87outfile = prefix+'_north_slack.mat'
 88
 89ts_dav, north_dav = dav_data(modelpath,filename,datasrc,dataset,outfile)
 90north_vmag = north_dav['magn']
 91north_elev = north_dav['surface_elev']
 92
 93ts_srf, north_srf = srf_data(modelpath,filename,datasrc,dataset,outfile)
 94north_vsrf = north_srf['magn']
 95
 96if ts_dav['crs'] in ['EPSG:32630']:
 97    north_loc = np.flip(np.asarray(utm.to_latlon(ts_dav['dataLoc'][0],ts_dav['dataLoc'][1],30,'N')))
 98else:
 99    north_loc = ts_dav['dataLoc']
100north_msl = np.mean(ts_dav['depth'])
101
102#-------------------------------------------
103# Display results for slack water
104#-------------------------------------------
105
106dstrt = dateNumFromDateStr('2023-09-10 00:00:00')
107dstop = dstrt + 1.0
108indx = np.where((dnum >= dstrt)&(dnum <= dstop))[0]
109
110mps2kts = 1.94384
111
112ts = (dnum[indx]-dnum[indx[0]])*24.0
113fig, (ax1,ax2) = plt.subplots(2)
114fig.set_size_inches([10,8])
115ax1.plot(ts,south_vsrf[indx]*mps2kts,'k',label='Surface Velocity Magnitude  [knots]')
116ax1.plot(ts,south_vmag[indx]*mps2kts,'--k',label='Depth Averaged Velocity Magnitude  [knots]')
117ax1.plot(ts,south_elev[indx],'g',label='Surface Elevation  [m]')
118ax1.set_xticks(np.arange(min(ts), max(ts)+1.0, 1.0))
119ax1.set_xlim([0,24])
120ax1.set_ylim([-2,7])
121ax1.grid('on')
122ax1.plot([0,24],[0,0],'k',linewidth=1)
123ax1.plot([0,24],[1,1],'--b',linewidth=1)
124ax1.plot([0,24],[2,2],'--b',linewidth=1)
125ax1.legend(ncols=2,loc='upper center',alignment='center')
126ax1.set_title('Southern Location (' +  
127              '{:+.6f}'.format(south_loc[1]) + 
128              r"$^{\circ}$"+'N, ' + 
129              '{:+.6f}'.format(south_loc[0]) + 
130              r"$^{\circ}$"+'E) , Mean Water Depth = ' + 
131              '{:.1f}'.format(south_msl)+' m')
132
133ax2.plot(ts,north_vsrf[indx]*mps2kts,'k',label='Surface Velocity Magnitude  [knot/s]')
134ax2.plot(ts,north_vmag[indx]*mps2kts,'--k',label='Depth Averaged Velocity Magnitude  [knots]')
135ax2.plot(ts,north_elev[indx],'g',label='Surface Elevation  [m]')
136ax2.set_xticks(np.arange(min(ts), max(ts)+1.0, 1.0))
137ax2.set_xlim([0,24])
138ax2.set_ylim([-2,7])
139ax2.grid('on')
140ax2.plot([0,24],[0,0],'k',linewidth=1)
141ax2.plot([0,24],[1,1],'--b',linewidth=1)
142ax2.plot([0,24],[2,2],'--b',linewidth=1)
143ax2.legend(ncols=2,loc='upper center',alignment='center')
144ax2.set_xlabel(dateStrFromDateNum(dstrt,timeFmt="%d %B %Y")+'  [ hrs ]')
145ax2.set_title('Northern Location (' + 
146              '{:+.6f}'.format(north_loc[1]) + 
147              r"$^{\circ}$"+'N, ' + 
148              '{:+.6f}'.format(north_loc[0]) + 
149              r"$^{\circ}$"+'E) , Mean Water Depth = ' + 
150              '{:.1f}'.format(north_msl)+' m')
151
152fig.show()
153
154print('South Slack Water Times:')
155for t in get_slack_water_times(dnum[indx], south_vmag[indx]):
156    print(t.split(' ')[1]+' UTC')
157print()
158
159lhrange,hlrange,lowTimes,highTimes,ltindx,htindx = get_tidal_range(dnum[indx],south_elev[indx])
160lowElev = south_elev[indx][ltindx]
161highElev = south_elev[indx][htindx]
162print_tide_times(lowTimes,highTimes,lowElev,highElev)
163
164print('North Slack Water Times:')
165for t in get_slack_water_times(dnum[indx], north_vmag[indx]):
166    print(t.split(' ')[1]+' UTC')
167print()
168
169lhrange,hlrange,lowTimes,highTimes,ltindx,htindx = get_tidal_range(dnum[indx],north_elev[indx])
170lowElev = north_elev[indx][ltindx]
171highElev = north_elev[indx][htindx]
172print_tide_times(lowTimes,highTimes,lowElev,highElev)
173
174input('Press ENTER to continue...')
175
176#--------------------------------------------------------------------------
177# INTERFACE
178#--------------------------------------------------------------------------
179

Tidal Analysis

Tidal reduction to a set of harmonics is a standard analysis applied to tidal data. This allows a time series to be represented in a spectral form, which in turn can be used to predict tides at other times. The accuracy of the predicted tides depends on the length of the original timeseries used in the harmonic reduction. The OCMW tools includes a :ref:1tidal <core.tidal>` module that uses UTide to apply harmonic reduction to timeseries of tidal data. The generation of 2D maps of tidal amplitude and phase for key harmonics requires the application of harmonic reduction to every model node.

The script shown in soi2dtideharmslisting demonstrates how tidal reduction can be applied to a multi-day timeseries of surface elevation data. The data generated can be used to map 2D phase and amplitude fields.

Listing 7 soi_free_surface_tidal_analysis.py
 1#!/usr/bin/env python3
 2# -*- coding: utf-8 -*-
 3"""
 4Example of applying tidal reduction to a time series of data extracted from
 5multiple OTM model output files for all nodes.
 6
 7Chris Old
 8IES, School of Engineering, University of Edinburgh
 9Nov 2023
10"""
11# ----------------------------------------------------------------------------
12#   IMPORTS
13# ----------------------------------------------------------------------------
14
15# Standard Python Dependencies
16# Non-Standard Python Dependencies
17import numpy as np
18# Local Module Dependencies
19from ocmw.core.timeFuncs import dateNumFromDateStr
20from ocmw.core.tidal import signalDecomposition, extractHarmValues
21from ocmw.dataproc.otm_extraction import extractMultiFile2DField
22from ocmw.dataproc.processModelData import locDataExists
23from ocmw.dataproc.ocmw_extract import load_ocmw_file, ocmwExtn
24# Other Dependencies
25
26
27# ----------------------------------------------------------------------------
28#   GLOBAL VARIABLES
29# ----------------------------------------------------------------------------
30
31
32# ----------------------------------------------------------------------------
33#   CLASS DEFINITIONS
34# ----------------------------------------------------------------------------
35
36
37# ----------------------------------------------------------------------------
38#   FUNCTION DEFINITIONS
39# ----------------------------------------------------------------------------
40
41
42# ----------------------------------------------------------------------------
43#   MAIN PROCESS
44# ----------------------------------------------------------------------------
45
46#-----------------------------------------------------------------------------
47# Load data
48#-----------------------------------------------------------------------------
49modelPath = 'D:/Models/soi_hr_msl_04/RESULTS'
50dataPath = modelPath.replace('RESULTS','EXTRACT')
51modstr = 'soi_hr_msl_04'
52startstr = '20230317'
53nfiles = 14
54varstr = 'FREE SURFACE'
55
56fileName = modstr+'_'+varstr.replace(' ','')+'_multi_day'+ocmwExtn
57
58if locDataExists(dataPath, fileName):
59    data = load_ocmw_file(dataPath,fileName)
60else:
61    data = extractMultiFile2DField(modelPath,modstr,startstr,nfiles,varstr)
62
63
64#-----------------------------------------------------------------------------
65# Extract time and elevation data
66#-----------------------------------------------------------------------------
67
68dnum = data['times']+dateNumFromDateStr(data['epoch'].replace('/','-'))
69elev = data['FREE_SURFACE']
70
71nNodes = np.max(elev.shape)
72
73#-----------------------------------------------------------------------------
74# Loop over all nodes applying harmonic reduction
75#-----------------------------------------------------------------------------
76
77for n in np.arange(nNodes):
78    elv = elev[:,n]
79    signal = signalDecomposition(dnum,elv)
80    harm, ampl, phase = extractHarmValues(signal['model'])
81    print(harm)
82    print(ampl)
83    print(phase)
84
85#--------------------------------------------------------------------------
86# INTERFACE
87#--------------------------------------------------------------------------
88
89
90

Map 2D model fields

Regional models provide detailed spatial information that can be used to generate 2D maps to visualise spatial variability. The model mesh definition is required to generate 2D maps, and to produce a smooth colourisation the node data should be interpolated onto the element centroids and the centroid value used to colour fill the element. The otm_extraction module includes functions to extract mesh data from a Selafin 2D output file or input geometry file, and the option to extend the mesh data which include the calculation of the element node centroid weighting factors for interpolating node data onto element centroids.

The script shown in Listing 8 gives an example of how to extract data from a 2D output Sealfin file and generate 2D maps of the bathymetry, velocity magnitude, and 2D vorticity field for a single model time step. The base plotting function plotMeshField in the graphics module. This returns handles to the figure, axes and colour bar objects so that the figure can be manually modified.

Listing 8 soi_map_2d_fields
  1#!/usr/bin/env python3
  2# -*- coding: utf-8 -*-
  3"""
  4Exmple of mapping 2D Open Telemac-Mascaret model fields where the data are
  5extracted directly from a Selafin output data file.
  6
  7Chris Old
  8IES, School of Engineering, University of Edinburgh
  9Feb 2024
 10"""
 11
 12# ----------------------------------------------------------------------------
 13#   IMPORTS
 14# ----------------------------------------------------------------------------
 15
 16# Standard Python Dependencies
 17# Non-Standard Python Dependencies
 18import numpy as np
 19# Local Module Dependencies
 20from ocmw.core.graphics import plotMeshField
 21from ocmw.core.physics import circulation2D
 22from ocmw.dataproc.otm_extraction import InitialiseFile, extractMesh
 23# Other Dependencies
 24
 25
 26# ----------------------------------------------------------------------------
 27#   GLOBAL VARIABLES
 28# ----------------------------------------------------------------------------
 29
 30
 31# ----------------------------------------------------------------------------
 32#   CLASS DEFINITIONS
 33# ----------------------------------------------------------------------------
 34
 35
 36# ----------------------------------------------------------------------------
 37#   FUNCTION DEFINITIONS
 38# ----------------------------------------------------------------------------
 39
 40
 41# ----------------------------------------------------------------------------
 42#   MAIN PROCESS
 43# ----------------------------------------------------------------------------
 44
 45#-----------------------------------------
 46# Open the Selafin file
 47#-----------------------------------------
 48datapath = 'D:/Models/SoI/soi_calib/RESULTS'
 49datafile = 'soi_calib_ext_v2_20231029_2D.slf'
 50
 51slf = InitialiseFile(datapath+'/'+datafile)
 52
 53#-----------------------------------------
 54# Extract the model mesh structure
 55#-----------------------------------------
 56mesh = extractMesh(slf,enhance=True,saveMesh=False)
 57
 58#-----------------------------------------
 59# Get data for first timestep
 60#-----------------------------------------
 61fields = slf.get_values(0)
 62varNames = slf.varnames
 63
 64#-----------------------------------------
 65# Plot bathymetry
 66#-----------------------------------------
 67print('bathymetry map')
 68if 'BOTTOM'.ljust(16) in varNames:
 69    indx = varNames.index('BOTTOM'.ljust(16))
 70    bathy = -fields[indx]
 71
 72meshx = mesh['meshx']
 73meshy = mesh['meshy']
 74nElem = mesh['nElems']
 75elems = mesh['elems']
 76cWghts = mesh['cWghts'].copy()
 77
 78print('  mapping node data to element centroids...')
 79bathy_cntr = np.zeros([nElem,],dtype=np.single)
 80for ielm in np.arange(nElem, dtype=np.int64):
 81    cw = cWghts[ielm]
 82    v = bathy[elems[ielm]]
 83    bathy_cntr[ielm] = np.sum((v*cw))
 84
 85print('  generating plot...')
 86fig_bthy = plotMeshField(mesh,bathy_cntr,figsize=(6,10),vmin=0.0,vmax=70,
 87                    titleStr="SoI: Bathymetry",
 88                    unitStr="m", 
 89                    show_cbar=True, 
 90                    cbar_width="4.0%", 
 91                    sub_reg=[303000.0,311000.0,6185000.0,6204000.0],
 92                    ismesh=False)
 93ax = fig_bthy[1]
 94ax.set_xlabel('Eastings  [ m ]')
 95ax.set_ylabel('Northings  [ m ]')
 96
 97#-----------------------------------------
 98# Plot velocity magnitude
 99#-----------------------------------------
100print('velocity magnitude map')
101indx = varNames.index('VELOCITY U'.ljust(16))
102velu = fields[indx]
103indx = varNames.index('VELOCITY V'.ljust(16))
104velv = fields[indx]
105
106vmag = np.sqrt(np.square(velu)+np.square(velv))
107
108print('  mapping node data to element centroids...')
109vmag_cntr = np.zeros([nElem,],dtype=np.single)
110for ielm in np.arange(nElem, dtype=np.int64):
111    cw = cWghts[ielm]
112    v = vmag[elems[ielm]]
113    vmag_cntr[ielm] = np.sum((v*cw))
114
115print('  generating plot...')
116fig_vmag = plotMeshField(mesh,vmag_cntr,figsize=(6,10),vmin=0.0,vmax=3.5,
117                    titleStr="SoI: Velocity Magnitude",
118                    unitStr="m/s", 
119                    show_cbar=True, 
120                    cbar_width="4.0%", 
121                    sub_reg=[303000.0,311000.0,6185000.0,6204000.0],
122                    ismesh=False)
123ax = fig_vmag[1]
124ax.set_xlabel('Eastings  [ m ]')
125ax.set_ylabel('Northings  [ m ]')
126
127#-----------------------------------------
128# Plot vorticity
129#-----------------------------------------
130print('vorticity map')
131print('  calculating vorticity...')
132vort = np.zeros([nElem,],dtype=np.single)
133
134for ielm in np.arange(nElem, dtype=np.int64):
135    x = meshx[elems[ielm,:]].tolist()
136    y = meshy[elems[ielm,:]].tolist()
137    vx = velu[elems[ielm,:]].tolist()
138    vy = velv[elems[ielm,:]].tolist()
139    vort[ielm] = circulation2D(x, y, vx, vy)/mesh['area'][ielm]
140
141print('  generating plot...')
142fig_vort = plotMeshField(mesh,vort,figsize=(6,10),vmin=-0.02,vmax=0.02,
143                    titleStr="SoI: 2D Vorticity",
144                    unitStr="$s^{-1}$", 
145                    show_cbar=True, 
146                    cbar_width="4.0%", 
147                    sub_reg=[303000.0,311000.0,6185000.0,6204000.0],
148                    ismesh=False)
149ax = fig_vort[1]
150ax.set_xlabel('Eastings  [ m ]')
151ax.set_ylabel('Northings  [ m ]')
152
153
154print('Process complete.\n')
155input('Press ENTER to exit...')
156
157
158#--------------------------------------------------------------------------
159# INTERFACE
160#--------------------------------------------------------------------------
161

Model Validation Against ADCP Data

The OCMW tools include a module for validataion of model velocity and elevation data against in situ ADCP data . The input data are co-located and contemporaneous timeseries of model and in situ dat, as generated using the otm_extract_by_location.py application and the processModelData module functions. The observation data are time interpoltated onto the model data timestamps, then a set of standard validatiaon statistics are calculated and saved to the ANALYSIS directory, and summarised in a graphic including a correlation plot.

The script shown in Listing 9 presents an example of validation of a Falls of warness, Orkneys model against ReDAPT ADCP data. in this example the PWRA velocity for a DeepGen IV turbine are used as the validation input data. These are calculated from both the model and the ADCP data. Validation statistics are calculated for the velocity magnitude, and graphics displayed.

Listing 9 fow_model_flow_validation.py
 1#!/usr/bin/env python3
 2# -*- coding: utf-8 -*-
 3"""
 4Example Script:  Validation of the FASTWATER FoW models against ReDAPT ADCP data
 5
 6
 7"""
 8
 9# ----------------------------------------------------------------------------
10#   IMPORTS
11# ----------------------------------------------------------------------------
12
13# Standard Python Dependencies
14# Non-Standard Python Dependencies
15import numpy as np
16# Local Module Dependencies
17from ocmw.core.tecSpecs import getTec
18from ocmw.dataproc.processModelData import hub_height_data
19from ocmw.calval.flowValidation import validateFlowParameter
20# Other Dependencies
21
22
23# ----------------------------------------------------------------------------
24#   GLOBAL VARIABLES
25# ----------------------------------------------------------------------------
26
27
28# ----------------------------------------------------------------------------
29#   CLASS DEFINITIONS
30# ----------------------------------------------------------------------------
31
32
33# ----------------------------------------------------------------------------
34#   FUNCTION DEFINITIONS
35# ----------------------------------------------------------------------------
36
37
38# ----------------------------------------------------------------------------
39#   MAIN PROCESS
40# ----------------------------------------------------------------------------
41# Define turbine to model
42tec = getTec('DeepGenIV')
43
44if tec.bottomMounted:
45    rstr = '_hab_'
46else:
47    rstr = '_dbs_'
48
49
50# Model data
51modDPath = 'D:/Models/fw_fow_case_study_v02'
52modFName = 'fow_vort_model_merged_td7_01.mat'
53modDataSrc = 'FoW Vorticity Refined Model'
54modDataSet = 'TD7_01'
55
56# Extract hub height PWRA data
57ts_mod, data_mod, pwra_mod = hub_height_data(modDPath, modFName,
58                                             modDataSrc, modDataSet, tec,
59                                             modFName)
60
61zhub = tec.hubHeight
62ofprts = modFName.split('.')
63pwrafile = ofprts[0]+rstr+str(np.int32(np.ceil(zhub))).zfill(2)+'m_PWRA.'+ofprts[1]
64modFile = modDPath+'/ANALYSIS/'+pwrafile
65
66# In situ ADCP data
67obsDPath = 'D:/Data/ReDAPT/NETCDF/ADCPTD7_01_Dep1'
68obsFName = 'ADCPTD7_01_Dep01_merged.mat'
69obsDataSrc = 'ReDAPT ADCP archive data'
70obsDataSet = 'ADCPTD7_01_Dep1'
71
72ts_obs, data_obs, pwra_obs = hub_height_data(obsDPath, obsFName,
73                                             obsDataSrc, obsDataSet, tec,
74                                             obsFName)
75
76ofprts = obsFName.split('.')
77pwrafile = ofprts[0]+rstr+str(np.int32(np.ceil(zhub))).zfill(2)+'m_PWRA.'+ofprts[1]
78obsFile = obsDPath+'/ANALYSIS/'+pwrafile
79
80
81# Calculate Validataion Statistics
82metrics, figs, info = validateFlowParameter(modFile,obsFile,['magnitude'],doFlowChar=False,varStr='PWRA Velocity')
83
84for fig in figs:
85    fig.show()
86
87print('Process complete.\n')
88input('Press ENTER to exit...')
89
90
91#--------------------------------------------------------------------------
92# INTERFACE
93#--------------------------------------------------------------------------
94

Convert ADCP data to netCDF

For integration with the OCMW tools and data archiving, raw ADCP data are converted into a standardised netCDF file format. The data are sorted into daily files to produce managable data fragments and to simplify data processing and archiving. The tools currrently support RDI Workhorse and Nortek Signature instruments. There are tools in-place to read Sonardyne Origin data, at some point conversion to netCDF will be included.

The script shown in Listing 10 shows the conversion of Nortek Signature ADCP data into daily netCDF files. The standardised netCDF files store depoyment and processing metadata in the global attributes for the netCDF file. The global attributes can be used to catalogue the files for databasing. Some of the metadata needs to be manually added, this done with the internal function populateMetadata. The remaining fields are populated using data from the raw ADCP files. A list of raw files is generated and sorted, then sequentially processed to extract the data into daily netCDF files. The functions for doing this are in the adcp2netcdf module.

Listing 10 convert_sig_to_netcdf
  1#!/usr/bin/env python3
  2# -*- coding: utf-8 -*-
  3"""
  4Created on Fri Oct 29 12:01:47 2021
  5
  6@author: cold2
  7"""
  8
  9# ----------------------------------------------------------------------------
 10#   IMPORTS
 11# ----------------------------------------------------------------------------
 12
 13# Standard Python Dependencies
 14import os
 15# Non-Standard Python Dependencies
 16import numpy as np
 17import utm
 18# Local Module Dependencies
 19from ocmw.core.fileTools import getListOfFiles, loadmat
 20from ocmw.core.timeFuncs import dateStrFromDateNum, matlab2std
 21from ocmw.dataman.adcp2netcdf import metaDataObj, build_adcp_netcdf
 22from ocmw.dataman.adcp2netcdf import populate_sig_core, populate_sig_fields
 23# Other Dependencies
 24
 25
 26# ----------------------------------------------------------------------------
 27#   GLOBAL VARIABLES
 28# ----------------------------------------------------------------------------
 29
 30
 31# ----------------------------------------------------------------------------
 32#   CLASS DEFINITIONS
 33# ----------------------------------------------------------------------------
 34
 35
 36# ----------------------------------------------------------------------------
 37#   FUNCTION DEFINITIONS
 38# ----------------------------------------------------------------------------
 39def populateMetadata():
 40    metadata = metaDataObj()
 41    metadata.title = 'RealTide UEDIN SIG1000 standardization'
 42    metadata.projects = 'RealTide'
 43    metadata.references = ''
 44    metadata.site = 'Fromveur Strait, France'
 45    metadata.deploy_id= ''
 46    metadata.deploy_mooring_type = 'Concrete bed frame'
 47    metadata.deploy_vessel = 'Olympic Challenger'
 48    metadata.deploy_vessel_operator = 'Olympic Subsea ASA'
 49    metadata.geospatial_easting = 0.0
 50    metadata.geospatial_northing = 0.0
 51    metadata.geospatial_easting = 0.0
 52    metadata.geospatial_northing = 0.0
 53    metadata.geospatial_crs = ''
 54    metadata.data_type = 'field measurement'
 55    metadata.processing_level = '1'
 56    metadata.quality_control_indicator = '1'
 57    
 58    # Instrument
 59    metadata.instrument_owner = 'Universidade do Algarve'
 60    metadata.instrument_maker = 'Nortek'
 61    metadata.instrument_type = 'Signature'
 62    metadata.instrument_freq = '1000kHz'
 63    metadata.instrument_serial_no = ''
 64    metadata.instrument_firmware_ver = ''
 65    metadata.instrument_head_hab = 0.0
 66    
 67    # File generation
 68    metadata.contact = ''
 69    metadata.institution = 'IES, School of Engineering, The University of Edinburgh, UK'
 70    metadata.data_assembly_center = ''
 71    metadata.date_created = ''
 72    metadata.date_updated = ''
 73    metadata.history = ''
 74    metadata.naming_authority = 'University of Edinburgh'
 75    metadata.data_format_version = '1.0'
 76    metadata.conventions = 'CF-1.8'
 77    metadata.netcdf_version = 'netCDF-4 classic model'
 78    metadata.license = ''
 79    metadata.citation = ''
 80    metadata.doi = ''
 81    return metadata
 82
 83# ----------------------------------------------------------------------------
 84#   MAIN PROCESS
 85# ----------------------------------------------------------------------------
 86
 87deploy = 'UEdin_Sig1000_Nov5'
 88
 89deploy_utm = [349601.7, 5368063.6]
 90deploy_loc = utm.to_latlon(deploy_utm[0],deploy_utm[1],30,'N')
 91
 92dpath = 'D:/Data/Fromveur_SIG1000_2019/measurements_SIG1000_2019/Converted'
 93files = getListOfFiles(dpath,'*.mat')
 94
 95metadata = populateMetadata()
 96
 97genfields = ['beam','ampl','corr','inst','earth','vert']
 98loadfields = ['beam','ampl','corr','inst','earth']
 99
100ncpath = 'D:/Data/Fromveur_SIG1000_2019/measurements_SIG1000_2019/NETCDF'
101findx = 0
102
103# open source deployment data file
104data_file = os.path.join(dpath,files[0]).replace('\\','/')
105data = loadmat(data_file)
106cfg = data['Config']
107# extract data dimensions
108num_beams = cfg['Burst_NBeams']
109num_bins = cfg['Burst_NCells']
110# set metadata
111metadata.instrument_serial_no = cfg['SerialNo']
112metadata.instrument_firmware_ver = str(cfg['FWVersion']/1000)
113metadata.instrument_beam_angle = 25.0
114metadata.instrument_head_hab = 0.8
115metadata.geospatial_easting = deploy_utm[0]
116metadata.geospatial_northing = deploy_utm[0]
117metadata.geospatial_latitude = deploy_loc[0]
118metadata.geospatial_longitude = deploy_loc[1]
119metadata.geospatial_crs = 'WGS84 / UTM Zone 30N'
120
121del cfg
122
123
124for file in files:
125
126    print(file)
127    dfile = os.path.join(dpath,file).replace('\\','/')
128    src = loadmat(dfile)
129    dn0 = np.double(np.floor(matlab2std(src['Data']['Burst_Time'][0])))
130    dn1 = np.double(np.floor(matlab2std(src['Data']['Burst_Time'][-1])))
131    
132    print(dn0,dn1)
133    
134    t = matlab2std(src['Data']['Burst_Time'])
135    
136    ndays = np.int64(dn1-dn0)+1
137    
138    for aday in range(ndays):
139        
140        print(aday)
141
142        tstr = dateStrFromDateNum(dn0+aday,timeFmt='%Y%m%d')
143        ncfile = deploy+'_'+tstr+'.nc'
144        print(ncfile)
145    
146        if not os.path.isfile(ncpath+'/'+ncfile):
147            build_adcp_netcdf(ncpath,ncfile,genfields,num_beams,num_bins)
148            populate_sig_core(dpath,file,ncpath,ncfile,metadata)
149
150        indices = np.where( (t >= dn0+aday) & (t < dn0+aday+1.0) )[0]
151        populate_sig_fields(dpath,dfile,ncpath,ncfile,loadfields,indices,metadata)
152
153
154#--------------------------------------------------------------------------
155# INTERFACE
156#--------------------------------------------------------------------------