#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Adcp instruments (:mod:`hydrac.instruments.adcp`)
=================================================
.. autoclass:: Adcp
:members:
:private-members:
.. autoclass:: ReadAdcpA
:members:
:private-members:
.. autoclass:: ReadAdcpP
:members:
:private-members:
"""
# from .instruments import instrument
from .acoustique import Acoustic
from hydrac.util import constants
import numpy as np
import struct as st
import os
import glob
import datetime
import time as ti
from tkinter import filedialog
from tkinter import Tk
from hydrac.util.parameters import Parameters
from hydrac.util.paramcontainer import ParamContainer
# %gui inline
# a=hydrac.Instruments.aquascat.aquascat('nnn','lll')
# a.param=hydrac.aquascat.aquascat.batch_read(a.filepath,'',a.param)
# test=a.param
# test[a.filepath[0]]['AbsRxFrequency']
[docs]class Adcp(Acoustic):
""" Adcp instruments class
Base class : :mod:`hydrac.instruments.acoustique.Acoustic`
The ADCP is a single frequency Doppler velocity profiler. It measures the
Doppler phase shifts induced by the moving targets (suspended particles) in
the water column. As well as measuring the multi-component velocity, it is
also capable of recording the backscattered intensity, which itself is of
interest in hydrac.
.. image:: ../_image/adcp.pdf
:width: 100%
:align: center
The Adcp class reads the raw binary files from ADCPs instruments and
stores the valuable information into the ``param`` modified dictionnary
with the common shape handled by hydrac
(see :mod:`hydrac.instruments.acoustique`).
The way the data are to be considered is similar to the example shown in
:mod:`hydrac.instruments.aquascat`.
Parameters
----------
name : str, {'Campaign_1'}
Tag name, for personal usage, for instance to add a campaign name."""
def __init__(self, name):
Acoustic.__init__(self)
self.name = name
self.instr_name = 'adcp'
self.adcp_presets = Parameters(constants.ADCP)
self.tag_M = 'mono'
self.filepath, self.tempdir, self.tag_adcpType = self.file_select()
self.load_ADCP()
# root = Tk()
# filez = filedialog.askopenfilenames(parent=root,title='Choose a file')
# root.withdraw() # use to hide tkinter window
# currdir = os.getcwd()
## self.tempdir = filedialog.askdirectory(initialdir=currdir, title='Please select a directory')
# self.filepath = filez[0]
# self.tempdir = os.path.dirname(os.path.abspath(filez[0]))
# self.tag_adcpType=os.path.split(os.path.abspath(filez[0]))[1][-5]
# if self.tag_adcpType=='p':
# read_ADCP_P(self)
# else:
# read_ADCP_A(self)
# print(str(self.filepath))
[docs] def _set_calib_adcp(self):
"""Calibration coefficients for each transducer are unique to each
transducer. As such, each instrument comes with a different set of
calibration coefficients Kt. Those are stored in xml file generated
(if non-existant) or updated (with :func:`_update_calib_aqa`) with the
instrument serial number and the Kts. This xml is accessed/generated
through the class :class:`hydrac.util.paramcontainer.ParamContainer`.
"""
rr = glob.glob(constants.dir_adcp_calib + '*.xml')
rr.sort()
try:
adcp = Parameters(ParamContainer(path_file=rr[-1]))
except(IndexError, AttributeError):
adcp = Parameters({})
return adcp
[docs] def _update_calib_adcp(self, a, b=None, remove_elem=False):
"""Update the calibration xml file either by adding a new instrument
or deleting the coefficients of a particular instrument. This xml is
accessed/generated through the class
:class:`hydrac.util.paramcontainer.ParamContainer`.
Parameters
----------
a : str
Instrument serial number read from binary files and stored in
``param``
b=None : dict
Contains the calibration coefficients, prompted by the user if the
instrument is used for the first time in hydrac.
remove_elem=False : bool
Removes the information relative to the serial number a in the xml.
"""
rr = glob.glob(constants.dir_adcp_calib + '*.xml')
rr.sort()
if remove_elem is False:
uiui = ParamContainer(path_file=rr[-1])
uiui._set_child(a, b)
uiui._save_as_xml(path_file=constants.dir_adcp_calib +
'adcp_calib.xml', find_new_name=False,
overwrite=True)
[docs] def file_select(self):
""" User input selection of the target files to read """
root = Tk()
filez = filedialog.askopenfilenames(parent=root, title='Choose a file')
root.withdraw() # use to hide tkinter window
filepath = list(filez)
filepath.sort()
filepath = tuple(filepath)
tempdir = os.path.dirname(os.path.abspath(filez[0]))
tag_adcpType = [os.path.split(os.path.abspath(filez[i]))[1][-5]
for i in range(0, len(filez))]
return filepath, tempdir, tag_adcpType
[docs] def param_shape_ADCP(self, inp):
"""Reshapes the modified dictionnary ``param`` to fit the standard key
names and data shape used throughout hydrac
"""
FNAME = ['pings_per_ensemble', 'cell_size', 'bin1_dist', 'n_cells',
'n_beams', 'beam_freq', 'TVG', 'intens', 'range',
'beam_angle', 'depth', 'temperature', 'salinity']
FNAME_out = ['NumAverage', 'BinLengthMM', 'StartBin', 'NumBins',
'NumChannels', 'Frequency', 'TVG', 'RawIntensity',
'BinRange', 'TransducerAngle', 'depth', 'temp', 'sal']
FNAME_meta = ['coord_sys', 'orientation', 'name', 'prof_mode',
'prog_ver']
inp2 = Parameters({'adcp': {}})
inpm = Parameters({})
[inp2.update({FNAME_out[tui]:inp[FNAME[tui]]})
if FNAME[tui] in inp.keys()
else inp2.update({FNAME_out[tui]:inp['config'][FNAME[tui]]})
for tui in range(0, len(FNAME))]
[inpm.update({FNAME_meta[tui]:inp[FNAME_meta[tui]]})
if FNAME_meta[tui] in inp.keys()
else inpm.update({FNAME_meta[tui]:inp['config'][FNAME_meta[tui]]})
for tui in range(0, len(FNAME_meta))]
inp2.FileName = os.path.split(self.currfilepath)[-1]
if self.tag_adcpType[self.Ncurrfilepath] == 'p':
FNAME_aux = ['intens_scale', 'absorption', 'avg_interval',
'avg_method', 'time_between_ping_groups',
'beam_pattern', 'blank', 'pitch', 'bt_x_disp',
'bt_y_disp', 'bt_east_vel', 'bt_error_vel',
'bt_east_vel', 'bt_north_vel', 'bt_vert_vel',
'nav_east_vel', 'nav_north_vel', 'latitude',
'longitude', 'bt_range', 'pitch', 'north_vel',
'heading', 'error_vel', 'east_vel', 'roll',
'vert_vel']
FNAME_meta = ['firm_ver', 'adcp_type', 'deploy_name',
'use_pitchroll']
inp2.time = np.array([datetime.datetime.timestamp(
datetime.datetime.fromordinal(int(inp['mtime'][i])-366) +
datetime.timedelta(days=inp['mtime'][i] % 1))
for i in range(0, len(inp['mtime']))])
inp2.TransducerKt = np.zeros((inp2.NumChannels))
inp2.PulseLength = np.zeros((inp2.NumChannels))
inp2.pressure = np.zeros((len(inp2.time), inp2.NumChannels))
inp2.depth = np.tile(inp2.depth, (len(inp2.BinRange),
inp2.NumChannels))
[inp2.adcp.update({FNAME_aux[tui]:inp[FNAME_aux[tui]]})
if FNAME_aux[tui] in inp.keys()
else inp2.adcp.update({FNAME_aux[tui]:inp['config'][FNAME_aux[tui]
]
}) for tui in range(0, len(FNAME_aux))]
[inpm.update({FNAME_meta[tui]: inp[FNAME_meta[tui]]})
if FNAME_meta[tui] in inp.keys()
else inpm.update({FNAME_meta[tui]: inp['config'][FNAME_meta[tui]]
}) for tui in range(0, len(FNAME_meta))]
else:
FNAME_aux = ['time_between_ping_groups', 'beam_pattern', 'blank',
'corr_threshold', 'fls_target_threshold', 'min_pgood',
'perc_good', 'pitch', 'pitch_std', 'sysbandwidth',
'bt_corr', 'bt_perc_good', 'bt_range', 'bt_vel',
'pitch_std', 'pitch', 'perc_good', 'north_vel',
'heading_std', 'heading', 'error_vel', 'east_vel',
'corr_raw', 'corr', 'bt_vel', 'roll', 'roll_std',
'vert_vel']
inp2.time = inp.rmtime
inp2.TransducerKt = np.zeros((inp2.NumChannels))
inp2.PulseLength = np.zeros((inp2.NumChannels))
[inp2.adcp.update({FNAME_aux[tui]: inp[FNAME_aux[tui]]})
if FNAME_aux[tui] in inp.keys()
else inp2.adcp.update({FNAME_aux[tui]:inp['config'][FNAME_aux[tui]
]
}) for tui in range(0, len(FNAME_aux))]
inpm.adcp_type = 'broadband'
inpm.deploy_name = ' '
inpm.serialnum = inp['config']['serialnum']
inp2.pressure = np.tile(inp['pressure'], (inp2.NumChannels, 1)).T
call_calib = self._set_calib_adcp()
inp2.adcp.adcp_presets = call_calib.RDI['f' + str(inp2.Frequency)]
inp2.TransducerRadius = np.zeros((inp2.NumChannels)) +\
inp2.adcp.adcp_presets.TransducerRadius
inp2.TransducerBeamWidth = np.zeros((inp2.NumChannels)) +\
inp2.adcp.adcp_presets.TransducerBeamWidth
inp2.TransducerAngle = np.tile(inp2.TransducerAngle,
(inp2.NumChannels))
inp2.RxGain = np.tile(0, (inp2.NumChannels))
inp2.BinLengthMM = np.tile(inp2.BinLengthMM, (inp2.NumChannels))
inp2.TxGain = np.tile(0, (inp2.NumChannels))
inp2.StartBin = np.tile(inp2.StartBin, (inp2.NumChannels))
inp2.BurstTime = inp2.time[0]
inp2.BinRange = np.tile(inp2.BinRange, (inp2.NumChannels, 1))
inp2.BinRange = inp2.BinRange.T
inp2.NumAverage = np.tile(inp2.NumAverage, (inp2.NumChannels))
inp2.NumBins = np.tile(inp2.NumBins, (inp2.NumChannels))
inp2.PingRate = 1 / (inp2.time[1] - inp2.time[0])
inp2.SampleRate = np.tile(1 / (inp2.time[1] - inp2.time[0]),
(inp2.NumChannels))
inp2.NumPingsTot = len(inp2.time)
inp2.TransducerTag = np.array(range(0, inp2.NumChannels)).T
inp2.TransducerName = np.tile(np.array(str(inp2.Frequency) +
' MHz', dtype=str),
(inp2.NumChannels))
inp2.NumProfilesSamples = np.tile(len(inp2.time), (inp2.NumChannels))
inp2.Frequency = np.tile(inp2.Frequency, (inp2.NumChannels)) * 1000
inp2.FrequencyRx = np.zeros(np.size((inp2.Frequency))) +\
inp2.adcp.adcp_presets.Frequency
inp2.MeanProfile = np.nanmean(inp2.RawIntensity,
axis=(np.where(np.array(
np.shape(inp2.RawIntensity
)) == len(inp2.time))
)[0][0])
inp2.battery = np.zeros(np.shape(inp2.temp))
inp2.NumSamples = inp2.NumPingsTot / inp2.NumAverage
inp2.RawIntensity = np.power(10, ((np.transpose(inp2.RawIntensity,
(0, 2, 1)) -
inp2.adcp.adcp_presets.EC0) *
inp2.adcp.adcp_presets.Kc +
inp2.adcp.adcp_presets.B -
inp2.adcp.adcp_presets.SL0)/20)
inp2.TransducerKt = [[a, b] for a, b in zip(list(inp2.FrequencyRx.T),
list(np.sqrt(3 / (16 *
np.pi)) +
np.zeros((inp2.NumChannels
)).T))]
# Ici AuxData ce sont les vitesses et tout et tout
inp2.adcp.time_aux = inp2.time
inp2.adcp.AuxNumSamples = inp2.NumSamples
inp2.adcp.AuxSampleRate = inp2.SampleRate
inp2.adcp.NumAuxChannels = len(inp2.adcp)
return inp2, inpm
[docs] def load_ADCP(self):
""" Launch decryption of the target files -
Two possibilities given the file extension :
ReadAdcpA (standard format) or ReadAdcpP (standard +
GPS tracking"""
if hasattr(self, 'filepath'):
for kk in range(len(self.filepath)):
if self.tag_adcpType[kk] == 'p':
self.Ncurrfilepath = kk
self.currfilepath = self.filepath[kk]
ReadAdcpP(self)
else:
self.Ncurrfilepath = kk
self.currfilepath = self.filepath[kk]
ReadAdcpA(self)
else:
raise AttributeError('No file selected')
class WrongFileTypeError(Exception):
pass
[docs]class ReadAdcpA(object):
""" Reads standard RDI fileformat.
Parameters
----------
obj : Adcp object :class:`hydrac.instruments.ascp.Adcp`
"""
def __init__(self, obj):
self.obj = obj
if self.obj.tag_adcpType[self.obj.Ncurrfilepath] == 'p':
raise WrongFileTypeError('The file extension is not adapted...')
adcp, metadcp = self.obj.param_shape_ADCP(self.read_ADCP(
self.obj.currfilepath, ''))
self.obj.param.update({"P" + str(self.obj.Ncurrfilepath): adcp})
self.obj.meta.update({"P" + str(self.obj.Ncurrfilepath): metadcp})
def checkheader(self, fd):
valid = 0
# Following the header bytes is numbytes and we move forward numbytes>0
numbytes, = st.unpack("h", fd.read(2))
if numbytes > 0:
try:
fd.seek(numbytes-2, 1)
tmp = [fd.read(1) for ui in range(1, 2+1)]
cfgid = np.array([st.unpack("B", tmp[ui])[0]
for ui in range(0, 1+1) if tmp[ui] != b''])
# will Skip the last ensemble (sloppy code)
if len(cfgid) == 2:
fd.seek(-numbytes-2, 1)
# and we have *another* 7F7F
if cfgid[0] == int('7F', 16) & cfgid[1] == int('7F', 16):
# ...ONLY THEN it is valid.
valid = 1
except IOError:
print("Crotte !")
fd.seek(-2, 1)
return valid
# -------------------------------------
def rd_hdr(self, fd):
# Read config data
# Changed by Matt Brennan to skip weird stuff at BOF (apparently
# can happen when switching from one flash card to another
# in moored ADCPs).
cfgid = np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 2+1)])
nread = 0
while (cfgid[0] != int('7F',16)) or (cfgid[1] != int('7F',16)) or (self.checkheader(fd) == 0):
nextbyte, = st.unpack("B", fd.read(1))
pos = fd.tell()
nread = nread + 1
# End of file
if not nextbyte:
print('EOF reached before finding valid cfgid')
hdr = np.NaN
break
cfgid[1] = cfgid[0]
cfgid[0] = nextbyte
if (pos % 1000) == 0:
print(('Still looking for valid cfgid at file position {} ...'
).format(pos))
pos = fd.tell() - 2
if nread > 0:
print(('Junk found at BOF...skipping {} bytes until '
).format(nread))
print(('cfgid={} {} at file pos {}'
).format(hex(cfgid[0]), hex(cfgid[1]), pos))
hdr, nbyte = self.rd_hdrseg(fd)
return hdr, pos
# -------------------------------------
def rd_fix(self, fd):
# Read config data
cfgid, = st.unpack("H", fd.read(2))
if cfgid != int('0000', 16):
print(('Fixed header ID incorrect - ' +
'data corrupted or not a BB/WH raw file?').format(cfgid))
cfg, nbyte = self.rd_fixseg(fd)
return cfg
# --------------------------------------
def rd_hdrseg(self, fd):
# Reads a Header
if 'hdr' is not locals() or 'hdr' is not globals():
hdr = {}
hdr.update({'nbyte': st.unpack("h", fd.read(2))[0]})
fd.seek(1, 1)
ndat, = st.unpack("b", fd.read(1))
hdr.update({'dat_offsets': np.array([st.unpack("h", fd.read(2))[0]
for ui in range(1, ndat+1)])})
nbyte = 4 + ndat * 2
return hdr, nbyte
# -------------------------------------
def getopt(self, val, *varargin):
# Returns one of a list (0=first in varargin, etc.)
if val > len(varargin):
opt = 'unknown'
else:
opt = varargin[val]
return opt
# -------------------------------------
def rd_fixseg(self, fd):
# Reads the configuration data from the fixed leader
# #disp(fread(fd,10,'uint8'))
# #fseek(fd,-10,'cof');
if 'cfg' is not locals() or 'cfg' is not globals():
cfg = {}
cfg.update({'name': 'wh-adcp'})
# default - depending on what data blocks are around
# we can modify this later in rd_buffer.
cfg.update({'sourceprog': 'instrument'})
cfg.update({'prog_ver': st.unpack("B", fd.read(1))[0] +
st.unpack("B", fd.read(1))[0] / 100})
# 8,9,16 - WH navigator
# 10 -rio grande
# 15, 17 - NB
# 19 - REMUS, or customer specific
# 11- H-ADCP
# 31 - Streampro
# 34 - NEMO
# 50 - WH, no bottom track (built on 16.31)
# 51 - WH, w/ bottom track
# 52 - WH, mariner
if np.fix(cfg['prog_ver']) == 4 or np.fix(cfg['prog_ver']) == 5:
cfg.update({'name': 'bb-adcp'})
elif np.fix(cfg['prog_ver']) == 8\
or np.fix(cfg['prog_ver']) == 9\
or np.fix(cfg['prog_ver']) == 10\
or np.fix(cfg['prog_ver']) == 16\
or np.fix(cfg['prog_ver']) == 50\
or np.fix(cfg['prog_ver']) == 51\
or np.fix(cfg['prog_ver']) == 52:
cfg.update({'name': 'wh-adcp'})
# phase 1 and phase 2
elif np.fix(cfg['prog_ver']) == 14 or np.fix(cfg['prog_ver']) == 23:
cfg.update({'name': 'os-adcp'})
else:
cfg.update({'name': 'unrecognized firmware version'})
# Coded stuff
config = np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 2+1)])
cfg.update({'config': np.base_repr(config[1], 2, 8) +
'-' +
np.base_repr(config[0], 2, 8)})
cfg.update({'beam_angle': self.getopt((config[1] & 3), 15, 20, 30)})
tag_ = [1 if x else 0 for x in [(config[0] & 16) == 16,
(config[0] & 8) == 8,
(config[0] & 128) == 128]]
cfg.update({'numbeams': self.getopt(tag_[0], 4, 5)})
cfg.update({'beam_freq': self.getopt((config[0] & 7), 75, 150, 300,
600, 1200, 2400, 38)})
# 1=convex, 0=concave
cfg.update({'beam_pattern': self.getopt(tag_[1], 'concave', 'convex')})
# 1=up,0=down
cfg.update({'orientation': self.getopt(tag_[2], 'down', 'up')})
# Flag for simulated data
cfg.update({'simflag': self.getopt(st.unpack("B", fd.read(1))[0],
'real',
'simulated')})
fd.seek(1, 1)
cfg.update({'n_beams': st.unpack("B", fd.read(1))[0]})
cfg.update({'n_cells': st.unpack("B", fd.read(1))[0]})
cfg.update({'pings_per_ensemble': st.unpack("H", fd.read(2))[0]})
cfg.update({'cell_size': st.unpack("H", fd.read(2))[0] * .01}) # meters
cfg.update({'blank': st.unpack("H", fd.read(2))[0] * .01}) # meters
cfg.update({'prof_mode': st.unpack("B", fd.read(1))[0]})
cfg.update({'corr_threshold': st.unpack("B", fd.read(1))[0]})
cfg.update({'n_codereps': st.unpack("B", fd.read(1))[0]})
cfg.update({'min_pgood': st.unpack("B", fd.read(1))[0]})
cfg.update({'evel_threshold': st.unpack("H", fd.read(2))[0]})
cfg.update({'time_between_ping_groups':
sum(np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 3+1)]) * np.array([60,
1,
.01]))
})
# Lots of bit-mapped info
coord_sys = st.unpack("B", fd.read(1))[0]
cfg.update({'coord': np.base_repr(coord_sys, 2, 8)})
tag_ = [1 if x else 0 for x in [(coord_sys & 4) == 4,
(coord_sys & 2) == 2,
(coord_sys & 1) == 1]]
cfg.update({'coord_sys': self.getopt(((coord_sys >> 3) & 3),
'beam', 'instrument',
'ship', 'earth')})
cfg.update({'use_pitchroll': self.getopt(tag_[0], 'no', 'yes')})
cfg.update({'use_3beam': self.getopt(tag_[1], 'no', 'yes')})
cfg.update({'bin_mapping': self.getopt(tag_[2], 'no', 'yes')})
# degrees
cfg.update({'xducer_misalign': st.unpack("h", fd.read(2))[0] * .01})
# degrees
cfg.update({'magnetic_var': st.unpack("h", fd.read(2))[0] * .01})
cfg.update({'sensors_src': np.base_repr(st.unpack("B", fd.read(1))[0],
2, 8)})
cfg.update({'sensors_avail': np.base_repr(st.unpack("B",
fd.read(1)
)[0], 2, 8)})
# meters
cfg.update({'bin1_dist': st.unpack("H", fd.read(2))[0] * .01})
# meters
cfg.update({'xmit_pulse': st.unpack("H", fd.read(2))[0] * .01})
cfg.update({'water_ref_cells':
np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 2+1)])})
cfg.update({'fls_target_threshold': st.unpack("B", fd.read(1))[0]})
fd.seek(1, 1)
cfg.update({'xmit_lag': st.unpack("H", fd.read(2))[0]*.01}) # meters
nbyte = 40
if np.fix(cfg['prog_ver']) == 8\
or np.fix(cfg['prog_ver']) == 10\
or np.fix(cfg['prog_ver']) == 16\
or np.fix(cfg['prog_ver']) == 50\
or np.fix(cfg['prog_ver']) == 51\
or np.fix(cfg['prog_ver']) == 52:
# Added CPU serial number with v8.14
if cfg['prog_ver'] >= 8.14:
cfg.update({'serialnum': np.array([st.unpack("B",
fd.read(1))[0]
for ui in range(1, 8+1)])})
nbyte = nbyte+8
# Added 2 more bytes with v8.24 firmware
if cfg['prog_ver'] >= 8.24:
cfg.update({'sysbandwidth': np.array([st.unpack("B",
fd.read(1))[0]
for ui in range(1, 2+1)
])})
nbyte = nbyte+2
# Added 1 more bytes with v16.05 firmware
if cfg['prog_ver'] >= 16.05:
cfg.update({'syspower': st.unpack("B", fd.read(1))[0]})
nbyte = nbyte+1
# Added bytes for REMUS, navigators, and HADCP
if cfg['prog_ver'] >= 16.27:
cfg.update({'navigator_basefreqindex':
st.unpack("B", fd.read(1))[0]})
nbyte = nbyte+1
cfg.update({'remus_serialnum':
np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 4+1)])})
nbyte = nbyte+4
cfg.update({'h_adcp_beam_angle': st.unpack("B",
fd.read(1))[0]})
nbyte = nbyte+1
elif np.fix(cfg['prog_ver']) == 9:
# Added CPU serial number with v8.14
if cfg['prog_ver'] >= 9.10:
cfg.update({'serialnum':
np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 8+1)])})
nbyte = nbyte+8
cfg.update({'sysbandwidth':
np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 2+1)])})
nbyte = nbyte+2
elif np.fix(cfg['prog_ver']) == 14 or np.fix(cfg['prog_ver']) == 23:
cfg.update({'serialnum':
np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 8+1)])})
nbyte = nbyte+8
# It is useful to have this precomputed.
cfg.update({'range': cfg['bin1_dist'] +
np.arange(0, cfg['n_cells']) * cfg['cell_size']})
if cfg['orientation'] == 1:
cfg['range'] = -cfg['range']
return cfg, nbyte
def find_all(self, a_str, sub):
start = 0
while True:
start = a_str.find(sub, start)
if start == -1:
return
yield start
# use start += 1 to find overlapping matches
start += len(sub)
# -----------------------------
def rd_buffer(self, fd, num_av):
# To save it being re-initialized every time.
# global ens hdr
# A fudge to try and read files not handled quite right.
# global FIXOFFSET self.SOURCE
# If num_av<0 we are reading only 1 element and initializing
if num_av < 0:
self.SOURCE = 0
# This reinitializes to whatever length of ens we want to average.
if (num_av < 0) or ('ens' is not locals()):
ens = {}
FIXOFFSET = 0
n = np.int(abs(num_av))
hdr, pos = self.rd_hdr(fd)
if not isinstance(hdr, dict):
ens = -1
cfg = np.NaN
return ens, hdr, cfg, pos
cfg = self.rd_fix(fd)
fd.seek(pos, 0)
del ens
ens = {}
ens.update({'number': np.zeros((n)),
'rtc': np.zeros((7, n), dtype=np.float),
'BIT': np.zeros((n)),
'ssp': np.zeros((n)), 'depth': np.zeros((n)),
'pitch': np.zeros((n)), 'roll': np.zeros((n)),
'heading': np.zeros((n)), 'temperature': np.zeros((n)),
'salinity': np.zeros((n)), 'mpt': np.zeros((n)),
'heading_std': np.zeros((n)),
'pitch_std': np.zeros((n)), 'roll_std': np.zeros((n)),
'adc': np.zeros((8, n)),
'error_status_wd': np.zeros((n)),
'pressure': np.zeros((n)),
'pressure_std': np.zeros((n)),
'east_vel': np.zeros((cfg['n_cells'], n)),
'north_vel': np.zeros((cfg['n_cells'], n)),
'vert_vel': np.zeros((cfg['n_cells'], n)),
'error_vel': np.zeros((cfg['n_cells'], n)),
'intens': np.zeros((cfg['n_cells'], 4, n),
dtype=np.float, order='C'),
'percent': np.zeros((cfg['n_cells'], 4, n),
dtype=np.float, order='C'),
'corr': np.zeros((cfg['n_cells'], 4, n),
dtype=np.float, order='C'),
'status': np.zeros((cfg['n_cells'], 4, n),
dtype=np.float, order='C'),
'bt_range': np.zeros((4, n)),
'bt_vel': np.zeros((4, n)),
'bt_corr': np.zeros((4, n)),
'bt_ampl': np.zeros((4, n)),
'bt_perc_good': np.zeros((4, n)),
'smtime': np.zeros((n)), 'emtime': np.zeros((n)),
'slatitude': np.zeros((n)),
'slongitude': np.zeros((n)),
'elatitude': np.zeros((n)),
'elongitude': np.zeros((n)),
'nmtime': np.zeros((n)), 'flags': np.zeros((n))})
num_av = np.int(abs(num_av))
k = 0
while k < num_av:
# This is in case junk appears in the middle of a file.
num_search = 6000
id1 = np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 2+1)])
search_cnt = 0
while (search_cnt < num_search) and (((id1[0] != int('7F',16)) or (id1[1] != int('7F', 16))) or self.checkheader(fd) == 0):
search_cnt = search_cnt+1
tmp = fd.read(1)
if tmp == b'':
nextbyte = []
else:
nextbyte, = st.unpack("B", fd.read(1))
# End of file
if not nextbyte:
print(('EOF reached after {} bytes searched' +
'for next valid ensemble start').format(search_cnt))
ens = -1
return ens, hdr, cfg, pos
id1[1] = id1[0]
id1[0] = nextbyte
if search_cnt == num_search:
print(('Searched {} entries...Not a workhorse/' +
'broadband file or bad data encountered: -> {}'
).format(search_cnt, id1))
elif search_cnt > 0:
print(('Searched {} bytes to find next ' +
'valid ensemble start').format(search_cnt))
# Starting position.
startpos = fd.tell() - 2
# Read the # data types.
[hdr, nbyte] = self.rd_hdrseg(fd)
byte_offset = nbyte+2
# Read all the data types.
for n in range(1, len(hdr['dat_offsets']) + 1):
id_, = st.unpack("H", fd.read(2))
# handle all the various segments of data.
# Note that since I read the IDs as a two
# byte number in little-endian order the
# high and low bytes are exchanged compared to
# the values given in the manual.
winrivprob = 0
kol = [int('2100', 16),
int('2101', 16),
int('2102', 16),
int('2103', 16),
int('2104', 16)]
if any(x == id_ for x in kol):
winrivprob = 1
if id_ == int('0000', 16): # Fixed leader
cfg, nbyte = self.rd_fixseg(fd)
nbyte = nbyte+2
elif id_ == int('0080', 16): # Variable Leader
k = k+1
ens['number'][k-1] = st.unpack("H", fd.read(2))[0]
ens['rtc'][:, k-1] = np.array([st.unpack("B", fd.read(1)
)[0]
for ui in range(1, 7+1)])
ens['number'][k-1] = ens['number'][k-1] +\
65536 * st.unpack("B", fd.read(1))[0]
ens['BIT'][k-1] = st.unpack("H", fd.read(2))[0]
ens['ssp'][k-1] = st.unpack("H", fd.read(2))[0]
# meters
ens['depth'][k-1] = st.unpack("H", fd.read(2))[0] * .1
# degrees
ens['heading'][k-1] = st.unpack("H", fd.read(2))[0] * .01
# degrees
ens['pitch'][k-1] = st.unpack("h", fd.read(2))[0] * .01
# degrees
ens['roll'][k-1] = st.unpack("h", fd.read(2))[0] * .01
# PSU
ens['salinity'][k-1] = st.unpack("h", fd.read(2))[0]
# Deg C
ens['temperature'][k-1] = st.unpack("h", fd.read(2)
)[0] * .01
# seconds
ens['mpt'][k-1] = sum(np.array([st.unpack("B", fd.read(1)
)[0]
for ui in range(1, 3+1)]
) * np.array([60, 1, .01]))
# degrees
ens['heading_std'][k-1] = st.unpack("B", fd.read(1))[0]
# degrees
ens['pitch_std'][k-1] = st.unpack("B", fd.read(1))[0] * .1
# degrees
ens['roll_std'][k-1] = st.unpack("B", fd.read(1))[0] * .1
ens['adc'][:, k-1] = np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 8+1)])
nbyte = 2+40
if cfg['name'] == 'bb-adcp':
if cfg['prog_ver'] >= 5.55:
# 14 zeros and one byte for number WM4 bytes
fd.seek(15, 1)
cent = st.unpack("B", fd.read(1))[0]
ens['rtc'][:, k-1] =\
np.array([st.unpack("B",
fd.read(1))[0]
for ui in range(1, 7+1)])
ens['rtc'][0, k-1] = ens['rtc'][0, k-1] +\
cent * 100
nbyte = nbyte + 15 + 8
# for WH versions
elif cfg['name'] == 'wh-adcp':
ens['error_status_wd'][k-1] = st.unpack("I",
fd.read(4))[0]
nbyte = nbyte+4
if np.fix(cfg['prog_ver']) == 8\
or np.fix(cfg['prog_ver']) == 10\
or np.fix(cfg['prog_ver']) == 16\
or np.fix(cfg['prog_ver']) == 50\
or np.fix(cfg['prog_ver']) == 51\
or np.fix(cfg['prog_ver']) == 52:
if cfg['prog_ver'] >= 8.13:
fd.seek(2, 1)
ens['pressure'][k-1] = st.unpack("I",
fd.read(4))[0]
ens['pressure_std'][k-1] = st.unpack("I",
fd.read(4)
)[0]
nbyte = nbyte+10
if cfg['prog_ver'] >= 8.24:
fd.seek(1, 1)
nbyte = nbyte+1
if (cfg['prog_ver'] >= 10.01
and cfg['prog_ver'] <= 10.99)\
or cfg['prog_ver'] >= 16.05:
cent = st.unpack("B", fd.read(1))[0]
ens['rtc'][:, k-1] =\
np.array([st.unpack("B",
fd.read(1))[0]
for ui in range(1, 7+1)])
ens['rtc'][0, k-1] = ens['rtc'][0, k-1] +\
cent * 100
nbyte = nbyte+8
elif np.fix(cfg['prog_ver']) == 9:
fd.seek(2, 1)
ens['pressure'][k-1] = st.unpack("I",
fd.read(4))[0]
ens['pressure_std'][k-1] = st.unpack("I",
fd.read(4))[0]
nbyte = nbyte+10
if cfg['prog_ver'] >= 9.10:
fd.seek(1, 1)
nbyte = nbyte+1
elif cfg['name'] == 'os-adcp':
fd.seek(16, 1)
yte = nbyte + 16
if cfg['prog_ver'] > 23:
fd.seek(2, 1)
nbyte = nbyte+2
# Velocities m/s
elif id_ == int('0100', 16):
vels = np.array([st.unpack("h", fd.read(2))[0]
for ui in range(1, (4 * cfg['n_cells'])
+ 1)
]).reshape((cfg['n_cells'], 4)) * .001
ens['east_vel'][:, k-1] = vels[:, 0]
ens['north_vel'][:, k-1] = vels[:, 1]
ens['vert_vel'][:, k-1] = vels[:, 2]
ens['error_vel'][:, k-1] = vels[:, 3]
nbyte = 2+4 * cfg['n_cells'] * 2
# Correlations
elif id_ == int('0200', 16):
ens['corr'][:, :, k-1] =\
np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1,
(4 * cfg['n_cells']) + 1)
]).reshape((cfg['n_cells'], 4))
nbyte = 2+4 * cfg['n_cells']
# Echo Intensities
elif id_ == int('0300', 16):
ens['intens'][:, :, k-1] =\
np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, (4 * cfg['n_cells']) + 1)
]).reshape((cfg['n_cells'], 4))
nbyte = 2+4 * cfg['n_cells']
# Percent good
elif id_ == int('0400', 16):
ens['percent'][:, :, k-1] =\
np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, (4 * cfg['n_cells']) + 1)
]).reshape((cfg['n_cells'], 4))
nbyte = 2+4 * cfg['n_cells']
# Status
elif id_ == int('0500', 16):
if cfg['name'] is 'os-adcp':
fd.seek(00, 1)
nbyte = 2+00
else:
# Note in one case with a 4.25 firmware SC-BB,
# it seems like this block was actually two bytes short
ens['status'][:, :, k-1] =\
np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1,
(4 * cfg['n_cells']) + 1)
]).reshape((cfg['n_cells'], 4))
nbyte = 2 + 4 * cfg['n_cells']
elif id_ == int('0600', 16): # Bottom track
# In WINRIVER GPS data is tucked
# into here in odd ways, as long as GPS is enabled.
if self.SOURCE == 2:
fd.seek(2, 1)
long1 = st.unpack("H", fd.read(2))[0]
fd.seek(6, 1)
cfac = 180 / (2 ** 31)
ens['slatitude'][k-1] = st.unpack("i", fd.read(4)
)[0] * cfac
if ens['slatitude'][k-1] == 0:
ens['slatitude'][k-1] = np.NaN
else:
# Skip over a bunch of stuff
fd.seek(14, 1)
ens['bt_range'][:, k-1] =\
np.array([st.unpack("H", fd.read(2))[0]
for ui in range(1, 4+1)]) * .01
ens['bt_vel'][:, k-1] =\
np.array([st.unpack("h", fd.read(2))[0]
for ui in range(1, 4+1)])
ens['bt_corr'][:, k-1] =\
np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 4+1)])
ens['bt_ampl'][:, k-1] =\
np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 4+1)])
ens['bt_perc_good'][:, k-1] =\
np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 4+1)])
if self.SOURCE == 2:
fd.seek(2, 1)
ens['slongitude'][k-1] = (long1+65536 * st.unpack("H",
fd.read(2))[0]) * cfac
if ens['slongitude'][k-1] > 180:
ens['slongitude'][k-1] =\
ens['slongitude'][k-1] - 360
if ens['slongitude'][k-1] == 0:
ens['slongitude'][k-1] = np.NaN
fd.seek(16, 1)
qual = st.unpack("B", fd.read(1))[0]
if qual == 0:
ens['slatitude'][k-1] = np.NaN
ens['slongitude'][k-1] = np.NaN
fd.seek(71-45-21, 1)
else:
fd.seek(71-45, 1)
nbyte = 2+68
if cfg['prog_ver'] >= 5.3:
fd.seek(78-71, 1)
ens['bt_range'][:, k-1] =\
ens['bt_range'][:, k-1] +\
np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 4+1)]) * 655.36
nbyte = nbyte+11
if cfg['name'] is 'wh-adcp':
if cfg['prog_ver'] >= 16.20:
fd.seek(4, 1)
nbyte = nbyte+4
# The raw files produced by VMDAS contain a binary navigation data
# block.
elif id_ == int('2000', 16):
cfg['sourceprog'] = 'VMDAS'
if self.SOURCE != 1:
print('\n***** Apparently a VMDAS file \n\n')
self.SOURCE = 1
utim = np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 4+1)])
mtime = datetime.date.toordinal(datetime.date(utim[2] +
utim[3] *
256,
utim[1],
utim[0])
) + 366
ens['smtime'][k-1] = mtime + st.unpack("I", fd.read(4)
)[0] / 8640000
fd.seek(4, 1) # PC clock offset from UTC
cfac = 180 / 2 ^ 31
ens['slatitude'][k-1] = st.unpack("i", fd.read(4)
)[0] * cfac
ens['slongitude'][k-1] = st.unpack("i", fd.read(4)
)[0] * cfac
ens['emtime'][k-1] = mtime+st.unpack("I", fd.read(4)
)[0] / 8640000
ens['elatitude'][k-1] = st.unpack("i", fd.read(4)
)[0] * cfac
ens['elongitude'][k-1] = st.unpack("i", fd.read(4)
)[0] * cfac
fd.seek(12, 1)
ens['flags'][k-1] = st.unpack("H", fd.read(2))[0]
fd.seek(6, 1)
utim = np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, 4+1)])
mtime = datetime.date.toordinal(datetime.date(utim[0] +
utim[1] *
256,
utim[3],
utim[2])
) + 366
ens['nmtime(k)'][k-1] = mtime + st.unpack("I", fd.read(4)
)[0] / 8640000
fd.seek(16, 1)
nbyte = 2+76
# New NMEA data block from WInRiverII
elif id_ == int('2022', 16):
cfg['sourceprog'] = 'WINRIVER2'
if self.SOURCE != 2:
print('\n***** Apparently a WINRIVER file -' +
'Raw NMEA data handler not yet implemented\n\n')
self.SOURCE = 2
specID = st.unpack("H", fd.read(2))[0]
msgsiz = st.unpack("h", fd.read(2))[0]
deltaT = [str(st.unpack("s", fd.read(1)))[3]
for ui in range(1, 8+1)]
nbyte = 2+12
fd.seek(msgsiz, 1)
nbyte = nbyte + msgsiz
if specID == 100:
pass
elif specID == 101:
pass
elif specID == 102:
pass
elif specID == 103:
pass
elif id_ == int('2100', 16):
cfg['sourceprog'] = 'WINRIVER'
if self.SOURCE != 2:
print('\n***** Apparently a WINRIVER file -' +
'Raw NMEA data handler not yet implemented\n\n')
self.SOURCE = 2
str_ = [str(st.unpack("s", fd.read(1)))[3]
for ui in range(1, 38+1)]
nbyte = 2+38
elif id_ == int('2101', 16):
cfg['sourceprog'] = 'WINRIVER'
if self.SOURCE != 2:
print('\n***** Apparently a WINRIVER file -' +
'Raw NMEA data handler not yet implemented\n\n')
self.SOURCE = 2
str_ = [str(st.unpack("s", fd.read(1)))[3]
for ui in range(1, 97+1)]
nbyte = 2+97
l = str_.find('$GPGGA', 0)
if l > -1:
ens['smtime'][k-1] = (np.double(str_[l+6: l+7]) +
(np.double(str_[l+8: l+9]) +
np.double(str_[l+10: l+11]) /
60) / 60) / 24
elif id_ == int('2102', 16):
cfg['sourceprog'] = 'WINRIVER'
if self.SOURCE != 2:
print('\n***** Apparently a WINRIVER file -' +
'Raw NMEA data handler not yet implemented\n\n')
self.SOURCE = 2
str_ = [str(st.unpack("s", fd.read(1)))[3]
for ui in range(1, 45+1)]
nbyte = 2+45
elif id_ == int('2103', 16):
cfg['sourceprog'] = 'WINRIVER'
if self.SOURCE != 2:
print('\n***** Apparently a WINRIVER file -' +
'Raw NMEA data handler not yet implemented\n\n')
self.SOURCE = 2
str_ = [str(st.unpack("s", fd.read(1)))[3]
for ui in range(1, 60+1)]
nbyte = 2+60
elif id_ == int('2104', 16):
cfg['sourceprog'] = 'WINRIVER'
if self.SOURCE != 2:
print('\n***** Apparently a WINRIVER file -' +
'Raw NMEA data handler not yet implemented\n\n')
self.SOURCE = 2
str_ = [str(st.unpack("s", fd.read(1)))[3]
for ui in range(1, 38+1)]
nbyte = 2+38
# Number of good pings
elif id_ == int('0701', 16):
fd.seek(4 * cfg['n_cells'], 1)
nbyte = 2+4 * cfg['n_cells']
# Sum of squared velocities
elif id_ == int('0702', 16):
fd.seek(4 * cfg['n_cells'], 1)
nbyte = 2+4 * cfg['n_cells']
# Sum of velocities
elif id_ == int('0703', 16):
fd.seek(4 * cfg['n_cells'], 1)
nbyte = 2+4 * cfg['n_cells']
# These blocks were implemented for 5-beam systems
# Beam 5 velocity (not implemented)
elif id_ == int('0A00', 16):
fd.seek(cfg['n_cells'], 1)
nbyte = 2 + cfg['n_cells']
# Beam 5 Number of good pings (not implemented)
elif id_ == int('0301', 16):
fd.seek(cfg['n_cells'], 1)
nbyte = 2 + cfg['n_cells']
# Beam 5 Sum of squared velocities (not implemented)
elif id_ == int('0302', 16):
fd.seek(cfg['n_cells'], 1)
nbyte = 2 + cfg.n_cells
# Beam 5 Sum of velocities (not implemented)
elif id_ == int('0303', 16):
fd.seek(cfg['n_cells'], 1)
nbyte = 2 + cfg['n_cells']
# Ambient sound profile (not implemented)
elif id_ == int('020C', 16):
fd.seek(4, 1)
nbyte = 2+4
# Fixed attitude data format for OS-ADCPs (not implemented)
elif id_ == int('3000', 16):
fd.seek(32, 1)
nbyte = 2+32
else:
if hex(id_)[2:3+1] == '30':
tmp_ = hex(id_)[2:]
nflds = sum(np.base_repr((int(tmp_[2:3+1], 16)
& int('3C', 16)), 2) == '1')
dfac = sum(np.base_repr((int(tmp_[2], 16)
& int('C', 16)), 2) == '1')
fd.seek(12 * nflds * dfac, 1)
nbyte = 2+12 * nflds * dfac
else:
print(('Unrecognized ID code: {}\n').format(id_))
nbyte = 2
ens = -1
return ens, hdr, cfg, pos
byte_offset = byte_offset + nbyte
if n < len(hdr['dat_offsets']):
if hdr['dat_offsets'][n] != byte_offset:
if winrivprob == 0:
print(('{}: Adjust location by {}\n'
).format(hex(id_)[2:],
hdr['dat_offsets'][n] -
byte_offset))
fd.seek(hdr['dat_offsets'][n] - byte_offset, 1)
byte_offset = hdr['dat_offsets'][n]
else:
if hdr['nbyte'] - 2 != byte_offset:
if winrivprob == 0:
print(('{}: Adjust location by {}\n'
).format(id, hdr['nbyte'] - 2 - byte_offset
))
fd.seek(hdr['nbyte'] - 2 - byte_offset, 1)
byte_offset = hdr['nbyte'] - 2
readbytes = fd.tell() - startpos
offset = (hdr['nbyte'] + 2) - byte_offset
if (offset != 4) and (FIXOFFSET == 0):
print('\n**************************************************\n')
if fd.tell() == file_size:
print(' EOF reached unexpectedly -' +
'discarding this last ensemble\n')
ens = -1
else:
print(('Adjust location by {} ' +
'(readbytes={}, hdr.nbyte={})\n'
).format(offset, readbytes, hdr['nbyte']))
print(' NOTE - If this appears at the beginning' +
'of the read, it is\n')
print(' is a program problem, possibly fixed by a fudge')
print('-If this appears at the end of the file it means')
print('The file is corrupted and only a partial' +
'record has been read')
print('*************************************\n')
FIXOFFSET = offset-4
fd.seek(4 + FIXOFFSET, 1)
big_err = (np.array([], dtype=int), np.array([], dtype=int))
# Blank out invalid data
ens['east_vel'][(ens['east_vel'] == -32.768)] = np.NaN
ens['east_vel'][big_err[0], big_err[1]] = np.NaN
ens['north_vel'][(ens['north_vel'] == -32.768)] = np.NaN
ens['north_vel'][big_err[0], big_err[1]] = np.NaN
ens['vert_vel'][(ens['vert_vel'] == -32.768)] = np.NaN
ens['vert_vel'][big_err[0], big_err[1]] = np.NaN
ens['error_vel'][(ens['error_vel'] == -32.768)] = np.NaN
ens['error_vel'][big_err[0], big_err[1]] = np.NaN
return ens, hdr, cfg, pos, self.SOURCE
# --------------------------------------
def nmedian(self, *args, **kargs):
# Copied from median but with handling of NaN different.
nargin = len(args)
x = args[0]
x = x.astype(float)
window = args[1]
size = np.array([np.shape(x)[op]
for op in range(0, len(np.shape(x)))])
if nargin == 2:
dim_ = min(size[size != 1])
dim = min(np.where(size != 1)[0])
if not dim:
dim = 0
elif nargin == 3:
dim = args[2] - 1
add_ = (dim+1) - len(size)
if add_ <= 0:
add = []
else:
add = np.ones(1, (dim+1) - len(size))[0]
siz = np.append(size, add)
n = np.size(x, dim)
# Permute and reshape so that DIM becomes
# the row dimension of a 2-D array
perm = np.append(np.arange(dim, max(len(size) - 1, dim) + 1),
np.arange(0, dim))
x = np.transpose(x, tuple(perm)).reshape((n, np.size(x) / n))
# Sort along first dimension
x = np.sort(x, 0)
size2 = np.array([np.shape(x)[op]
for op in range(0, len(np.shape(x)))])
n1 = size2[0]
n2 = size2[1]
if n1 == 1:
y = x
else:
if n2 == 1:
# kk=sum(finite(x),1)
kk = sum(np.isfinite(x), 0)
if kk > 0:
x1 = x[np.fix((kk-1) / 2).astype(int)]
x2 = x[np.fix(kk / 2).astype(int)]
x[abs(x - (x1 + x2) / 2) > window] = np.NaN
x = np.sort(x, 1)
kk = sum(np.isfinite(x), 1)
x[np.isnan(x)] = 0
y = np.NaN
if kk > 0:
y = sum(x) / kk
else:
kk = sum(np.isfinite(x), 0)
ll = kk < n1-2
kk[ll] = 0
x[:, ll] = x[:, ll] * np.NaN
xx = np.transpose(x).reshape((n1 * n2, 1))
x1 = xx[np.fix((kk-1) / 2).astype(int) +
1 + np.array(np.arange(0, n2)) * n1-1]
x2 = xx[np.fix(kk / 2).astype(int) +
1 + np.array(np.arange(0, n2)) * n1-1]
x[(abs(x - np.dot(np.ones((n1, 1)),
np.transpose((x1 + x2) / 2.0))) > window)
] = np.NaN
x = np.sort(x, 0)
kk = sum(np.isfinite(x), 0)
x[np.isnan(x)] = 0
y = np.NaN + np.ones((1, n2))[0]
if any(kk):
tmp_ = x.reshape((n2, n1))
y[np.where(kk > 0)[0]] = np.sum(x[:, np.where(kk > 0)[0]],
axis=0) / kk[kk > 0]
siz[dim] = 1
inverseorder = perm * 0
inverseorder[perm-1] = np.arange(0, len(perm+1))
if len(np.shape(y)) <= 1:
return y
else:
y = y.reshape(tuple(siz[perm]))
y = np.transpose(y, tuple(inverseorder))
return y
# --------------------------------------
def nmean(self, *args, **kargs):
nargin = len(args)
x = args[0]
x = x.astype(float)
kk = np.isfinite(x)
x[kk is False] = 0
size = np.array([np.shape(x)[op]
for op in range(0, len(np.shape(x)))])
if nargin == 1:
# Determine which dimension SUM will use
dim_ = min(size[size != 1])
dim = min(np.where(size != 1)[0])
if not dim:
dim = 0
elif nargin == 2:
dim = args[1] - 1
if dim+1 > len(size):
y = x
else:
ndat = np.sum(kk, dim)
indat = ndat == 0
ndat[indat] = 1
y = np.sum(x, dim) / ndat
y[indat] = np.NaN
return y
def read_ADCP(self, filename, pathname, numav=1, nens=-1, bas=None,
despike=None):
varargin = [filename, numav, nens, 'baseyear', bas, 'despike', despike]
self.SOURCE = 0
name = varargin[0]
varargin.remove(varargin[0])
if ((varargin[0] is None)or (varargin[1] is None)) and\
((varargin[0] is not None) or (varargin[1] is not None)) :
raise ValueError("both first two args must be " +
"simultaneously None or not None")
elif (varargin[0] and varargin[1]) is None:
[varargin.remove(varargin[x]) for x in [1, 0]]
froom = [i for i, v in enumerate(varargin) if v is None]
froomp1 = [y-1 for y in froom]
froom = froom + froomp1
froom.sort(reverse=True)
[varargin.remove(varargin[x]) for x in froom]
print(varargin)
century = 2000
vels = 'no' # Default to simple averaging
lv = len(varargin)
if lv >= 1 and isinstance(varargin[0], str) is False:
num_av = varargin[0]
varargin.remove(varargin[0])
lv = lv-1
if lv >= 1 and isinstance(varargin[0], str) is False:
nens = varargin[0]
varargin.remove(varargin[0])
lv = lv-1
# Read optional args
while len(varargin) > 0:
if varargin[0][0:3] == 'bas':
century = varargin[1]
elif varargin[0][0:3] == 'des':
if isinstance(varargin[1], str):
if varargin[1] == 'no':
vels = 'no'
else:
vels = [.3, .3, .3]
else:
vels = varargin[1]
else:
print(('Unknown command line option -> {}'
).format(varargin[0]))
del varargin[0: 2]
# Check file information first
tutu = os.stat(name)
naminfo = {'name': name, 'date': ti.gmtime(tutu[-2]), 'bytes': tutu[6]}
if not naminfo:
print(('ERROR******* Can''t find file {}\n').format(name))
print(('\nOpening file {}\n\n').format(name))
fd = open(name, 'rb')
fd.seek(0, 2)
file_size = fd.tell()
fd.seek(0)
# Read first ensemble to initialize parameters
# Initialize and read first two records
[ens, hdr, cfg, pos, self.SOURCE] = self.rd_buffer(fd, -2)
# from datetime import datetime, timedelta
# def datetime2matlabdn(self,dt):
# ord = dt.toordinal()
# mdn = dt + datetime.timedelta(days = 366)
# frac = (dt-datetime.datetime(dt.year,
# dt.month,dt.day,0,0,0)).seconds / (24.0 * 60.0 * 60.0)
# return mdn.toordinal() + frac
if not isinstance(ens, dict) and ens == -1:
print('No Valid data found')
adcp = {}
return adcp
fd.seek(pos, 0)
if (cfg['prog_ver'] < 16.05
and cfg['prog_ver'] > 5.999) \
or cfg['prog_ver'] < 5.55:
print(('***** Assuming that the century begins' +
'year {} (info not in this firmware ' +
'version)').format(century))
else:
century = 0 # century included in clock.
ens['rtc'] = ens['rtc'].astype(int)
dats_ = [datetime.datetime(century + ens['rtc'][0, ui],
ens['rtc'][1, ui],
ens['rtc'][2, ui],
ens['rtc'][3, ui],
ens['rtc'][4, ui],
ens['rtc'][5, ui],
ens['rtc'][6, ui])
for ui in range(0, int(np.shape(ens['rtc'])[1]))]
dats = [ti.mktime(dats_[i].timetuple())
for i in range(0, int(np.shape(ens['rtc'])[1]))]
t_int = np.diff(dats)
print(('Record begins at {}\n'
).format(ti.strftime('%Y-%m-%d %H:%M:%S',
ti.localtime(dats[0]))))
print(('Ping interval appears to be {}\n'
).format(ti.strftime('%H:%M:%S', ti.gmtime(t_int))))
extrabytes = 0
nensinfile = np.fix(naminfo['bytes'] / (hdr['nbyte'] + 2 + extrabytes))
print(('\nEstimating {} ensembles in this file\n').format(nensinfile))
if len([nens]) == 1:
if nens == -1:
nens = nensinfile
print(('Reading {} ensembles, reducing by a factor of {}\n'
).format(nens, num_av))
else:
print(('Reading ensembles {}-{}, reducing by a factor of {}\n'
).format(nens, num_av))
fd.seek((hdr['nbyte'] + 2 + extrabytes) * (nens[0] - 1), 'cof')
nens = int(np.diff(nens) + 1)
nens = int(nens)
# Number of records after averaging.
n = int(np.fix(nens / num_av))
print(('Final result {} values\n').format(n))
if num_av > 1:
if isinstance(vels, str):
print('\n Simple mean used for ensemble averaging\n')
else:
print(('\n Averaging after outlier rejection with' +
'parameters [{}]\n').format(vels))
# Structure to hold all ADCP data
# Note that I am not storing all the data
# contained in the raw binary file, merely
# things I think are useful.
if cfg['sourceprog'] == 'WINRIVER':
adcp = {'config': cfg, 'mtime': np.zeros((n)),
'rmtime': np.zeros((nens)), 'number': np.zeros((n)),
'pitch': np.zeros((n)), 'roll': np.zeros((n)),
'heading': np.zeros((n)), 'pitch_std': np.zeros((n)),
'roll_std': np.zeros((n)), 'heading_std': np.zeros((n)),
'depth': np.zeros((n)), 'temperature': np.zeros((n)),
'salinity': np.zeros((n)), 'pressure': np.zeros((n)),
'pressure_std': np.zeros((n)),
'east_vel': np.zeros((cfg['n_cells'], n)),
'north_vel': np.zeros((cfg['n_cells'], n)),
'vert_vel': np.zeros((cfg['n_cells'], n)),
'error_vel': np.zeros((cfg['n_cells'], n)),
'corr': np.zeros((cfg['n_cells'], 4, n),
dtype=np.float, order='C'),
'status': np.zeros((cfg['n_cells'], 4, n),
dtype=np.float, order='C'),
'intens': np.zeros((cfg['n_cells'], 4, n),
dtype=np.float, order='C'),
'bt_range': np.zeros((4, n)), 'bt_vel': np.zeros((4, n)),
'bt_corr': np.zeros((4, n)), 'bt_ampl': np.zeros((4, n)),
'bt_perc_good': np.zeros((4, n)),
'nav_mtime': np.zeros((n)),
'nav_longitude': np.zeros((n)),
'nav_latitude': np.zeros((n)),
'corr_raw': np.zeros((cfg['n_cells'], 4, nens),
dtype=np.float, order='C'),
'status_raw': np.zeros((cfg['n_cells'], 4, nens),
dtype=np.float, order='C'),
'intens_raw': np.zeros((cfg['n_cells'], 4, nens),
dtype=np.float, order='C')}
elif cfg['sourceprog'] == 'VMDAS':
adcp = {'config': cfg, 'mtime': np.zeros((n)),
'rmtime': np.zeros((nens)), 'number': np.zeros((n)),
'pitch': np.zeros((n)), 'roll': np.zeros((n)),
'heading': np.zeros((n)), 'pitch_std': np.zeros((n)),
'roll_std': np.zeros((n)), 'heading_std': np.zeros((n)),
'depth': np.zeros((n)), 'temperature': np.zeros((n)),
'salinity': np.zeros((n)), 'pressure': np.zeros((n)),
'pressure_std': np.zeros((n)),
'east_vel': np.zeros((cfg['n_cells'], n)),
'north_vel': np.zeros((cfg['n_cells'], n)),
'vert_vel': np.zeros((cfg['n_cells'], n)),
'error_vel': np.zeros((cfg['n_cells'], n)),
'corr': np.zeros((cfg['n_cells'], 4, n),
dtype=np.float, order='C'),
'status': np.zeros((cfg['n_cells'], 4, n),
dtype=np.float, order='C'),
'intens': np.zeros((cfg['n_cells'], 4, n),
dtype=np.float, order='C'),
'bt_range': np.zeros((4, n)), 'bt_vel': np.zeros((4, n)),
'bt_corr': np.zeros((4, n)), 'bt_ampl': np.zeros((4, n)),
'bt_perc_good': np.zeros((4, n)),
'nav_smtime': np.zeros((n)), 'nav_emtime': np.zeros((n)),
'nav_slongitude': np.zeros((n)),
'nav_elongitude': np.zeros((n)),
'nav_slatitude': np.zeros((n)),
'nav_elatitude': np.zeros((n)),
'nav_mtime': np.zeros((n)),
'corr_raw': np.zeros((cfg['n_cells'], 4, nens),
dtype=np.float, order='C'),
'status_raw': np.zeros((cfg['n_cells'], 4, nens),
dtype=np.float, order='C'),
'intens_raw': np.zeros((cfg['n_cells'], 4, nens),
dtype=np.float, order='C')}
else:
adcp = {'config': cfg, 'mtime': np.zeros((n)),
'rmtime': np.zeros((int(nens))),
'number': np.zeros((n)),
'pitch': np.zeros((n)), 'roll': np.zeros((n)),
'heading': np.zeros((n)), 'pitch_std': np.zeros((n)),
'roll_std': np.zeros((n)), 'heading_std': np.zeros((n)),
'depth': np.zeros((n)), 'temperature': np.zeros((n)),
'salinity': np.zeros((n)), 'pressure': np.zeros((n)),
'pressure_std': np.zeros((n)),
'east_vel': np.zeros((cfg['n_cells'], n)),
'north_vel': np.zeros((cfg['n_cells'], n)),
'vert_vel': np.zeros((cfg['n_cells'], n)),
'error_vel': np.zeros((cfg['n_cells'], n)),
'corr': np.zeros((cfg['n_cells'], 4, n),
dtype=np.float, order='C'),
'status': np.zeros((cfg['n_cells'], 4, n),
dtype=np.float, order='C'),
'intens': np.zeros((cfg['n_cells'], 4, n),
dtype=np.float, order='C'),
'bt_range': np.zeros((4, n)), 'bt_vel': np.zeros((4, n)),
'bt_corr': np.zeros((4, n)), 'bt_ampl': np.zeros((4, n)),
'bt_perc_good': np.zeros((4, n)),
'corr_raw': np.zeros((cfg['n_cells'], 4, nens),
dtype=np.float, order='C'),
'status_raw': np.zeros((cfg['n_cells'], 4, nens),
dtype=np.float, order='C'),
'intens_raw': np.zeros((cfg['n_cells'], 4, nens),
dtype=np.float, order='C')}
# Calibration factors for backscatter data
del ens
# Loop for all records
for k in range(1, int(n)):
# Gives display so you know something is going on...
if k % 50 == 0:
print(('\n{}').format(k * num_av))
# Read an ensemble
ens = self.rd_buffer(fd, num_av)[0]
if not isinstance(ens, dict): # If aborting...
print(('Only {} records found..' +
'suggest re-running RDRADCP ' +
'using this parameter\n').format((k-1) * num_av))
print(('(If this message preceded by a POSSIBLE PROGRAM ' +
'PROBLEM message, re-run using {})\n'
).format((k-1) * num_av-1))
break
ens['rtc'] = ens['rtc'].astype(int)
dats_ = [datetime.datetime(century + ens['rtc'][0, ui],
ens['rtc'][1, ui],
ens['rtc'][2, ui],
ens['rtc'][3, ui],
ens['rtc'][4, ui],
ens['rtc'][5, ui],
ens['rtc'][6, ui])
for ui in range(0, int(np.shape(ens['rtc'])[1]))]
dats = [ti.mktime(dats_[i].timetuple())
for i in range(0, int(np.shape(ens['rtc'])[1]))]
adcp['mtime'][k-1] = np.median(dats)
adcp['rmtime'][(k-1) * num_av: k * num_av] = dats
adcp['number'][k-1] = ens['number'][0]
adcp['heading'][k-1] = np.mean(ens['heading'])
adcp['pitch'][k-1] = np.mean(ens['pitch'])
adcp['roll'][k-1] = np.mean(ens['roll'])
adcp['heading_std'][k-1] = np.mean(ens['heading_std'])
adcp['pitch_std'][k-1] = np.mean(ens['pitch_std'])
adcp['roll_std'][k-1] = np.mean(ens['roll_std'])
adcp['depth'][k-1] = np.mean(ens['depth'])
adcp['temperature'][k-1] = np.mean(ens['temperature'])
adcp['salinity'][k-1] = np.mean(ens['salinity'])
adcp['pressure'][k-1] = np.mean(ens['pressure'])
adcp['pressure_std'][k-1] = np.mean(ens['pressure_std'])
if isinstance(vels, str):
adcp['east_vel'][:, k-1] = np.nanmean(ens['east_vel'], axis=1)
adcp['north_vel'][:, k-1] = np.nanmean(ens['north_vel'],
axis=1)
adcp['vert_vel'][:, k-1] = np.nanmean(ens['vert_vel'], axis=1)
adcp['error_vel'][:, k-1] = np.nanmean(ens['error_vel'],
axis=1)
else:
adcp['east_vel'][:, k-1] = self.nmedian(ens['east_vel'],
vels[0], 2)
adcp['north_vel'][:, k-1] = self.nmedian(ens['north_vel'],
vels[0], 2)
adcp['vert_vel'][:, k-1] = self.nmedian(ens['vert_vel'],
vels[1], 2)
adcp['error_vel'][:, k-1] = self.nmedian(ens['error_vel'],
vels[2], 2)
adcp['corr'][:, :, k-1] = np.nanmean(ens['corr'], axis=2)
adcp['status'][:, :, k-1] = np.nanmean(ens['status'], axis=2)
adcp['intens'][:, :, k-1] = np.nanmean(ens['intens'], axis=2)
adcp['corr_raw'][:, :, (k-1) * num_av: (k) * num_av] = ens['corr']
adcp['status_raw'][:, :, (k-1) * num_av: (k) * num_av] =\
ens['status']
adcp['intens_raw'][:, :, (k-1) * num_av: (k) * num_av] =\
ens['intens']
if 'perc_good' not in adcp:
adcp.update({'perc_good': np.zeros((cfg['n_cells'], 4, n),
dtype=np.float, order='C')})
adcp['perc_good'][:, :, k-1] = np.nanmean(ens['percent'], axis=2)
adcp['bt_range'][:, k-1] = np.nanmean(ens['bt_range'], axis=1)
adcp['bt_vel'][:, k-1] = np.nanmean(ens['bt_vel'], axis=1)
adcp['bt_corr'][:, k-1] = np.nanmean(ens['bt_corr'], axis=1)
adcp['bt_ampl'][:, k-1] = np.nanmean(ens['bt_ampl'], axis=1)
adcp['bt_perc_good'][:, k-1] = np.nanmean(ens['bt_perc_good'],
axis=1)
intens_cal = np.dot(np.transpose([np.log10(adcp['config']['range'])
]), [np.array([1, 1, 1, 1]) * 20
]) * 0
adcp['TVG'] = intens_cal
if cfg['sourceprog'] == 'WINRIVER':
adcp['nav_mtime'][k-1] = np.nanmean(ens['smtime'], axis=0)
adcp['nav_longitude'][k-1] = np.nanmean(ens['slongitude'],
axis=0)
adcp['nav_latitude'][k-1] = np.nanmean(ens['slatitude'],
axis=0)
elif cfg['sourceprog'] == 'VMDAS':
adcp['nav_smtime'][k-1] = ens['smtime[0]']
adcp['nav_emtime'][k-1] = ens['emtime[0]']
adcp['nav_slatitude'][k-1] = ens['slatitude[0]']
adcp['nav_elatitude'][k-1] = ens['elatitude[0]']
adcp['nav_slongitude'][k-1] = ens['slongitude[0]']
adcp['nav_elongitude'][k-1] = ens['elongitude[0]']
adcp['nav_mtintenime'][k-1] = np.nanmean(ens['nmtime'], axis=0)
print('\n')
print(('Read to byte {} in a file of ' +
'size {} bytes\n').format(fd.tell(), naminfo['bytes']))
if fd.tell() + hdr['nbyte'] < naminfo['bytes']:
print(('-->There may be another {} ensembles unread\n'
).format(np.fix((naminfo['bytes'] - fd.tell()
) / (hdr['nbyte'] + 2))))
fd.close()
return Parameters(adcp)
[docs]class ReadAdcpP(object):
""" Reads -p RDI fileformat meant for en route acquisition.
Parameters
----------
obj : Adcp object :class:`hydrac.instruments.ascp.Adcp`
"""
def __init__(self, obj):
self.obj = obj
if self.obj.tag_adcpType[self.obj.Ncurrfilepath] != 'p':
raise WrongFileTypeError('The file extension' +
'does not contain navigation data...')
adcp, metadcp =\
self.obj.param_shape_ADCP(self.read_ADCP(self.obj.currfilepath,
'', 1, -1, 'none', 'yes'))
self.obj.param.update({"P" + str(self.obj.Ncurrfilepath): adcp})
self.obj.meta.update({"P" + str(self.obj.Ncurrfilepath): metadcp})
# --------------------------------------
def read_hdr(self, fd, adcp_type, n_cells):
# Reads a Header
hdr = {}
hdrid, = st.unpack("H", fd.read(2)) # Head ID
if hdrid != int('7e7f', 16):
# This corrects a case where the .000 file seemed to have
# one byte too many (prog_ver 2.72)
if fd.tell() > 3: # Not the very first segment
backnum = -10
twin = 20
fd.seek(backnum, 1)
bb = np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, twin+1)])
ff = (np.where(bb[0: twin-2] == 127) ==
np.where(bb[1: twin-1] == 126))
if ff is True:
ff = np.where(bb[0: twin-2] == 127)[0][0]
if ff >= 0:
fd.seek(- twin + (ff-1) + 2, 1)
print(('Warning - segment length wrong,' +
'correction {}\n').format(backnum + ff-1+2))
else:
print(('File ID is {} not 7E7F - data corrupted or ' +
'not a P-file?').format(hex(hdrid)))
else:
print(('File ID is {} not 7E7F - data corrupted or ' +
'not a P-file?').format(hex(hdrid)))
nbyte, = st.unpack("H", fd.read(2))
hdr.update({'nbyte': nbyte-1})
n_deps, = st.unpack("B", fd.read(1))
hdr.update({'n_deps': n_deps})
if adcp_type == 'narrowband':
hdr.update({'n_deps': n_cells})
ndat, = st.unpack("B", fd.read(1))
if ndat == 0:
dat_offsets = []
else:
dat_offsets = [st.unpack("H", fd.read(2))[0]
for ui in range(1, ndat+1)]
hdr.update({'dat_offsets': dat_offsets})
hdr['nbyte'] = hdr['nbyte'] - 6 - ndat * 2
return hdr
# -------------------------------------
def read_cfgseg(self, fd):
# Read config data
hdr = self.read_hdr(fd, 0, 0)
cfgid, = st.unpack("H", fd.read(2))
if cfgid != int('000a', 16):
print(('Cnfig ID {} is incorrect - ' +
'data corrupted or not a P-file?').format(hex(cfgid)))
cfg = self.read_cfg(hdr, fd)
fd.seek(2, 1)
return cfg
# -------------------------------------
def getopt(self, val, *varargin):
# Returns one of a list (0=first in varargin, etc.)
if val > len(varargin):
opt = 'unknown'
else:
opt = varargin[val]
return opt
# -------------------------------------
def read_cfg(self, hdr, fd):
# Reads the configuration data
# global INFO
cfg = {}
fptr = fd.tell() - 2 # Pointer to beginning of segment
# typeNames = {
# 'int8' :'b',
# 'uint8' :'B',
# 'int16' :'h',
# 'uint16' :'H',
# 'int32' :'i',
# 'uint32' :'I',
# 'int64' :'q',
# 'uint64' :'Q',
# 'float' :'f',
# 'double' :'d',
# 'char' :'s'}
tmp_ = fd.read(1)
tmp, = st.unpack("B", tmp_)
cfg.update({'adcp_type': self.getopt(tmp, 'narrowband', 'broadband')})
if cfg['adcp_type'] is 'broadband': # BB P-file format
cfg.update({'prog_ver': st.unpack("B", fd.read(1))[0] +
st.unpack("B", fd.read(1))[0] / 100})
cfg.update({'firm_ver': st.unpack("H", fd.read(2))[0] / 100})
cfg.update({'n_beams': st.unpack("B", fd.read(1))[0]})
cfg.update({'beam_angle': st.unpack("H", fd.read(2))[0]})
cfg.update({'beam_freq': st.unpack("H", fd.read(2))[0]})
cfg.update({'prof_mode': st.unpack("B", fd.read(1))[0]})
cfg.update({'coord_sys': self.getopt(st.unpack("H", fd.read(2))[0],
'beam', 'earth', 'ship',
'instrument')})
cfg.update({'orientation': self.getopt(st.unpack("H",
fd.read(2))[0],
'up', 'down')})
cfg.update({'beam_pattern': self.getopt(st.unpack("H",
fd.read(2))[0],
'convex', 'concave')})
cfg.update({'n_cells': st.unpack("B", fd.read(1))[0]})
# seconds
cfg.update({'time_between_ping_groups':
st.unpack("I", fd.read(4))[0] * .01})
cfg.update({'pings_per_ensemble': st.unpack("H", fd.read(2))[0]})
# meters
cfg.update({'cell_size':
st.unpack("H", fd.read(2))[0] * .01})
# meters
cfg.update({'blank': st.unpack("H", fd.read(2))[0]*.01})
# meters
cfg.update({'depth': st.unpack("I", fd.read(4))[0]*.01})
cfg.update({'avg_method': self.getopt(st.unpack("B", fd.read(1)
)[0], 'time',
'space')})
cfg.update({'avg_interval':
st.unpack("I", fd.read(4))[0]*.01}) # sec. or meters
cfg.update({'magnetic_var':
st.unpack("i", fd.read(4))[0]*.001}) # degrees
cfg.update({'compass_offset':
st.unpack("i", fd.read(4))[0]*.001}) # degrees
cfg.update({'xducer_misalign':
st.unpack("i", fd.read(4))[0]*.001}) # degrees
cfg.update({'intens_scale':
st.unpack("H", fd.read(2))[0]*.001}) # db/count
cfg.update({'absorption':
st.unpack("H", fd.read(2))[0]*.001}) # db/m
cfg.update({'salinity':
st.unpack("H", fd.read(2))[0]*.001}) # ppt
cfg.update({'ssp':
st.unpack("i", fd.read(4))[0]*.01}) # m/s
cfg.update({'ssp_use':
self.getopt(st.unpack("B", fd.read(1))[0],
'yes', 'no')})
cfg.update({'use_pitchroll':
self.getopt(st.unpack("B", fd.read(1))[0],
'yes', 'no')})
fd.seek(20, 1) # Bunch of stuff for discharge calcs.
cfg.update({'bin1_dist':
st.unpack("H", fd.read(2))[0]*.01}) # meters
cfg.update({'xmit_pulse':
st.unpack("H", fd.read(2))[0]*.01}) # meters
fd.seek(4, 1)
nchar, = st.unpack("B", fd.read(1))
cfg.update({'deploy_name':
[str(st.unpack("s", fd.read(1)))[3]
for ui in range(1, nchar+1)]})
fd.seek(2, 1)
s1, = st.unpack("B", fd.read(1))
tmp2_ = [str(st.unpack("s", fd.read(1)))[3]
for ui in range(1, s1+1)]
cfg.update({'cmd_set1': [str(st.unpack("s", fd.read(1)))[3]
for ui in range(1, s1+1)]})
s2, = st.unpack("B", fd.read(1))
cfg.update({'cmd_set2': [str(st.unpack("s", fd.read(1)))[3]
for ui in range(1, s2+1)]})
cfg.update({'comment': ' '*80})
print('CFG data')
offset = hdr['nbyte'] - (96 + nchar + s1 + s2)
if offset != 0:
print((' Header: excess bytes {}').format(offset))
fd.seek(offset, 1) # get to end of segment
# #########ICI############
elif cfg['adcp_type'] is 'narrowband': # NB p-format
cfg.update({'prog_ver':
st.unpack("B", fd.read(1))[0] +
st.unpack("B", fd.read(1))[0]/100})
cfg.update({'serialnum': st.unpack("H", fd.read(2))[0]})
cfg.update({'firm_ver': st.unpack("H", fd.read(2))[0]/100})
cfg.update({'n_beams': st.unpack("B", fd.read(1))[0]})
cfg.update({'beam_angle': st.unpack("B", fd.read(1))[0]})
cfg.update({'beam_freq': st.unpack("H", fd.read(2))[0]})
cfg.update({'range_switch': st.unpack("B", fd.read(1))[0]})
cfg.update({'coord_sys':
self.getopt(st.unpack("B", fd.read(1))[0],
'beam', 'earth')})
cfg.update({'orientation':
self.getopt(st.unpack("B", fd.read(1))[0],
'up', 'down')})
cfg.update({'beam_pattern':
self.getopt(st.unpack("B", fd.read(1))[0],
'convex', 'concave')})
cfg.update({'n_cells': st.unpack("B", fd.read(1))[0]})
cfg.update({'time_between_ping_groups':
st.unpack("I", fd.read(4))[0]*.01}) # seconds
cfg.update({'pings_per_ensemble': st.unpack("H", fd.read(2))[0]})
cfg.update({'cell_size': st.unpack("B", fd.read(1))[0]}) # meters
cfg.update({'xmit_pulse':
st.unpack("H", fd.read(2))[0]}) # TR pulse length (m)
cfg.update({'blank': st.unpack("H", fd.read(2))[0]}) # meters
cfg.update({'depth': st.unpack("I", fd.read(4))[0]}) # meters
cfg.update({'avg_method':
self.getopt(st.unpack("B", fd.read(1))[0],
'time', 'space')})
cfg.update({'avg_interval':
st.unpack("I", fd.read(4))[0]*.01}) # sec; or meters
cfg.update({'magnetic_var':
st.unpack("i", fd.read(4))[0]*.001}) # degrees
cfg.update({'compass_offset':
st.unpack("i", fd.read(4))[0]*.001}) # degrees
cfg.update({'xducer_misalign':
st.unpack("i", fd.read(4))[0]*.001}) # degrees
cfg.update({'intens_scale':
st.unpack("H", fd.read(2))[0]*.001}) # db/count
cfg.update({'absorption':
st.unpack("H", fd.read(2))[0]*.001}) # db/m
cfg.update({'salinity':
st.unpack("H", fd.read(2))[0]*.001}) # ppt
cfg.update({'ssp':
st.unpack("I", fd.read(4))[0]*.001}) # m/s
cfg.update({'ssp_use':
self.getopt(st.unpack("B", fd.read(1))[0],
'yes', 'no')})
cfg.update({'use_pitchroll':
self.getopt(st.unpack("B", fd.read(1))[0],
'yes', 'no')})
fd.seek(20, 1)
fd.seek(8, 1)
cfg.update({'bin1_dist':
cfg['depth'] + cfg['blank'] + cfg['cell_size'] / 2})
nchar, = st.unpack("B", fd.read(1))
cfg.update({'deploy_name': [str(st.unpack("s", fd.read(1)))[3]
for ui in range(1, nchar+1)]})
fd.seek(2, 1)
s1, = st.unpack("B", fd.read(1))
cfg.update({'cmd_set1': [str(st.unpack("s", fd.read(1)))[3]
for ui in range(1, s1+1)]})
cfg.update({'comment': ' '*80})
print('CFG data')
offset = hdr['nbyte'] - (94 + nchar + s1)
if offset != 0:
print(('Header: excess bytes {}').format(offset))
fd.seek(offset, 1) # get to end of segment
else:
print('Unknown adcp type')
# Calibration factors for backscatter data
cfg['range'] = cfg['bin1_dist'] +\
np.arange(0, cfg['n_cells'], 1, dtype=np.float) * cfg['cell_size']
if self.iNFO['create'] == 1: # Set flag to store next time
self.iNFO['create'] = 2
self.iNFO['fileptr'] = np.append(self.iNFO['fileptr'],
fptr-6-len(hdr['dat_offsets'])*2)
self.iNFO['numrec'] = np.append(self.iNFO['numrec'], 0)
return cfg
# -----------------------------
def read_ensemble(self, fd, num_av, cfg, ens=None):
# To save it being re-initialized every time.
# If num_av<0 we are reading only 1 element and initializing
if num_av < 0:
ens = {}
n = np.abs(num_av)
pos = fd.tell()
hdr = self.read_hdr(fd, cfg['adcp_type'], cfg['n_cells'])
fd.seek(pos, 0)
ens.update({'day': np.zeros((n), dtype=np.float),
'time': np.zeros((n)),
'number': np.zeros((n)),
'pitch': np.zeros((n)),
'roll': np.zeros((n)),
'heading': np.zeros((n)),
'temperature': np.zeros((n)),
'bt_east_vel': np.zeros((n)),
'bt_north_vel': np.zeros((n)),
'bt_vert_vel': np.zeros((n)),
'bt_error_vel': np.zeros((n)),
'bt_range': np.zeros((4, n)),
'bt_x_disp': np.zeros((n)),
'bt_y_disp': np.zeros((n)),
'elapsed_time': np.zeros((n)),
'bt_path_len': np.zeros((n)),
'num_averaged': np.zeros((n)),
'east_vel': np.zeros((hdr['n_deps'], n)),
'north_vel': np.zeros((hdr['n_deps'], n)),
'vert_vel': np.zeros((hdr['n_deps'], n)),
'error_vel': np.zeros((hdr['n_deps'], n)),
'intens': np.zeros((hdr['n_deps'], 4, n),
dtype=np.float, order='C'),
'percent': np.zeros((hdr['n_deps'], n)),
'nav_east_vel': np.zeros((n)),
'nav_north_vel': np.zeros((n)),
'nav_x_disp': np.zeros((n)),
'nav_y_disp': np.zeros((n)),
'nav_path_len': np.zeros((n)),
'latitude': np.zeros((n)),
'longitude': np.zeros((n)),
'filterwidth': np.zeros((n))})
num_av = 1
ens['day'][0] = 0
if cfg['adcp_type'] is 'narrowband':
bad_val = 32767
else:
bad_val = -32768
bad_vald1000 = bad_val * 0.001
k = 0
while k < num_av:
hdr = self.read_hdr(fd, cfg['adcp_type'], cfg['n_cells'])
k = k+1
for n in range(0, len(hdr['dat_offsets'])):
id_, = st.unpack("H", fd.read(2))
if id_ == int('000B', 16):
# case '000B'
if cfg['adcp_type'] is 'broadband': # BB format
# day from beginning of century
ens['day'][k-1] = st.unpack("I", fd.read(4))[0]
# secs from beginning of day
ens['time'][k-1] = st.unpack("I", fd.read(4))[0]*.01
ens['number'][k-1] = st.unpack("H", fd.read(2))[0]
# in degrees
ens['pitch'][k-1] = st.unpack("h", fd.read(2))[0]*.01
# degrees
ens['roll'][k-1] = st.unpack("h", fd.read(2))[0]*.01
# degrees
ens['heading'][k-1] = st.unpack("H", fd.read(2))[0]*.01
# degrees celsius
ens['temperature'][k-1] = st.unpack("h", fd.read(2)
)[0]*.01
# next 4 in m/s
ens['bt_east_vel'][k-1] = st.unpack("h", fd.read(2)
)[0]*.001
ens['bt_north_vel'][k-1] = st.unpack("h", fd.read(2)
)[0]*.001
ens['bt_vert_vel'][k-1] = st.unpack("h", fd.read(2)
)[0]*.001
ens['bt_error_vel'][k-1] = st.unpack("h", fd.read(2)
)[0]*.001
fd.seek(8, 1)
# m/s
ens['bt_range'][:, k-1] =\
np.array([st.unpack("H", fd.read(2))[0]
for ui in range(1, 4+1)]) * .01
# next 2 in m
ens['bt_x_disp'][k-1] = st.unpack("i", fd.read(4)
)[0]*.01
ens['bt_y_disp'][k-1] = st.unpack("i", fd.read(4)
)[0]*.01
# sec
ens['elapsed_time'][k-1] = st.unpack("I", fd.read(4)
)[0]*.01
# m
ens['bt_path_len'][k-1] = st.unpack("I", fd.read(4)
)[0]*.01
fd.seek(14, 1)
ens['num_averaged'][k-1] = st.unpack("H", fd.read(2)
)[0]
else: # see previous lines for units if not indicated
# seconds since Jan1
ens['time'][k-1] = st.unpack("I", fd.read(4))[0]
ens['number'][k-1] = st.unpack("H", fd.read(2))[0]
ens['pitch'][k-1] = st.unpack("h", fd.read(2))[0]*.01
ens['roll'][k-1] = st.unpack("h", fd.read(2))[0]*.01
ens['heading'][k-1] = st.unpack("H", fd.read(2))[0]*.01
ens['temperature'][k-1] = st.unpack("h", fd.read(2)
)[0]*.01
ens['bt_east_vel'][k-1] = st.unpack("h", fd.read(2)
)[0]*.001
ens['bt_north_vel'][k-1] = st.unpack("h", fd.read(2)
)[0]*.001
ens['bt_vert_vel'][k-1] = st.unpack("h", fd.read(2)
)[0]*.001
ens['bt_error_vel'][k-1] = st.unpack("h", fd.read(2)
)[0]*.001
ens['bt_range'][:, k-1] =\
np.array([st.unpack("H", fd.read(2))[0]
for ui in range(1, 4+1)]) * .01
ens['bt_x_disp'][k-1] = st.unpack("i", fd.read(4)
)[0]*.01
ens['bt_y_disp'][k-1] = st.unpack("i", fd.read(4)
)[0]*.01
ens['elapsed_time'][k-1] = st.unpack("I", fd.read(4)
)[0]*.01
ens['bt_path_len'][k-1] = st.unpack("I", fd.read(4)
)[0]*.01
fd.seek(14, 1)
ens['num_averaged'][k-1] = st.unpack("H", fd.read(2)
)[0]
if self.iNFO['create'] > 0:
self.iNFO['numrec'][-1] = self.iNFO['numrec'][-1] + 1
if self.iNFO['create'] == 2: # Store the time
if cfg['adcp_type'] is 'broadband':
self.iNFO['mtime'] =\
np.append(self.iNFO['mtime'],
datetime.date.toordinal(datetime.date(progyear, 1, 1))+366+ens['day'][k-1]+ens['time'][k-1]/86400)
self.iNFO['create'] = 1
else:
self.iNFO['mtime'] =\
np.append(self.iNFO['mtime'],
datetime.date.toordinal(datetime.date(progyear, 1, 1))+366+ens['time'][k-1]/86400)
self.iNFO['create'] = 1
# Velocities
elif id_ == int('0001', 16):
# m/s
vels = np.array([[st.unpack("h", fd.read(2))[0]
for ui in range(1, 4+1)]
for iu in range(1, hdr['n_deps']+1)])*.001
ens['east_vel'][:, k-1] = vels[:, 0]
ens['north_vel'][:, k-1] = vels[:, 1]
ens['vert_vel'][:, k-1] = vels[:, 2]
ens['error_vel'][:, k-1] = vels[:, 3]
elif id_ == int('0003', 16): # Echo Intensities
ens['intens'][:, :, k-1] =\
np.array([[st.unpack("B", fd.read(1))[0]
for ui in range(1, 4+1)]
for iu in range(1, hdr['n_deps']+1)])
# Percent good
elif id_ == int('0104', 16):
ens['percent'][:, k-1] =\
np.array([st.unpack("B", fd.read(1))[0]
for ui in range(1, hdr['n_deps']+1)])
# Discharge
elif id_ == int('000C', 16):
if cfg['adcp_type'] is 'narrowband':
if cfg['prog_ver'] < 1.85:
fd.seek(4*hdr['n_deps'], 1)
else:
fd.seek(4*(hdr['n_deps']-1)+2, 1)
else:
if cfg['prog_ver'] > 6:
fd.seek(4*hdr['n_deps'], 1)
else:
fd.seek(4*(hdr['n_deps']-1)+2, 1)
elif id_ == int('000E', 16): # Navigation
# 2 next in m/s
ens['nav_east_vel'][k-1] = st.unpack("h", fd.read(2)
)[0]*.001
ens['nav_north_vel'][k-1] = st.unpack("h", fd.read(2)
)[0]*.001
# 3 next in m
ens['nav_x_disp'][k-1] = st.unpack("i", fd.read(4))[0]*.01
ens['nav_y_disp'][k-1] = st.unpack("i", fd.read(4))[0]*.01
ens['nav_path_len'][k-1] = st.unpack("I", fd.read(4)
)[0]*.01
# deg latitude
ens['latitude'][k-1] = st.unpack("i", fd.read(4)
)[0]*.001/3600
# deg long
ens['longitude'][k-1] = st.unpack("i", fd.read(4)
)[0]*.001/3600
ens['filterwidth'][k-1] = st.unpack("H", fd.read(2))[0]
elif id_ == int('000A', 16):
self.read_cfg(hdr, fd)
else:
print(('Unrecognized ID code: {}').format(hex(id_)))
if n+1 < len(hdr['dat_offsets']):
skp = hdr['dat_offsets'][n+1]-hdr['dat_offsets'][n]
else:
skp = hdr['nbyte']-hdr['dat_offsets'][n]
print((' - Skipping forward {} bytes\n').format(skp))
fd.seek(skp, 1)
fd.seek(2, 1)
if cfg['adcp_type'] is 'broadband':
big_err = np.where(abs(ens['bt_error_vel']) > abs(bad_vald1000))
else:
big_err = (np.array([], dtype=int), np.array([], dtype=int))
ens['bt_east_vel'][(ens['bt_east_vel'] == bad_vald1000)] = np.NaN
ens['bt_east_vel'][big_err] = np.NaN
ens['bt_north_vel'][ens['bt_north_vel'] == bad_vald1000] = np.NaN
ens['bt_north_vel'][big_err] = np.NaN
ens['bt_vert_vel'][ens['bt_vert_vel'] == bad_vald1000] = np.NaN
ens['bt_vert_vel'][big_err] = np.NaN
ens['bt_error_vel'][ens['bt_error_vel'] == bad_vald1000] = np.NaN
ens['bt_error_vel'][big_err] = np.NaN
if cfg['adcp_type'] is 'broadband':
ens['bt_range'][:, np.isnan(ens['bt_error_vel'])] = np.NaN
big_err = 0
ens['east_vel'][(ens['east_vel'] == bad_vald1000)] = np.NaN
ens['east_vel'][big_err] = np.NaN
ens['north_vel'][ens['north_vel'] == bad_vald1000] = np.NaN
ens['north_vel'][big_err] = np.NaN
ens['vert_vel'][ens['vert_vel'] == bad_vald1000] = np.NaN
ens['vert_vel'][big_err] = np.NaN
ens['error_vel'][ens['error_vel'] == bad_vald1000] = np.NaN
ens['error_vel'][big_err] = np.NaN
ens['nav_east_vel'][ens['nav_east_vel'] == bad_vald1000] = np.NaN
ens['nav_east_vel'][big_err] = np.NaN
ens['nav_north_vel'][ens['nav_north_vel'] == bad_vald1000] = np.NaN
ens['nav_north_vel'][big_err] = np.NaN
ens['longitude'][ens['longitude'] == 2147483647*.001/3600] = np.NaN
ens['latitude'][ens['latitude'] == 2147483647*.001/3600] = np.NaN
return ens, hdr
# --------------------------------------
def nmedian(self, *args, **kargs):
# Copied from median but with handling of NaN different.
nargin = len(args)
x = args[0]
if len(x) == 1:
y = x
return y
x = x.astype(float)
window = args[1]
size = np.array([np.shape(x)[op] for op in range(0, len(np.shape(x)))])
if nargin == 2:
dim_ = min(size[size != 1])
dim = min(np.where(size != 1)[0])
if not dim:
dim = 0
elif nargin == 3:
dim = args[2]-1
add_ = (dim+1) - len(size)
if add_ <= 0:
add = []
else:
add = np.ones(1, (dim+1)-len(size))[0]
siz = np.append(size, add)
n = np.int(np.size(x, dim))
# Permute and reshape so that DIM becomes the
# row dimension of a 2-D array
perm = np.append(np.arange(dim, max(len(size)-1, dim)+1),
np.arange(0, dim))
x = np.transpose(x, tuple(perm)).reshape((int(n), int(np.size(x)/n)))
# Sort along first dimension
x = np.sort(x, 0)
size2 = np.array([np.shape(x)[op]
for op in range(0, len(np.shape(x)))])
n1 = size2[0]
n2 = size2[1]
if n1 == 1:
y = x
else:
if n2 == 1:
# kk=sum(finite(x),1)
kk = sum(np.isfinite(x), 0)
if kk > 0:
x1 = x[np.fix((kk-1)/2).astype(int)]
x2 = x[np.fix(kk/2).astype(int)]
x[abs(x-(x1+x2)/2) > window] = np.NaN
x = np.sort(x, 1)
kk = sum(np.isfinite(x), 1)
x[np.isnan(x)] = 0
y = np.NaN
if kk > 0:
y = sum(x)/kk
else:
kk = sum(np.isfinite(x), 0)
ll = kk < n1-2
kk[ll] = 0
x[:, ll] = x[:, ll] * np.NaN
xx = np.transpose(x).reshape((n1*n2, 1))
x1 = xx[np.fix((kk-1) / 2).astype(int) +
1 + np.array(np.arange(0, n2)) * n1-1]
x2 = xx[np.fix(kk / 2).astype(int) +
1 + np.array(np.arange(0, n2)) * n1-1]
x[(abs(x - np.dot(np.ones((n1, 1)),
np.transpose((x1 + x2) / 2.0)
)) > window)] = np.NaN
x = np.sort(x, 0)
kk = sum(np.isfinite(x), 0)
x[np.isnan(x)] = 0
y = np.NaN + np.ones((1, n2))[0]
if any(kk):
tmp_ = x.reshape((n2, n1))
y[np.where(kk > 0)[0]] = np.sum(x[:, np.where(kk > 0)[0]],
axis=0)/kk[kk > 0]
siz[dim] = 1
inverseorder = perm * 0
inverseorder[perm-1] = np.arange(0, len(perm+1))
siz = np.array([int(x) for x in siz])
perm = np.array([int(x) for x in perm])
inverseorder = np.array([int(x) for x in inverseorder])
if len(np.shape(y)) <= 1:
return y
else:
y = y.reshape(tuple(siz[perm]))
y = np.transpose(y, tuple(inverseorder))
return y
# --------------------------------------
def nmean(self, *args, **kargs):
# R_NMEAN Computes the mean of matrix ignoring NaN values
# R_NMEAN(X,DIM) takes the mean along the dimension DIM of X.
nargin = len(args)
x = args[0]
x = x.astype(float)
kk = np.isfinite(x)
x[kk is False] = 0
size = np.array([np.shape(x)[op]
for op in range(0, len(np.shape(x)))])
if nargin == 1:
# Determine which dimension SUM will use
dim_ = min(size[size != 1])
dim = min(np.where(size != 1)[0])
if not dim:
dim = 0
elif nargin == 2:
dim = args[1]-1
if dim+1 > len(size):
y = x
else:
ndat = np.sum(kk, dim)
indat = ndat == 0
ndat[indat] = 1
y = np.sum(x, dim) / ndat
y[indat] = np.NaN
return y
def read_ADCP(self, filename, pathname, num_av=1, num_ping=-1,
ref='none', despike='yes'):
varargin = [filename, num_av, num_ping, 'ref', ref, 'despike', despike]
# Parse input
readbyfile = 1
self.iNFO = {}
self.iNFO.update({'create': 0})
self.iNFO.update({'mtime': []})
self.iNFO.update({'fileptr': []})
self.iNFO.update({'numrec': []})
num_av = 1
nens = -1
ref = 'best'
vels = [0.3, 0.3, 0.3, 2, 50]
# Assumed year
progyear = 2000
if isinstance(varargin[0], str):
readbyfile = 1
name = varargin[0]
varargin.remove(varargin[0])
else:
readbyfile = 0
tlims = varargin[0]
varargin.remove(varargin[0])
lv = len(varargin)
if (lv > 0) &\
(isinstance(varargin[lv-1], str)) &\
(varargin[lv-1] == 'info'):
print('Adding file info to database\n')
self.iNFO['create'] = 1
varargin.remove(varargin[lv-1])
lv = lv-1
if lv >= 1 & (type(varargin[0]) is not str):
num_av = varargin[0]
varargin.remove(varargin[0])
lv = lv-1
if lv >= 1 & (type(varargin[0]) is not str):
nens = varargin[0]
varargin.remove(varargin[0])
lv = lv-1
while len(varargin) > 0:
if varargin[0][0:3] == 'bas':
progyear = varargin[1]
elif varargin[0][0:3] == 'ref':
ref = varargin[1]
elif varargin[0][0:3] == 'des':
if isinstance(varargin[1], str):
if varargin[1] == 'no':
vels = 'no'
else:
vels = [0.3, 0.3, 0.3, 2, 50]
else:
vels = varargin[1]
else:
print('(Unknown command line ' +
'option -> {})'.format(varargin[0]))
del varargin[0: 2]
ref = ref[0: 3]
shipdisp = 10 * 3 * num_av / 2 * 10
print('(Assuming base year is {} (for NB-ADCP))'.format(progyear))
if readbyfile == 1:
naminfo = dir(name)
tutu = os.stat(name)
naminfo = {'name': name,
'date': ti.gmtime(tutu[-2]),
'bytes': tutu[6]}
print(('Opening file {}').format(name))
fd = open(name, 'rb')
# Read first segment with configuration data
cfg = self.read_cfgseg(fd)
pos = fd.tell()
[ens, hdr] = self.read_ensemble(fd, -num_av, cfg)
bytes_ = fd.tell()-pos
fd.seek(pos, 0)
if nens == -1:
nens = np.fix((naminfo['bytes']-pos)/bytes_)-1
print(('\nEstimating {} ensembles in this file\nReducing' +
'by a factor of {}\n').format(nens, num_av))
else:
print(('\nReading {} ensembles in this file\nReducing by' +
'a factor of {}\n').format(nens, num_av))
n = np.fix(nens / num_av)
n_per_file = n
else:
import pickle
with open('Pdatabase', 'rb') as f:
Pdatabase = pickle.load(f)
kstart = sum(Pdatabase['stime'] <= tlims[0])
if sum(Pdatabase.etime <= tlims(1)) == kstart:
kstart = kstart+1
kend = sum(Pdatabase['etime'] <= tlims[1])+1
if sum(Pdatabase['stime'] <= tlims[1])+1 == kend:
kend = kend-1
if kend == kstart:
kkstart = sum(Pdatabase['cfgtimes'][kstart-1] <= tlims[0])
kkend = sum(Pdatabase['cfgtimes'][kstart-1] <= tlims[1])
nens = sum(Pdatabase['numrec'][kstart-1][kkstart-1: kkend-1])
print(('\nOpening file {}\n'
).format(Pdatabase['filename'][kstart-1]))
print(('Using data blocks {} - {}\n'
).format(kkstart, kkend))
else:
kkstart = max(1, sum(Pdatabase['cfgtimes'][kstart-1] <=
tlims[0]))
nens = sum(Pdatabase['numrec'][kstart-1][kkstart-1: -1])
for kk in range(kstart+1, kend+1):
kkend = sum(Pdatabase['cfgtimes'][kk-1] <= tlims[1])
nens = np.append(nens,
sum(Pdatabase['numrec'][kk-1][0: kkend-1])
)
print(('\nOpening first file {}\n'
).format(Pdatabase['filename'][kstart-1]))
print(('Using data blocks {} - {}\n'
).format(kkstart, len(Pdatabase['cfgtimes'][kstart-1])))
n_per_file = np.cumsum(np.fix(nens/num_av))
n = np.int(n_per_file[-1])
print(('\nReading {} ensembles\nReducing by a factor of {}\n'
).format(sum(nens), num_av))
fd = open(Pdatabase['filename'][kstart-1], 'rb')
fd.seek(Pdatabase['cfgptr'][kstart-1][kkstart-1], 0)
cfg = self.read_cfgseg(fd)
pos = fd.tell()
[ens, hdr] = self.read_ensemble(fd, -num_av, cfg)
fd.seek(pos, 0)
if num_av > 1:
if isinstance(vels, str):
print('\n Simple mean used for ensemble averaging\n')
else:
print(('\n Averaging after outlier rejection with' +
'parameters {}\n').format(vels))
nens = int(nens)
n = int(n)
adcp = {}
# Structure to hold all ADCP data (config. data goes to cfg)
adcp.update({'name': ' ', 'bb-adcp': ' ', 'config': cfg,
'day': np.zeros((n), dtype=np.float),
'time': np.zeros((n), dtype=np.float),
'mtime': np.zeros((n), dtype=np.float),
'number': np.zeros((n), dtype=np.float),
'pitch': np.zeros((n), dtype=np.float),
'roll': np.zeros((n), dtype=np.float),
'heading': np.zeros((n), dtype=np.float),
'temperature': np.zeros((n), dtype=np.float),
'bt_east_vel': np.zeros((n), dtype=np.float),
'bt_north_vel': np.zeros((n), dtype=np.float),
'bt_vert_vel': np.zeros((n), dtype=np.float),
'bt_error_vel': np.zeros((n), dtype=np.float),
'bt_range': np.zeros((4, n), dtype=np.float),
'bt_x_disp': np.zeros((n), dtype=np.float),
'bt_y_disp': np.zeros((n), dtype=np.float),
'east_vel': np.zeros((hdr['n_deps'], n)),
'north_vel': np.zeros((hdr['n_deps'], n)),
'vert_vel': np.zeros((hdr['n_deps'], n)),
'error_vel': np.zeros((hdr['n_deps'], n)),
'intens': np.zeros((hdr['n_deps'], 4, n),
dtype=np.float, order='C'),
'nav_east_vel': np.zeros((n), dtype=np.float),
'nav_north_vel': np.zeros((n), dtype=np.float),
'latitude': np.zeros((n), dtype=np.float),
'longitude': np.zeros((n), dtype=np.float),
'comment': ' '*80,
'intens_raw': np.zeros((cfg['n_cells'], 4, nens),
dtype=np.float, order='C'),
'rday': np.zeros((nens), dtype=np.float),
'rtime': np.zeros((nens), dtype=np.float),
'rmtime': np.zeros((nens), dtype=np.float), })
if cfg['adcp_type'] is 'broadband':
adcp['name'] = 'bb-adcp'
elif cfg['adcp_type'] is 'narrowband':
adcp['name'] = 'nb-adcp'
else:
adcp['name'] = 'unknown'
intens_cal = np.dot(np.transpose([np.log10(cfg['range'])]),
[np.array([1, 1, 1, 1])*20]
) + np.dot(np.transpose([cfg['range']]),
[np.array([1, 1, 1, 1]) *
cfg['absorption'] * 2])
# Loop for all records
for k in range(1, int(n) + 1):
if k > int(n_per_file):
kstart = kstart+1
n_per_file[0] = []
print(('\nOpening next file {}\n'
).format(Pdatabase['filename'][kstart-1]))
fd = open(Pdatabase['filename'][kstart-1], 'rb')
fd.seek(Pdatabase['cfgptr'][kstart-1][0], 0)
cfg = self.read_cfgseg(fd)
pos = fd.tell()
[ens, hdr] = self.read_ensemble(fd, -num_av, cfg)
fd.seek(pos, 0)
if k % 50 == 0:
print(('\n{}').format(k * num_av))
print('.')
ens, hdr = self.read_ensemble(fd, num_av, cfg, ens)
if ens['day'][num_av-1] != ens['day'][0]:
adcp['time'][k-1] = np.median(ens['time'] +
(ens['day'] > ens['day'][0]) *
86400)
adcp['rtime'][(k-1) * num_av: k * num_av] =\
ens['time'] + (ens['day'] > ens['day'][0]) * 86400
if adcp['time'][k-1] > 86400:
adcp['time'][k-1] = adcp['time'][k-1] - 86400
adcp['rtime'][(k-1) * num_av: k * num_av] =\
adcp['rtime'][(k-1) * num_av: k * num_av] - 86400
adcp['day'][k-1] = ens['day'][num_av-1]
adcp['rday'][(k-1) * num_av: k * num_av] = ens['day']
else:
adcp['day'][k-1] = ens['day'][0]
adcp['rday'][(k-1) * num_av: k * num_av] = ens['day']
else:
adcp['time'][k-1] = np.median(ens['time'])
adcp['rtime'][(k-1) * num_av: k * num_av] = ens['time']
adcp['day'][k-1] = ens['day'][0]
adcp['rday'][(k-1) * num_av: k * num_av] = ens['day']
if cfg['adcp_type'] is 'narrowband':
adcp['mtime'][k-1] =\
datetime.date.toordinal(datetime.date(progyear, 1, 1)) +\
366 + adcp['time'][k-1] / 86400
adcp['rmtime'][(k-1) * num_av: k * num_av] =\
datetime.date.toordinal(datetime.date(progyear, 1, 1)) +\
366 + adcp['rtime'][(k-1) * num_av: k * num_av] / 86400
else:
adcp['mtime'][k-1] =\
datetime.date.toordinal(datetime.date(progyear, 1, 1)) +\
366 + adcp['day'][k-1] + adcp['time'][k-1] / 86400
adcp['rmtime'][(k-1) * num_av: k * num_av] =\
datetime.date.toordinal(datetime.date(progyear, 1, 1)) +\
366 + adcp['day'][k-1] + adcp['time'][k-1] / 86400
adcp['number'][k-1] = ens['number'][0]
adcp['pitch'][k-1] = np.mean(ens['pitch'])
adcp['roll'][k-1] = np.mean(ens['roll'])
adcp['heading'][k-1] = np.mean(ens['heading'])
adcp['temperature'][k-1] = np.mean(ens['temperature'])
if isinstance(vels, str):
adcp['bt_east_vel'][k-1] = np.nanmean(ens['bt_east_vel'],
axis=0)
adcp['bt_north_vel'][k-1] = np.nanmean(ens['bt_north_vel'],
axis=0)
adcp['bt_vert_vel'][k-1] = np.nanmean(ens['bt_vert_vel'],
axis=0)
adcp['bt_error_vel'][k-1] = np.nanmean(ens['bt_error_vel'],
axis=0)
adcp['bt_range'][:, k-1] = np.nanmean(ens['bt_range'], axis=1)
adcp['bt_x_disp'][k-1] = np.nanmean(ens['bt_x_disp'], axis=0)
adcp['bt_y_disp'][k-1] = np.nanmean(ens['bt_y_disp'], axis=0)
adcp['nav_east_vel'][k-1] = np.nanmean(ens['nav_east_vel'],
axis=0)
adcp['nav_north_vel'][k-1] = np.nanmean(ens['nav_north_vel'],
axis=0)
adcp['latitude'][k-1] = np.nanmean(ens['latitude'], axis=0)
adcp['longitude'][k-1] = np.nanmean(ens['longitude'], axis=0)
ol = np.ones((np.shape(ens['east_vel'])[0], 1))
if ((ref == 'bes') or (ref == 'bot')) &\
np.isfinite(adcp['bt_east_vel'][k-1]):
adcp['east_vel'][:, k-1] =\
np.nanmean(ens['east_vel'] -
ens['bt_east_vel'][ol-1, :], axis=1)
adcp['north_vel'][:, k-1] = \
np.nanmean(ens['north_vel'] -
ens['bt_north_vel'][ol-1, :], axis=1)
adcp['vert_vel'][:, k-1] =\
np.nanmean(ens['vert_vel'] -
ens['bt_vert_vel'][ol-1, :], axis=1)
adcp['error_vel'][:, k-1] =\
np.nanmean(ens['error_vel'] -
ens['bt_error_vel'][ol-1, :], axis=1)
elif ((ref == 'bes') or (ref == 'nav')) &\
np.isfinite(adcp['nav_east_vel'][k-1]):
if ref == 'bes':
print('No Bottom track - using Nav to correct' +
'velocities')
adcp['east_vel'][:, k-1] =\
np.nanmean(ens['east_vel'] +
ens['nav_east_vel'][ol-1, :], axis=1)
adcp['north_vel'][:, k-1] =\
np.nanmean(ens['north_vel'] +
ens['nav_north_vel'][ol-1, :], axis=1)
adcp['vert_vel'][:, k-1] =\
np.nanmean(ens['vert_vel'], axis=1)
adcp['error_vel'][:, k-1] = np.nanmean(ens['error_vel'],
axis=1)
elif ref == 'non':
adcp['east_vel'][:, k-1] = np.nanmean(ens['east_vel'],
axis=1)
adcp['north_vel'][:, k-1] = np.nanmean(ens['north_vel'],
axis=1)
adcp['vert_vel'][:, k-1] = np.nanmean(ens['vert_vel'],
axis=1)
adcp['error_vel'][:, k-1] = np.nanmean(ens['error_vel'],
axis=1)
else:
print('Velocities uncorrected - set to NaN')
adcp['east_vel'][:, k-1] = adcp['east_vel'][:, k-1]*np.NaN
adcp['north_vel'][:, k-1] = adcp['north_vel' +
''][:, k-1]*np.NaN
adcp['vert_vel'][:, k-1] = adcp['vert_vel'][:, k-1]*np.NaN
adcp['error_vel'][:, k-1] = adcp['error_vel' +
''][:, k-1]*np.NaN
else:
adcp['bt_east_vel'][k-1] = self.nmedian(ens['bt_east_vel'],
vels[3])
adcp['bt_north_vel'][k-1] = self.nmedian(ens['bt_north_vel'],
vels[3])
adcp['bt_vert_vel'][k-1] = self.nmedian(ens['bt_vert_vel'],
vels[1])
adcp['bt_error_vel'][k-1] = self.nmedian(ens['bt_error_vel'],
vels[2])
adcp['bt_range'][:, k-1] = self.nmedian(ens['bt_range'],
vels[4], 2)
adcp['bt_x_disp'][k-1] = self.nmedian(ens['bt_x_disp'],
shipdisp)
adcp['bt_y_disp'][k-1] = self.nmedian(ens['bt_y_disp'],
shipdisp)
adcp['nav_east_vel'][k-1] = self.nmedian(ens['nav_east_vel'],
vels[3])
adcp['nav_north_vel'][k-1] = self.nmedian(ens['nav_north_vel'],
vels[3])
adcp['latitude'][k-1] = self.nmedian(ens['latitude'],
shipdisp/111e3)
adcp['longitude'][k-1] = self.nmedian(ens['longitude'],
shipdisp/111e3)
ol = np.ones((np.shape(ens['east_vel'])[0], 1))
if ((ref == 'bes') | (ref == 'bot')) &\
np.isfinite(adcp['bt_east_vel'][k-1]):
adcp['east_vel'][:, k-1] =\
self.nmedian(ens['east_vel'] -
ens['bt_east_vel'][ol-1, :], vels[0], 2)
adcp['north_vel'][:, k-1] =\
self.nmedian(ens['north_vel'] -
ens['bt_north_vel'][ol-1, :], vels[0], 2)
adcp['vert_vel'][:, k-1] =\
self.nmedian(ens['vert_vel'] -
ens['bt_vert_vel'][ol-1, :],
vels[1], 2)
adcp['error_vel'][:, k-1] =\
self.nmedian(ens['error_vel'] -
ens['bt_error_vel'][ol-1, :],
vels[2], 2)
elif ((ref == 'bes') or (ref == 'nav')) &\
(np.isfinite(adcp['nav_east_vel'][k-1])):
if ref == 'bes':
print('No Bottom track - using Nav to correct' +
'velocities')
adcp['east_vel'][:, k-1] =\
self.nmedian(ens['east_vel'] +
ens['nav_east_vel'][ol-1, :],
vels[0], 2)
adcp['north_vel'][:, k-1] =\
self.nmedian(ens['north_vel'] +
ens['nav_north_vel'][ol-1, :],
vels[0], 2)
adcp['vert_vel'][:, k-1] =\
self.nmedian(ens['vert_vel'], vels[1], 2)
adcp['error_vel'][:, k-1] = self.nmedian(ens['error_vel'],
vels[2], 2)
elif ref == 'non':
adcp['east_vel'][:, k-1] = self.nmedian(ens['east_vel'],
vels[0], 2)
adcp['north_vel'][:, k-1] = self.nmedian(ens['north_vel'],
vels[0], 2)
adcp['vert_vel'][:, k-1] = self.nmedian(ens['vert_vel'],
vels[1], 2)
adcp['error_vel'][:, k-1] = self.nmedian(ens['error_vel'],
vels[2], 2)
else:
print('Velocities uncorrected - set to NaN')
adcp['east_vel'][:, k-1] =\
adcp['east_vel'][:, k-1]*np.NaN
adcp['north_vel'][:, k-1] =\
adcp['north_vel'][:, k-1]*np.NaN
adcp['vert_vel'][:, k-1] =\
adcp['vert_vel'][:, k-1]*np.NaN
adcp['error_vel'][:, k-1] =\
adcp['error_vel'][:, k-1]*np.NaN
adcp['intens'][:, :, k-1] = np.nanmean(ens['intens'], axis=2) *\
cfg['intens_scale']
adcp['intens_raw'][:, :, (k-1)*num_av: (k)*num_av] =\
ens['intens'] * cfg['intens_scale']
adcp['TVG'] = intens_cal
fd.close()
# Make or update the database
if self.iNFO['create'] == 1:
if cfg['adcp_type'] is 'broadband':
etime = adcp['mtime'][k-1]
else:
etime =\
datetime.date.toordinal(datetime.date(progyear, 1, 1)) +\
adcp['time'][k-1]/86400
if os.path.isfile('Pdatabase_' + name[0: -4]) is False:
Pdatabase = {}
Pdatabase.update({'cfg': cfg,
'filename': name,
'stime': self.iNFO['mtime'][0],
'etime': etime,
'cfgtimes': self.iNFO['mtime'],
'cfgptr': self.iNFO['fileptr'],
'numrec': self.iNFO['numrec']})
with open('Pdatabase_'+name[0: -4], 'wb') as f:
pickle.dump(Pdatabase, f)
else:
with open('Pdatabase_' + name[0: -4], 'rb') as f:
Pdatabase = pickle.load(f)
try:
tmp_ = [index
for index, value in enumerate(Pdatabase['stime'])
if value <= self.iNFO['mtime'][0]]
k = sum(tmp_)
except TypeError:
k = 1
return Parameters(adcp)