1#!/usr/bin/env python
+ 2# -*- coding: utf-8 -*-
+ 3"""Module velosearaptor.madcp with functions for moored ADCPs.
+ 4
+ 5### General Use
+ 6Processing a raw ADCP file from a moored deployment is a two-step process.
+ 7- Instantiate a processing object with `ProcessADCP` or `ProcessADCPyml`. The
+ 8 former expects the path to the raw data and a number of dictionaries with
+ 9 processing parameters as input. The latter reads a .yml file containing the
+ 10 processing parameters. See the respective docstrings for more information.
+ 11- Process raw pings and possibly run a ping-averaging method on the data.
+ 12 Options here are `ProcessADCP.process_pings`,
+ 13 `ProcessADCP.average_ensembles`, and `ProcessADCP.burst_average_ensembles`.
+ 14
+ 15### Notes
+ 16Some general notes for this module.
+ 17
+ 18#### Depth Gridding
+ 19The depth vector for the ADCP raw data (in instrument coordinates) is
+ 20calculated in :meth:`pycurrents.rdiraw.FileBBWHOS` as <br>
+ 21`dep = np.arange(NCells) * CellSize + Bin1Dist`
+ 22The depth vector thus points to the center of each bin.
+ 23
+ 24From the RDI manual *WorkHorse Monitor, Sentinel, Mariner, Quartermaster, and
+ 25Long Ranger ADCPs Commands and Output Data Format*:
+ 26
+ 27> This [Bin1Dist] field contains the distance to the middle of the first depth
+ 28> cell (bin). This distance is a function of depth cell length (WS), the
+ 29> profiling mode (WM), the blank after transmit distance (WF), and speed of
+ 30> sound.
+ 31
+ 32"""
+ 33
+ 34import datetime
+ 35import logging
+ 36import os
+ 37import pathlib
+ 38from pathlib import Path
+ 39from shutil import which
+ 40from subprocess import PIPE, Popen # for magdec
+ 41from warnings import warn
+ 42
+ 43import gsw
+ 44import matplotlib.pyplot as plt
+ 45import numpy as np
+ 46import scipy as sp
+ 47import xarray as xr
+ 48from pycurrents.adcp.rdiraw import Multiread
+ 49from pycurrents.adcp.transform import Transform, rdi_xyz_enu
+ 50from pycurrents.codas import to_date, to_day
+ 51from pycurrents.data import seawater
+ 52from pycurrents.num import interp1
+ 53from pycurrents.num.nptools import rangeslice
+ 54from pycurrents.system import Bunch
+ 55from tqdm import tqdm
+ 56
+ 57from . import io
+ 58
+ 59# Standard logging
+ 60logger = logging.getLogger(__name__)
+ 61
+ 62
+ 63class ProcessADCP:
+ 64 """Moored ADCP Processing.
+ 65
+ 66 An instance of ProcessADCP is initialized by providing raw data and
+ 67 processing parameters. Once initialized, the instance method
+ 68 :meth:`average_ensembles` can be used to average over just a few or
+ 69 all pings.
+ 70
+ 71 Magnetic declination is automatically calculated via a call to the command
+ 72 line tool
+ 73 [magdec](https://currents.soest.hawaii.edu/hgstage/geomag/file/tip) that
+ 74 must be installed. Once calculated, the magnetic declination is stored in
+ 75 `magdec` and automatically applied to the data.
+ 76
+ 77 Prior to time averaging, data are edited based on correlation and error
+ 78 velocity thresholds. The variable `pg` is calculated based on the number of
+ 79 excluded pings, i.e. it is the numbur of good pings divided by the number
+ 80 of total pings within one depth bin and one time bin. It may be used to
+ 81 further filter the data.
+ 82
+ 83 A pdf document on ADCP data collection and processing principles can be
+ 84 downloaded from RDI
+ 85 [here](https://www.comm-tec.com/Docs/Manuali/RDI/BBPRIME.pdf).
+ 86
+ 87 Time-averaged data are grouped together in an `xarray.Dataset` in the
+ 88 instance attribute `ds` for easy access and convenient output to netcdf
+ 89 format.
+ 90
+ 91 Parameters
+ 92 ----------
+ 93 raw_data : str or list or Path
+ 94 Location(s) of raw data.
+ 95 meta_data : dict
+ 96 Dictionary with meta data. At a minimum entries for `lon` and `lat` are
+ 97 needed. If `mooring` and `sn` provide mooring name and serial number
+ 98 then these will be used to name the log file produced during
+ 99 processing.
+ 100 driftparams : dict, optional
+ 101 Time drift parameters. See notes below.
+ 102 tgridparams : dict, optional
+ 103 Time gridding parameters. See notes below.
+ 104 dgridparams : dict, optional
+ 105 Depth gridding parameters. See notes below.
+ 106 editparams : dict, optional
+ 107 Editing parameters. See notes below.
+ 108 ibad : int, optional
+ 109 Mark beam with bad data (zero based). Defaults to None.
+ 110 logdir : str, optional
+ 111 Log file directory. Defaults to `log/`.
+ 112 verbose : bool, optional
+ 113 Output more processing info to screen.
+ 114 magdec : float, optional
+ 115 Magnetic declination in degrees.
+ 116 pressure : xr.DataArray, optional
+ 117 Pressure time series in dbar as xr.DataAarray with coordinate `time`.
+ 118
+ 119 Attributes
+ 120 ----------
+ 121 files : list
+ 122 List pointing to raw data file(s).
+ 123 dday_start : float
+ 124 Start time of ADCP time series, determined either through `t0` in
+ 125 `tgridparams` or the start of the time series once at depth.
+ 126 dday_end : float
+ 127 End time of ADCP time series, determined either through `t1` in
+ 128 `tgridparams` or the end of the time series at depth.
+ 129 start_ddays : list
+ 130 Start times of averaging intervals
+ 131 dday_mid : float
+ 132 Time stamp for ping average as determined in :meth:`make_start_ddays`.
+ 133 dt : float
+ 134 Time that is inclusive of one average. For regular averaging, this is
+ 135 `dt_hours` in `tgridparams` converted to Julian days. For
+ 136 burst-averagint, this is a time interval that is inclusive of all pings
+ 137 in one bursts (but goes slightly beyond that into the time between
+ 138 bursts).
+ 139 time_drift_rate : float
+ 140 Clock drift calculated from `driftparams`.
+ 141 orientation : str
+ 142 Instrument orientation `up` or `down`.
+ 143 magdec : float
+ 144 Magnetic declination.
+ 145 m : pycurrents.adcp.rdiraw.Multiread
+ 146 Multiread instance.
+ 147 tsdat : Bunch
+ 148 Auxiliary data.
+ 149 raw : xarray.Dataset
+ 150 Raw ADCP data.
+ 151 ds : xarray.Dataset
+ 152 Time-averaged dataset added by running :meth:`average_ensembles`.
+ 153
+ 154 Notes
+ 155 -----
+ 156 Various parameters are passed to the instance through dictionaries. Their
+ 157 specifics are described below. Gridding and editing parameters can be
+ 158 updated after creating the ProcessADCP instance via `parse_dgridparams`,
+ 159 `parse_tgridparams`, and `parse_editparams`.
+ 160
+ 161 **Time drift parameters**
+ 162 Provide clock drift parameters via `driftparams`. Accepted entries are
+ 163 - `end_adcp` : Time of ADCP at data download
+ 164 - `end_pc` : UTC time at data download
+ 165
+ 166 The difference between `end_adcp` and `end_pc` is used to linearly correct
+ 167 for instrument clock drift.
+ 168
+ 169 **Time gridding parameters**
+ 170 Provide time gridding parameters via `tgridparams`. Accepted entries are
+ 171 - `dt_hours` : Time grid interval. Defaults to 0.5h.
+ 172 - `t0` : Start time for gridding. Determined from data if not provided.
+ 173 - `t1` : End time for gridding. Determined from data if not provided.
+ 174 - burst_average : bool
+ 175 Set ensemble averaging to act on burst sampling scheme. Defaults to False.
+ 176
+ 177 **Depth gridding parameters**
+ 178 Provide depth gridding parameters via `dgridparams`. Accepted entries are
+ 179 - `dtop` : Shallow depth in m.
+ 180 - `dbot` : Deep depth in m.
+ 181 - `dinterval` : Vertical grid size in m. Defaults to 5m.
+ 182
+ 183 Values for `dbot` and `dtop` are generated if not provided.
+ 184
+ 185 **Editing parameters**
+ 186 Provide editing parameters via `editparams`.
+ 187 - `max_e`=0.2, # absolute max e
+ 188 - `max_e_deviation`=2, # max in terms of sigma
+ 189 - `min_correlation`=64, # 64 is RDI default
+ 190 - `maskbins` : Array with booleans indexing into the ADCP bins. Use the
+ 191 convenience method `generate_binmask`.
+ 192 - `pg_limit` : float or int or None.
+ 193 Percent good limit applied prior to interpolating to the universal
+ 194 depth grid in `burst_average_ensembles`.
+ 195
+ 196 """
+ 197
+ 198 # Default editing parameters.
+ 199 _editparams = dict(
+ 200 max_e=0.2, # absolute max e
+ 201 max_e_deviation=2, # max in terms of sigma
+ 202 min_correlation=64, # 64 is RDI default
+ 203 maskbins=None, # do not mask any bins
+ 204 pg_limit=50, # percent good limit applied in `burst_average_ensembles`
+ 205 )
+ 206
+ 207 def __init__(
+ 208 self,
+ 209 raw_data,
+ 210 meta_data,
+ 211 driftparams=None,
+ 212 tgridparams=None,
+ 213 dgridparams=None,
+ 214 editparams=None,
+ 215 ibad=None,
+ 216 logdir="log",
+ 217 verbose=False,
+ 218 plot=False,
+ 219 pressure_scale_factor=1,
+ 220 magdec=None,
+ 221 pressure=None,
+ 222 ):
+ 223 self.meta_data = Bunch(meta_data.copy())
+ 224 self.ibad = ibad
+ 225 self.logdir = logdir
+ 226 self.verbose = verbose
+ 227
+ 228 self._magdec_provided = magdec
+ 229 self._magdec = magdec
+ 230
+ 231 self._pressure_provided = pressure
+ 232 self._pressure_scale_factor = pressure_scale_factor
+ 233
+ 234 self._raw = None
+ 235 self._default_dgridparams = None
+ 236
+ 237 self.parse_file_locations(raw_data)
+ 238 self._initiate_data_reader()
+ 239 self._read_auxiliary_data()
+ 240 self._parse_meta_data()
+ 241
+ 242 self._set_up_logger()
+ 243
+ 244 self.parse_driftparams(driftparams)
+ 245 self._parse_sysconfig()
+ 246 self.parse_dgridparams(dgridparams)
+ 247 self.parse_tgridparams(tgridparams)
+ 248 self.parse_editparams(editparams)
+ 249
+ 250 self.make_start_ddays()
+ 251
+ 252 if plot:
+ 253 self.plot_pressure()
+ 254
+ 255 def parse_file_locations(self, raw_data, min_file_size=1e4):
+ 256 """Parse input for raw data files.
+ 257
+ 258 Input can either be a single file name as a str, a single file as a
+ 259 Path instance, a list of either of these, or a Path instance pointing to a
+ 260 directory with raw ADCP files. In the latter case, files that are
+ 261 smaller than a threshold will not be included in the processing.
+ 262
+ 263 Outputs to attribute `files`. The output can be fed to Multiread instances.
+ 264
+ 265 Parameters
+ 266 ----------
+ 267 raw_data : str or list or Path
+ 268 Location(s) of raw data.
+ 269 min_file_size : int
+ 270 Minimum size for file to be included. Defaults to 1e4 which
+ 271 corresponds to about 10kB and is a good value for excluding small
+ 272 files without any actual data.
+ 273
+ 274 """
+ 275
+ 276 def list_dir(dir, min_file_size):
+ 277 # List all raw files.
+ 278 all_raw_files = list(sorted(dir.glob("*.00*")))
+ 279 # only files larger than about 10kB
+ 280 files = [
+ 281 file.as_posix()
+ 282 for file in all_raw_files
+ 283 if file.stat().st_size > min_file_size
+ 284 ]
+ 285 return files
+ 286
+ 287 input_type = type(raw_data)
+ 288 if input_type is list:
+ 289 if type(raw_data[0]) is str:
+ 290 self.files = raw_data
+ 291 elif type(raw_data[0]) is pathlib.PosixPath:
+ 292 self.files = [file.as_posix() for file in raw_data]
+ 293 elif input_type is str:
+ 294 if Path(raw_data).is_dir():
+ 295 self.files = list_dir(Path(raw_data), min_file_size)
+ 296 else:
+ 297 self.files = [raw_data]
+ 298 elif input_type is pathlib.PosixPath:
+ 299 if raw_data.is_dir():
+ 300 self.files = list_dir(raw_data, min_file_size)
+ 301 else:
+ 302 self.files = [raw_data.as_posix()]
+ 303
+ 304 def _initiate_data_reader(self):
+ 305 """Initiate a Multiread data reader.
+ 306
+ 307 Adds attribute `m`.
+ 308
+ 309 """
+ 310 self.m = Multiread(self.files, sonar="wh", ibad=self.ibad)
+ 311 # Make some more meta data realily available by reading a single ping
+ 312 # from the raw data.
+ 313 ping = self.m.read(start=0, stop=1)
+ 314 self.meta_data.Bin1Dist = ping.FL.Bin1Dist / 100.0
+ 315 self.meta_data.NCells = ping.FL.NCells
+ 316 self.meta_data.CellSize = ping.FL.CellSize / 100.0
+ 317
+ 318 def _generate_external_pressure_interpolator(self, dat):
+ 319 """Interpolate external pressure to pycurrents time vector.
+ 320
+ 321 Parameters
+ 322 ----------
+ 323 dat : pycurrents.adcp.rdiraw.Bunch
+ 324 Data structure.
+ 325
+ 326 Returns
+ 327 -------
+ 328 Interpolator
+ 329
+ 330 """
+ 331
+ 332 # We already have a function to convert the pycurrents time to
+ 333 # datetime64. Interpolate in this time domain.
+ 334 time = io.yday0_to_datetime64(dat.yearbase, dat.dday)
+ 335 p_interpolated = self._pressure_provided.interp(time=time).data
+ 336
+ 337 # Beginning and end may have NaN's if the ADCP was started before
+ 338 # the external pressure sensor. Make sure this happens only when
+ 339 # outside the water and then replace with atmospheric pressure.
+ 340 if np.any(np.isnan(p_interpolated)):
+ 341 nan_mask = np.isnan(p_interpolated)
+ 342 i_nan = np.flatnonzero(nan_mask)
+ 343
+ 344 first_few_good_median = np.median(p_interpolated[~nan_mask][:20])
+ 345 last_few_good_median = np.median(p_interpolated[~nan_mask][-20:])
+ 346
+ 347 if np.isnan(p_interpolated[0]):
+ 348 if first_few_good_median < 1:
+ 349 i_divide = np.flatnonzero(np.diff(i_nan) - 1) + 1
+ 350 if i_divide.size == 0:
+ 351 p_interpolated[nan_mask] = first_few_good_median
+ 352 else:
+ 353 p_interpolated[0:i_divide] = first_few_good_median
+ 354
+ 355 if np.isnan(p_interpolated[-1]):
+ 356 if last_few_good_median < 1:
+ 357 i_divide = np.flatnonzero(np.diff(i_nan) - 1) + 1
+ 358 if i_divide.size == 0:
+ 359 p_interpolated[nan_mask] = last_few_good_median
+ 360 else:
+ 361 p_interpolated[i_divide:] = last_few_good_median
+ 362
+ 363 # Now generate an interpolation function that will take dday as input
+ 364 # for later per-ensemble interpolation.
+ 365 self._external_pressure_interpolator = sp.interpolate.interp1d(
+ 366 dat.dday,
+ 367 p_interpolated,
+ 368 bounds_error=False,
+ 369 fill_value="extrapolate",
+ 370 )
+ 371
+ 372 def _external_pressure_to_dat(self, dat):
+ 373 """Interpolate external pressure to pycurrents time vector.
+ 374
+ 375 Parameters
+ 376 ----------
+ 377 dat : pycurrents.adcp.rdiraw.Bunch
+ 378 Data structure.
+ 379
+ 380 Returns
+ 381 -------
+ 382 array-like
+ 383 Pressure
+ 384
+ 385 """
+ 386 return self._external_pressure_interpolator(dat.dday)
+ 387
+ 388 def _scale_pycurrents_pressure(self, dat):
+ 389 # Initial pressure units: 10 Pa (about 1 mm or 0.001 decibar).
+ 390 # Converting to decibars.
+ 391 return dat.VL["Pressure"] / 1000.0 * self._pressure_scale_factor
+ 392
+ 393 def _read_auxiliary_data(self):
+ 394 """Read auxiliary data.
+ 395
+ 396 Adds attribute `tsdat`.
+ 397
+ 398 """
+ 399 tsdat = self.m.read(varlist=["VariableLeader"])
+ 400 tsdat.temperature = tsdat.VL["Temperature"] / 100.0
+ 401 # Replace pressure if provided from external sensor
+ 402 if self._pressure_provided is not None:
+ 403 self._generate_external_pressure_interpolator(tsdat)
+ 404 tsdat.pressure = self._external_pressure_to_dat(tsdat)
+ 405 else:
+ 406 tsdat.pressure = self._scale_pycurrents_pressure(tsdat)
+ 407 self.tsdat = tsdat
+ 408
+ 409 def _parse_meta_data(self):
+ 410 """Parse meta data.
+ 411
+ 412 - Add essential meta data to attributes. Will throw a KeyError if no
+ 413 lon/lat provided.
+ 414
+ 415 - Read serial number from raw data and complain if it does not match
+ 416 the meta data SN.
+ 417 """
+ 418 essential_meta_data = ["lon", "lat"]
+ 419
+ 420 [
+ 421 self._safely_add_attribute_from_params(k, self.meta_data)
+ 422 for k in essential_meta_data
+ 423 ]
+ 424
+ 425 # Check SN
+ 426 sn_internal = int.from_bytes(self.tsdat.FL.Inst_SN, "little")
+ 427
+ 428 if "sn" not in self.meta_data:
+ 429 self.meta_data.sn = sn_internal
+ 430
+ 431 # Check internal SN matches user set one
+ 432 if sn_internal != self.meta_data.sn:
+ 433 warn(
+ 434 f"Serial number in file, {sn_internal}, is different from that set by user, {self.meta_data.sn}. Keeping user value."
+ 435 )
+ 436
+ 437 @property
+ 438 def default_dgridparams(self):
+ 439 """Determine default depth gridding parameters.
+ 440
+ 441 The grid is centered on the median depth of the ADCP (plus distance to
+ 442 the center of the first bin) to avoid unnecessary binning into
+ 443 neighboring depth cells. The default size of the depth bins mimicks the
+ 444 size of ADCP bins.
+ 445 """
+ 446 if self._default_dgridparams is None:
+ 447 # Only use pressure at depth, not on deck
+ 448 ii = np.flatnonzero(self.tsdat.pressure > 15)
+ 449 p = self.tsdat.pressure[ii]
+ 450 # Determine limits of the pressure distribution but leave out the
+ 451 # top 5 and bottom 2 percent of data points. This way we are hoping
+ 452 # to avoid any outliers and a possible pressure record of ascent
+ 453 # and/or descent.
+ 454 p_top, p_bot = np.round(
+ 455 seawater.depth2(np.percentile(p, [5, 98]), self.lat)
+ 456 )
+ 457 self.p_median = np.median(p)
+ 458 pdep_median = np.round(seawater.depth2(self.p_median, self.lat))
+ 459 n = self.meta_data.NCells + 2
+ 460 d_interval = self.meta_data.CellSize
+ 461 distance_to_first_bin = np.round(self.meta_data.Bin1Dist)
+ 462 distance_to_last_bin = distance_to_first_bin + n * d_interval
+ 463
+ 464 if self.sysconfig["up"]:
+ 465 # Set minimum grid depth level. Anything shallower than 10m
+ 466 # will be garbage anyways so let's throw this out.
+ 467 dtop = p_top - distance_to_last_bin
+ 468 if dtop < 10:
+ 469 dtop = 10
+ 470 dbot = pdep_median - distance_to_first_bin
+ 471 # Successively add bins until we reach maximum pressure. This
+ 472 # will take care of mooring knockdowns.
+ 473 while dbot < p_bot:
+ 474 dbot += d_interval
+ 475 else:
+ 476 dtop = pdep_median + distance_to_first_bin
+ 477 while dtop > p_top:
+ 478 dtop -= d_interval
+ 479 dbot = p_bot + distance_to_last_bin
+ 480
+ 481 self._default_dgridparams = dict(
+ 482 dtop=dtop,
+ 483 dbot=dbot,
+ 484 d_interval=d_interval,
+ 485 )
+ 486
+ 487 return self._default_dgridparams
+ 488
+ 489 def parse_dgridparams(self, dgridparams):
+ 490 """Parse depth gridding parameters.
+ 491
+ 492 See top level class notes for more info.
+ 493
+ 494 Parameters
+ 495 ----------
+ 496 dgridparams : dict
+ 497
+ 498 """
+ 499 self.dgridparams = Bunch(self.default_dgridparams)
+ 500 if dgridparams is not None:
+ 501 self.dgridparams.update_values(dgridparams, strict=True)
+ 502 else:
+ 503 logger.warning(
+ 504 "No depth gridding parameters provided, using default values."
+ 505 )
+ 506 # Default depth grid parameters are based on the median pressure to
+ 507 # avoid binning into neighboring grid cells as much as possible.
+ 508 # Therefore, we start assembling the depth grid from the bottom up for
+ 509 # an uplooker and from the top down for a downlooker.
+ 510 if self.orientation == "up":
+ 511 self.dgrid = np.arange(
+ 512 self.dgridparams.dbot,
+ 513 self.dgridparams.dtop,
+ 514 -self.dgridparams.d_interval,
+ 515 dtype=float,
+ 516 )
+ 517 elif self.orientation == "down":
+ 518 self.dgrid = np.arange(
+ 519 self.dgridparams.dtop,
+ 520 self.dgridparams.dbot,
+ 521 self.dgridparams.d_interval,
+ 522 dtype=float,
+ 523 )
+ 524
+ 525 def parse_tgridparams(self, tgridparams):
+ 526 """Parse time gridding parameters.
+ 527
+ 528 See top level class notes for more info.
+ 529
+ 530 Parameters
+ 531 ----------
+ 532 tgridparams : dict
+ 533
+ 534 """
+ 535 # Find time at depth to determine default time grid parameters.
+ 536 # Differentiate between time series only in the water and time series
+ 537 # including the overshoot on mooring deployment.
+ 538 p = self.tsdat.pressure
+ 539 if ~np.any(p < 10):
+ 540 t0 = self.dday[0]
+ 541 t1 = self.dday[-1]
+ 542 else:
+ 543 at_depth = np.nonzero(p > self.p_median)[0][0]
+ 544 t0 = self.dday[at_depth]
+ 545 in_water = np.nonzero(p > self.p_median / 2)[0][-1]
+ 546 t1 = self.dday[in_water]
+ 547
+ 548 # Generate a set of default time gridding parameters and then update
+ 549 # from the input parameters provided.
+ 550 default_tgridparams = dict(dt_hours=0.5, t0=t0, t1=t1, burst_average=False)
+ 551 self.tgridparams = Bunch(default_tgridparams)
+ 552 if tgridparams is not None:
+ 553 self.tgridparams.update_values(tgridparams, strict=True)
+ 554 else:
+ 555 logger.warning(
+ 556 "No time gridding parameters provided, using default values."
+ 557 )
+ 558
+ 559 def parse_editparams(self, editparams):
+ 560 """Parse editing parameters.
+ 561
+ 562 See top level class notes for more info.
+ 563
+ 564 Parameters
+ 565 ----------
+ 566 editparams : dict
+ 567
+ 568 """
+ 569 self.editparams = Bunch(self._editparams)
+ 570 if editparams is not None:
+ 571 self.editparams.update_values(editparams, strict=True)
+ 572 else:
+ 573 logger.warning("No edit parameters provided, using default values.")
+ 574
+ 575 def parse_driftparams(self, driftparams):
+ 576 """Parse time drift parameters.
+ 577
+ 578 See top level class notes for more info.
+ 579
+ 580 Parameters
+ 581 ----------
+ 582 driftparams : dict
+ 583
+ 584 """
+ 585 driftparams = dict() if driftparams is None else driftparams
+ 586 self.driftparams = driftparams
+ 587 self.yearbase = self.m.yearbase
+ 588 t0 = self.tsdat.dday[0]
+ 589 self.t0 = t0
+ 590 t1_adcp = driftparams.get("end_adcp", None)
+ 591 if t1_adcp is not None:
+ 592 t1_pc = to_day(self.m.yearbase, *driftparams["end_pc"])
+ 593 t1_adcp = to_day(self.m.yearbase, *driftparams["end_adcp"])
+ 594 self.time_drift_rate = (t1_pc - t0) / (t1_adcp - t0)
+ 595 else:
+ 596 logger.warning(
+ 597 "No time drift parameters provided, not applying any clock correction."
+ 598 )
+ 599 self.time_drift_rate = 1
+ 600
+ 601 self.dday = self._correct_dday(self.tsdat.dday)
+ 602
+ 603 def _correct_dday(self, dday_orig):
+ 604 """Apply linear correction for clock drift.
+ 605
+ 606 Parameters
+ 607 ----------
+ 608 dday_orig : array-like
+ 609 Origial time vector.
+ 610
+ 611 Returns
+ 612 -------
+ 613 array-like
+ 614 Corrected time vector
+ 615
+ 616 """
+ 617
+ 618 return self.t0 + self.time_drift_rate * (dday_orig - self.t0)
+ 619
+ 620 def _parse_sysconfig(self):
+ 621 # We have to get the up/down reading from sysconfig for a time when the
+ 622 # instrument was in the water. Look at pressure and determine ensembles
+ 623 # during time at depth.
+ 624 ii = np.flatnonzero(self.tsdat.pressure > 15)
+ 625 depth_ii = (ii[0] + ii[-1]) // 2
+ 626 try:
+ 627 self.tsdat.pressure[depth_ii] > 15
+ 628 except ValueError:
+ 629 print("could not determine ensemble index at depth for reading sysconfig")
+ 630 at_depth = self.m.read(
+ 631 varlist=["VariableLeader"], start=depth_ii, stop=depth_ii + 1
+ 632 )
+ 633 self.orientation = "up" if at_depth.sysconfig.up else "down"
+ 634 self.sysconfig = at_depth.sysconfig
+ 635
+ 636 def make_start_ddays(self):
+ 637 """Generate time stamps for ping averaging.
+ 638
+ 639 Notes
+ 640 -----
+ 641 If turning on burst averaging, other input values will be ignored and
+ 642 the averaging interval will be determined from the burst sampling
+ 643 scheme apparent in the ping pattern.
+ 644
+ 645 """
+ 646 dday_start = self.tgridparams.t0
+ 647 dday_end = self.tgridparams.t1
+ 648 dt_hours = self.tgridparams.dt_hours
+ 649 is_burst_average = self.tgridparams.burst_average
+ 650
+ 651 # Save whether we are averaging over bursts or not.
+ 652 self.is_burst_average = is_burst_average
+ 653 # Generate time stamps and stuff.
+ 654 if not is_burst_average:
+ 655 print("no burst average")
+ 656 self.dday_start = dday_start
+ 657 self.dday_end = dday_end
+ 658 self.dt = dt_hours / 24.0
+ 659 self.start_ddays = np.arange(dday_start, dday_end, dt_hours / 24.0)
+ 660 # Time stamps for the averages. Midpoints of averaging intervals.
+ 661 self.dday_mid = self.start_ddays + self.dt / 2
+ 662 else:
+ 663 dday_diff = np.diff(self.dday)
+ 664 # Determine ping interval within burst and time between bursts.
+ 665 burst_dt = np.median(dday_diff)
+ 666 print(f"time between pings within burst: {burst_dt * 24 * 60 * 60:1.1f} s")
+ 667 # It seems safe to assume that the time between bursts is at least
+ 668 # four times as long as the time between individual pings within a
+ 669 # burst.
+ 670 inter_burst_dt = np.median(dday_diff[dday_diff > burst_dt * 4])
+ 671 print(f"time between bursts: {inter_burst_dt * 24 * 60:1.1f} min")
+ 672 # Find starting points of all bursts.
+ 673 start_indices = np.flatnonzero(dday_diff > burst_dt * 4)
+ 674 # Increase index so we are at the end of the larger time differences.
+ 675 start_indices += 1
+ 676 # Include the very beginning.
+ 677 start_indices = np.insert(start_indices, 0, 0)
+ 678 self.start_ddays = self.dday[start_indices]
+ 679 self.dday_start = self.start_ddays[0]
+ 680 # Generate a dt that is inclusive of one burst. We know the number
+ 681 # of pings in a burst from the difference between the
+ 682 # start_indices. Let's go a bit beyond the time needed (3 more ping
+ 683 # intervals chosen here).
+ 684 pings_per_burst = np.int32(np.median(np.diff(start_indices)))
+ 685 print(f"{pings_per_burst} pings per burst")
+ 686 print(f"{start_indices.shape[0]} bursts total")
+ 687 self.dt = burst_dt * pings_per_burst + burst_dt * 3
+ 688
+ 689 self.dday_start = self.start_ddays[0]
+ 690 self.dday_end = dday_end
+ 691
+ 692 # Time stamps in the middle of the burst
+ 693 self.dday_mid = self.start_ddays + pings_per_burst * burst_dt / 2
+ 694
+ 695 def read_ensemble(self, iens):
+ 696 """Read ensembles (several individual pings grouped together).
+ 697
+ 698 Parameters
+ 699 ----------
+ 700 iens : int
+ 701 Index into ensemble start times (they are generated in
+ 702 :meth:`make_start_ddays`).
+ 703
+ 704 Returns
+ 705 -------
+ 706 dat : pycurrents.adcp.rdiraw.Bunch
+ 707 Dictionary with data.
+ 708
+ 709 Raises
+ 710 ------
+ 711 ValueError
+ 712 If `iens` is too large to index into `start_ddays`.
+ 713
+ 714 """
+ 715
+ 716 if iens > len(self.start_ddays) - 1:
+ 717 raise ValueError("ens num %d is out of range" % iens)
+ 718
+ 719 # Get indices within the interval.
+ 720 sl = rangeslice(
+ 721 self.dday, self.start_ddays[iens], self.start_ddays[iens] + self.dt
+ 722 )
+ 723 # Use the indices to extract data.
+ 724 dat = self.m.read(start=sl.start, stop=sl.stop)
+ 725 if dat is None:
+ 726 return None
+ 727 dat.dday_orig = dat.dday
+ 728 dat.dday = self._correct_dday(dat.dday_orig)
+ 729 if self._pressure_provided is not None:
+ 730 # Replace pressure if provided
+ 731 dat.pressure = self._external_pressure_to_dat(dat)
+ 732 else:
+ 733 dat.pressure = self._scale_pycurrents_pressure(dat)
+ 734 sign = -1 if self.orientation == "up" else 1
+ 735 pdepth = seawater.depth2(dat.pressure, self.lat)
+ 736 dat.depth = pdepth[:, np.newaxis] + sign * dat.dep
+ 737 return dat
+ 738
+ 739 @property
+ 740 def magdec(self):
+ 741 """Magnetic declination.
+ 742
+ 743 If not provided as input argument magdec is calculated using
+ 744 [magdec](https://currents.soest.hawaii.edu/hgstage/geomag/file/tip)
+ 745 (must be installed) based on `lon` and `lat`.
+ 746
+ 747 """
+ 748 if self._magdec is None:
+ 749 if self.lat is None:
+ 750 logger.warning(
+ 751 "No lon/lat provided, cannot calculate magnetic declination."
+ 752 )
+ 753 self._magdec = 0
+ 754 else:
+ 755 # Look for magdec executable
+ 756 magdec_found = True
+ 757 magdec_path = which("magdec")
+ 758
+ 759 if magdec_path is None:
+ 760 magdec_found = False
+ 761 package_dir = os.path.dirname(__file__)
+ 762
+ 763 if not magdec_found:
+ 764 # Try this package directory
+ 765 magdec_path = os.path.join(package_dir, "magdec")
+ 766 magdec_found = os.path.isfile(magdec_path)
+ 767
+ 768 if not magdec_found:
+ 769 # Try the magdec installation directory
+ 770 magdec_path = os.path.abspath(
+ 771 os.path.join(package_dir, "../geomag/magdec")
+ 772 )
+ 773 magdec_found = os.path.isfile(magdec_path)
+ 774
+ 775 if not magdec_found:
+ 776 raise FileNotFoundError(
+ 777 "Cannot find program magdec on the system path or paths within velosearaptor."
+ 778 )
+ 779
+ 780 logger.info(f"magdec found at {magdec_path}")
+ 781
+ 782 n = len(self.start_ddays)
+ 783 dday_mid = self.start_ddays[n // 2]
+ 784 y, m, d = to_date(self.yearbase, dday_mid)[:3]
+ 785
+ 786 output = Popen(
+ 787 [
+ 788 magdec_path,
+ 789 str(self.lon),
+ 790 str(self.lat),
+ 791 str(y),
+ 792 str(m),
+ 793 str(d),
+ 794 ],
+ 795 stdout=PIPE,
+ 796 ).communicate()[0]
+ 797 output = output.strip()
+ 798 logger.info("magdec output is: %s", output)
+ 799 self._magdec = float(output.split()[0])
+ 800 elif self._magdec_provided is not None:
+ 801 logger.info(f"magdec {self._magdec_provided} provided.")
+ 802 return self._magdec
+ 803
+ 804 @property
+ 805 def raw(self):
+ 806 """Raw ADCP data."""
+ 807 if self._raw is None:
+ 808 print("Reading raw data...")
+ 809 self._raw = io.read_raw_rdi(self.files)
+ 810 self._raw.coords["bin"] = (("z"), np.arange(self._raw.z.size))
+ 811 return self._raw
+ 812
+ 813 def _edit(self, ens):
+ 814 """Apply editing to xyze."""
+ 815 ep = self.editparams
+ 816 cond = (ens.cor < ep.min_correlation).any(axis=-1)
+ 817 ens.xyze[cond] = np.ma.masked
+ 818 e = ens.xyze[:, :, 3]
+ 819 max_e = min(ep.max_e, e.std() * ep.max_e_deviation)
+ 820 ens.max_e_applied = max_e
+ 821 cond = np.abs(e) > max_e
+ 822 ens.xyze[cond] = np.ma.masked
+ 823 if ep.maskbins is not None:
+ 824 ens.xyze[:, ep.maskbins, :] = np.ma.masked
+ 825
+ 826 def _to_enu(self, ens):
+ 827 """
+ 828 add enu
+ 829 GV: enu is east, north, up, errvel (optional) whereas xyz are
+ 830 instrument coordinates.
+ 831 """
+ 832 ens.enu = rdi_xyz_enu(
+ 833 ens.xyze,
+ 834 ens.heading + self.magdec,
+ 835 ens.pitch,
+ 836 ens.roll,
+ 837 orientation=self.orientation,
+ 838 )
+ 839
+ 840 def _burst_average_depth(self, ens):
+ 841 """Depth-average within a burst.
+ 842
+ 843 Average the depth vectors if doing burst-averages. Otherwise we just
+ 844 return all depth vectors of this ensemble.
+ 845
+ 846 Parameters
+ 847 ----------
+ 848 ens : Bunch
+ 849 Ensemble data.
+ 850
+ 851 Returns
+ 852 -------
+ 853 depth : array-like
+ 854 Depth vector for each ping.
+ 855
+ 856 """
+ 857 if self.is_burst_average:
+ 858 depth_mean = ens.depth.mean(axis=0)
+ 859 depth = np.tile(depth_mean, (ens.dday.size, 1))
+ 860 else:
+ 861 depth = ens.depth
+ 862 return depth
+ 863
+ 864 def _regrid_enu(self, ens, method="linear"):
+ 865 """Depth-grid enu velocities."""
+ 866 shape = (ens.dday.size, self.dgrid.size, ens.enu.shape[-1])
+ 867 enu_grid = np.ma.zeros(shape)
+ 868 enu_grid[:] = np.ma.masked
+ 869 # Average the depth vectors if doing burst-averages. Otherwise we just
+ 870 # return all depth vectors of this ensemble.
+ 871 # TO DO: If calculating averages over regular time intervals, we need
+ 872 # to low pass filter the pressure time series prior to creating the
+ 873 # depth vectors.
+ 874 depth = self._burst_average_depth(ens)
+ 875 for i in range(ens.dday.size):
+ 876 enu_grid[i] = interp1(
+ 877 depth[i], ens.enu[i], self.dgrid, axis=0, method=method
+ 878 )
+ 879 ens.enu_grid = enu_grid
+ 880
+ 881 def _regrid_amp(self, ens, method="linear"):
+ 882 """Depth-grid amplitudes (averaged over all 4 beams)."""
+ 883 shape = (ens.dday.size, self.dgrid.size)
+ 884 amp_grid = np.ma.zeros(shape)
+ 885 amp_grid[:] = np.ma.masked
+ 886 depth = self._burst_average_depth(ens)
+ 887 for i in range(ens.dday.size):
+ 888 amp_grid[i] = interp1(
+ 889 depth[i],
+ 890 ens.amp[i].mean(axis=-1),
+ 891 self.dgrid,
+ 892 axis=0,
+ 893 method=method,
+ 894 )
+ 895 ens.amp_grid = amp_grid
+ 896
+ 897 def _binmap_one_beam(self, ens, beam_number):
+ 898 """Binmap single ping data for a single beam by linear interpolation.
+ 899
+ 900 Mapping is applied to velocity, amplitude and correlation.
+ 901
+ 902 Currently only works for 4 beam ADCP.
+ 903
+ 904 Parameters
+ 905 ----------
+ 906 ens : Bunch
+ 907 An ADCP dataset read by Multiread.
+ 908 beam_number : int
+ 909 An integer from 1 to 4 representing the beam number.
+ 910
+ 911 Returns
+ 912 -------
+ 913 veli : ndarray
+ 914 Mapped velocity.
+ 915 ampi : ndarray
+ 916 Mapped amplitude.
+ 917 cori : ndarray
+ 918 Mapped correlation.
+ 919
+ 920 """
+ 921
+ 922 if beam_number not in [1, 2, 3, 4]:
+ 923 raise ValueError("Beam number must be 1, 2, 3 or 4.")
+ 924
+ 925 tba = np.tan(np.deg2rad(ens.sysconfig.angle)) # Tangent of beam angle
+ 926 pitch = np.deg2rad(ens.pitch)
+ 927 roll = np.deg2rad(ens.roll)
+ 928
+ 929 # The true bin distances
+ 930 if beam_number == 1:
+ 931 dep = (
+ 932 ens.dep[None, :]
+ 933 * ((np.cos(roll) - tba * np.sin(roll)) * np.cos(pitch))[:, None]
+ 934 ) # None adds a new axis.
+ 935 elif beam_number == 2:
+ 936 dep = (
+ 937 ens.dep[None, :]
+ 938 * ((np.cos(roll) + tba * np.sin(roll)) * np.cos(pitch))[:, None]
+ 939 )
+ 940 elif beam_number == 3:
+ 941 dep = (
+ 942 ens.dep[None, :]
+ 943 * ((np.cos(pitch) + tba * np.sin(pitch)) * np.cos(roll))[:, None]
+ 944 )
+ 945 elif beam_number == 4:
+ 946 dep = (
+ 947 ens.dep[None, :]
+ 948 * ((np.cos(pitch) - tba * np.sin(pitch)) * np.cos(roll))[:, None]
+ 949 )
+ 950
+ 951 vel = ens.vel[..., beam_number - 1]
+ 952 amp = ens.amp[..., beam_number - 1]
+ 953 cor = ens.cor[..., beam_number - 1]
+ 954
+ 955 # Calculate interpolating weights (this hogs RAM!)
+ 956 dz = np.diff(dep, axis=1)
+ 957 w = np.clip((ens.dep - dep[:, :-1, None]) / dz[:, :, None], 0, 1)
+ 958
+ 959 # Determine data above or below the deepest bins
+ 960 above = (w == 1.0).all(axis=1)
+ 961 below = (w == 0.0).all(axis=1)
+ 962 bad = above | below
+ 963
+ 964 # Calculate differences
+ 965 dvel = np.diff(vel, axis=1)
+ 966 damp = np.diff(amp, axis=1)
+ 967 dcor = np.diff(cor, axis=1)
+ 968
+ 969 veli = vel[:, [0]] + np.sum(w * dvel[:, :, None], axis=1)
+ 970 veli[bad] = np.nan
+ 971
+ 972 ampi = amp[:, [0]] + np.sum(w * damp[:, :, None], axis=1)
+ 973 ampi[bad] = np.nan
+ 974
+ 975 cori = cor[:, [0]] + np.sum(w * dcor[:, :, None], axis=1)
+ 976 cori[bad] = np.nan
+ 977
+ 978 return veli, ampi, cori
+ 979
+ 980 def _binmap_all_beams(self, ens):
+ 981 """Binmap single ping data for all beams."""
+ 982
+ 983 for beam_number in [1, 2, 3, 4]:
+ 984 veli, ampi, cori = self._binmap_one_beam(ens, beam_number)
+ 985 ens.vel[..., beam_number - 1] = veli
+ 986 ens.amp[..., beam_number - 1] = ampi
+ 987 ens.cor[..., beam_number - 1] = cori
+ 988
+ 989 ens[f"vel{beam_number}"] = veli
+ 990 ens[f"amp{beam_number}"] = ampi
+ 991 ens[f"cor{beam_number}"] = cori
+ 992
+ 993 def _calculate_xyze(self, ens, ibad=None):
+ 994 """Calculate xyze from along-beam data."""
+ 995 if ens.sysconfig.convex:
+ 996 geom = "convex"
+ 997 else:
+ 998 geom = "concave"
+ 999
+1000 trans = Transform(angle=ens.sysconfig.angle, geometry=geom)
+1001
+1002 ens.xyze = trans.beam_to_xyz(ens.vel, ibad=ibad)
+1003
+1004 def process_pings(self, start=None, stop=None, binmap=False, ens_size=50000):
+1005 """Process single ping data without averaging.
+1006
+1007 Adds results as dictionary under `ave` and as `xarray.Dataset` under `ds`.
+1008
+1009 Writes processing parameters to the log file.
+1010
+1011 Parameters
+1012 ----------
+1013 start : int, optional
+1014 Start processing at this ping number.
+1015 stop : int, optional
+1016 Stop processing at this ping number.
+1017 binmap : bool, optional
+1018 Do binmapping of along-beam data.
+1019 ens_size : int, optional
+1020 Pings are processed in ensembles to reduce memory usage.
+1021 This parameter sets how many pings are in an ensemble. The default is 50000.
+1022
+1023 """
+1024 if start is None:
+1025 idx_start = np.searchsorted(self.dday, self.dday_start)
+1026 if stop is None:
+1027 idx_stop = np.searchsorted(self.dday, self.dday_end)
+1028
+1029 ens_idxs = np.hstack((np.arange(idx_start, idx_stop, ens_size), idx_stop))
+1030 write_idxs = ens_idxs - ens_idxs[0] # Arrays we write to start at index 0
+1031 npings = idx_stop - idx_start
+1032 nens = ens_idxs.size - 1
+1033 ndgrid = self.tsdat.dep.size
+1034
+1035 logger.info("Processing all pings")
+1036 logger.info(f"Binmapping is {binmap}")
+1037
+1038 uvwe = np.ma.zeros((npings, ndgrid, 4), dtype=np.float32)
+1039
+1040 pg = np.zeros((npings, ndgrid), dtype=np.int8)
+1041 amp = np.ma.zeros((npings, ndgrid), dtype=np.float32)
+1042
+1043 # temperature = np.ma.zeros((npings,), dtype=np.float32)
+1044 # pressure = np.ma.zeros((npings,), dtype=np.float32)
+1045
+1046 temperature = self.tsdat.temperature[idx_start:idx_stop]
+1047 pressure = self.tsdat.pressure[idx_start:idx_stop]
+1048
+1049 dday = self.dday[idx_start:idx_stop]
+1050
+1051 # Loop over ensembles
+1052 for i in tqdm(range(nens)):
+1053 ens = self.m.read(start=ens_idxs[i], stop=ens_idxs[i + 1])
+1054 idx0 = write_idxs[i]
+1055 idx1 = write_idxs[i + 1]
+1056
+1057 if ens is not None:
+1058 # I pulled this out of read_ensembles because we need to calculate depth for _ave2nc to work.
+1059 ens.dday_orig = ens.dday
+1060 ens.dday = self._correct_dday(ens.dday_orig)
+1061 if self._pressure_provided is not None:
+1062 # Replace pressure if provided
+1063 ens.pressure = self._external_pressure_to_dat(ens)
+1064 else:
+1065 ens.pressure = self._scale_pycurrents_pressure(ens)
+1066 sign = -1 if self.orientation == "up" else 1
+1067 pdepth = seawater.depth2(ens.pressure, self.lat)
+1068 ens.depth = pdepth[:, np.newaxis] + sign * ens.dep
+1069
+1070 if binmap:
+1071 self._binmap_all_beams(ens)
+1072 # Now we have to recalculate xyze with the binmapped data.
+1073 self._calculate_xyze(ens)
+1074
+1075 self._edit(ens) # modifies xyze
+1076 self._to_enu(ens) # transform to earth coords (east, north, up)
+1077
+1078 else:
+1079 uvwe[idx0:idx1] = np.ma.masked
+1080 amp[idx0:idx1] = np.ma.masked
+1081 # pressure[idx0:idx1] = np.ma.masked
+1082 # temperature[idx0:idx1] = np.ma.masked
+1083 continue
+1084
+1085 uvwe[idx0:idx1] = ens.enu
+1086
+1087 # pgi = 100 * ens.enu_grid[..., 0].count(axis=0) // nprofs
+1088 # pg[i] = pgi.astype(np.int8)
+1089 amp[idx0:idx1] = ens.amp.mean(axis=-1) # Average over beams... why?
+1090
+1091 self.ave = Bunch(
+1092 u=uvwe[..., 0],
+1093 v=uvwe[..., 1],
+1094 w=uvwe[..., 2],
+1095 e=uvwe[..., 3],
+1096 pg=pg,
+1097 amp=amp,
+1098 temperature=temperature,
+1099 pressure=pressure,
+1100 # npings=npings,
+1101 dday=dday,
+1102 yearbase=self.yearbase,
+1103 dep=ens.dep, # <<<<----- BAD!!! This a fudge because I don't want to calculated a depth vector. We use the depth vector from the last ensemble.
+1104 editparams=self.editparams,
+1105 tgridparams=self.tgridparams,
+1106 # dgridparams=self.dgridparams,
+1107 magdec=self.magdec,
+1108 lon=self.lon,
+1109 lat=self.lon,
+1110 )
+1111
+1112 self._ave2nc()
+1113 self._add_meta_data_to_ds()
+1114 self._log_processing_params()
+1115
+1116 def average_ensembles(self, start=None, stop=None):
+1117 """Time-averaging.
+1118
+1119 Adds results as dictionary under `ave` and as `xarray.Dataset` under `ds`.
+1120
+1121 Writes processing parameters to the log file.
+1122
+1123 Parameters
+1124 ----------
+1125 start : int
+1126 Range start for averaging. Index into start times of averaging
+1127 intervals.
+1128 stop : int
+1129 Range start for averaging. Index into start times of averaging
+1130 intervals.
+1131
+1132 """
+1133
+1134 nens_orig = len(self.start_ddays)
+1135 indices_orig = np.arange(nens_orig)
+1136 indices = indices_orig[start:stop]
+1137 if start is None and stop is None:
+1138 logger.info("Averaging all ensembles")
+1139 else:
+1140 logger.info(f"Averaging ensembles {indices[0]} to {indices[-1]}")
+1141 nens = len(indices)
+1142 ndgrid = len(self.dgrid)
+1143 uvwe = np.ma.zeros((nens, ndgrid, 4), dtype=np.float32)
+1144 uvwe_std = np.ma.zeros((nens, ndgrid, 4), dtype=np.float32)
+1145
+1146 pg = np.zeros((nens, ndgrid), dtype=np.int8)
+1147 amp = np.ma.zeros((nens, ndgrid), dtype=np.float32)
+1148
+1149 temperature = np.ma.zeros((nens,), dtype=np.float32)
+1150 pressure = np.ma.zeros((nens,), dtype=np.float32)
+1151 pressure_std = np.ma.zeros((nens,), dtype=np.float32)
+1152 pressure_max = np.ma.zeros((nens,), dtype=np.float32)
+1153
+1154 npings = np.zeros((nens,), dtype=np.int16)
+1155
+1156 dday = self.dday_mid[start:stop]
+1157
+1158 for i, iens in enumerate(tqdm(indices)):
+1159 ens = self.read_ensemble(iens)
+1160 if ens is not None:
+1161 self._edit(ens) # modifies xyze
+1162 self._to_enu(ens) # transform to earth coords (east, north, up)
+1163 self._regrid_enu(ens)
+1164 self._regrid_amp(ens)
+1165
+1166 nprofs = ens.enu_grid.shape[0]
+1167 else:
+1168 nprofs = 0
+1169 npings[i] = nprofs
+1170 if nprofs < 2:
+1171 uvwe[i] = np.ma.masked
+1172 uvwe_std[i] = np.ma.masked
+1173 # (pg is not a masked array)
+1174 amp[i] = np.ma.masked
+1175 pressure[i] = np.ma.masked
+1176 pressure_std[i] = np.ma.masked
+1177 pressure_max[i] = np.ma.masked
+1178 temperature[i] = np.ma.masked
+1179 continue
+1180
+1181 uvwe[i] = ens.enu_grid.mean(axis=0)
+1182 uvwe_std[i] = ens.enu_grid.std(axis=0)
+1183
+1184 pgi = 100 * ens.enu_grid[..., 0].count(axis=0) // nprofs
+1185 pg[i] = pgi.astype(np.int8)
+1186 amp[i] = ens.amp_grid.mean(axis=0)
+1187
+1188 pressure[i] = ens.pressure.mean()
+1189 pressure_std[i] = ens.pressure.std()
+1190 pressure_max[i] = ens.pressure.max()
+1191 temperature[i] = ens.temperature.mean()
+1192
+1193 self.ave = Bunch(
+1194 u=uvwe[..., 0],
+1195 v=uvwe[..., 1],
+1196 w=uvwe[..., 2],
+1197 e=uvwe[..., 3],
+1198 u_std=uvwe_std[..., 0],
+1199 v_std=uvwe_std[..., 1],
+1200 w_std=uvwe_std[..., 2],
+1201 e_std=uvwe_std[..., 3],
+1202 pg=pg,
+1203 amp=amp,
+1204 temperature=temperature,
+1205 pressure=pressure,
+1206 pressure_std=pressure_std,
+1207 pressure_max=pressure_max,
+1208 npings=npings,
+1209 dday=dday,
+1210 yearbase=self.yearbase,
+1211 dep=self.dgrid,
+1212 editparams=self.editparams,
+1213 tgridparams=self.tgridparams,
+1214 dgridparams=self.dgridparams,
+1215 magdec=self.magdec,
+1216 lon=self.lon,
+1217 lat=self.lon,
+1218 )
+1219
+1220 self._ave2nc()
+1221 self._add_meta_data_to_ds()
+1222 self._log_processing_params()
+1223
+1224 def burst_average_ensembles(self, start=None, stop=None, interpolate_bin=None):
+1225 """Time-averaging prior to depth-gridding.
+1226
+1227 Uses pre-defined editing parameters that can be updated with
+1228 `parse_editparams`.
+1229
+1230 Adds results as dictionary under `ave` and as `xarray.Dataset` under `ds`.
+1231
+1232 Writes processing parameters to the log file.
+1233
+1234 Parameters
+1235 ----------
+1236 start : int, optional
+1237 Range start for averaging. Index into start times of averaging
+1238 intervals. Defaults to None (start at beginning).
+1239 stop : int, optional
+1240 Range start for averaging. Index into start times of averaging
+1241 intervals. Defaults to None (start at beginning).
+1242 interpolate_bin : int or None, optional
+1243 Interpolate over a single, previously masked, bin. Defaults to None (no interpolation).
+1244
+1245 """
+1246 pg_condition = self.editparams.pg_limit
+1247 nens_orig = len(self.start_ddays)
+1248 indices_orig = np.arange(nens_orig)
+1249 indices = indices_orig[start:stop]
+1250 if start is None and stop is None:
+1251 logger.info("Averaging all ensembles")
+1252 else:
+1253 logger.info(f"Averaging ensembles {indices[0]} to {indices[-1]}")
+1254 nens = len(indices)
+1255 ndgrid = len(self.dgrid)
+1256 uvwe = np.ma.zeros((nens, ndgrid, 4), dtype=np.float32)
+1257 uvwe_std = np.ma.zeros((nens, ndgrid, 4), dtype=np.float32)
+1258
+1259 pg = np.zeros((nens, ndgrid), dtype=np.int8)
+1260 amp = np.ma.zeros((nens, ndgrid), dtype=np.float32)
+1261
+1262 temperature = np.ma.zeros((nens,), dtype=np.float32)
+1263 pressure = np.ma.zeros((nens,), dtype=np.float32)
+1264 pressure_std = np.ma.zeros((nens,), dtype=np.float32)
+1265 pressure_max = np.ma.zeros((nens,), dtype=np.float32)
+1266
+1267 npings = np.zeros((nens,), dtype=np.int16)
+1268
+1269 dday = self.dday_mid[start:stop]
+1270
+1271 for i, iens in enumerate(tqdm(indices)):
+1272 ens = self.read_ensemble(iens)
+1273 if ens is not None:
+1274 self._edit(ens) # modifies xyze
+1275 self._to_enu(ens) # transform to earth coords (east, north, up)
+1276 self._regrid_enu(ens)
+1277 self._regrid_amp(ens)
+1278
+1279 nprofs = ens.enu_grid.shape[0]
+1280 else:
+1281 nprofs = 0
+1282 npings[i] = nprofs
+1283 if nprofs < 2:
+1284 uvwe[i] = np.ma.masked
+1285 uvwe_std[i] = np.ma.masked
+1286 # (pg is not a masked array)
+1287 amp[i] = np.ma.masked
+1288 pressure[i] = np.ma.masked
+1289 pressure_std[i] = np.ma.masked
+1290 pressure_max[i] = np.ma.masked
+1291 temperature[i] = np.ma.masked
+1292 continue
+1293
+1294 # Depth vector for interpolation
+1295 depth = self._burst_average_depth(ens)
+1296 depth = depth[0, :]
+1297
+1298 # Calculate burst-average on instrument-relative depth grid
+1299 uvwe_inst = ens.enu.mean(axis=0)
+1300 uvwe_std_inst = ens.enu.std(axis=0)
+1301
+1302 pgi_inst = 100 * ens.enu[..., 0].count(axis=0) // nprofs
+1303
+1304 if pg_condition is not None:
+1305 pgi_index = pgi_inst < pg_condition
+1306 uvwe_inst[pgi_index, :] = np.ma.masked
+1307 uvwe_std_inst[pgi_index, :] = np.ma.masked
+1308
+1309 if interpolate_bin is not None:
+1310 zi = interpolate_bin
+1311 neighbors = [zi - 2, zi - 1, zi + 1, zi + 2]
+1312 dtmp = depth[neighbors]
+1313 tmp = interp1(
+1314 dtmp,
+1315 uvwe_inst[neighbors, :],
+1316 depth[zi],
+1317 axis=0,
+1318 method="linear",
+1319 )
+1320 uvwe_inst[zi, :] = tmp
+1321
+1322 # Interpolate burst-average to universal depth grid.
+1323 uvwe_grid = interp1(depth, uvwe_inst, self.dgrid, axis=0, method="linear")
+1324 uvwe_std_grid = interp1(
+1325 depth, uvwe_std_inst, self.dgrid, axis=0, method="linear"
+1326 )
+1327 uvwe[i] = uvwe_grid
+1328 uvwe_std[i] = uvwe_std_grid
+1329
+1330 # Interpolate pg to universal depth grid. Not overly satisfying but
+1331 # seems like that's what we need to do here.
+1332 pgi_grid = interp1(depth, pgi_inst, self.dgrid, axis=0, method="linear")
+1333 pg[i] = pgi_grid.astype(np.int8)
+1334
+1335 # Not changed to averaging in instrument-relative coordinates first.
+1336 amp[i] = ens.amp_grid.mean(axis=0)
+1337
+1338 pressure[i] = ens.pressure.mean()
+1339 pressure_std[i] = ens.pressure.std()
+1340 pressure_max[i] = ens.pressure.max()
+1341 temperature[i] = ens.temperature.mean()
+1342
+1343 self.ave = Bunch(
+1344 u=uvwe[..., 0],
+1345 v=uvwe[..., 1],
+1346 w=uvwe[..., 2],
+1347 e=uvwe[..., 3],
+1348 u_std=uvwe_std[..., 0],
+1349 v_std=uvwe_std[..., 1],
+1350 w_std=uvwe_std[..., 2],
+1351 e_std=uvwe_std[..., 3],
+1352 pg=pg,
+1353 amp=amp,
+1354 temperature=temperature,
+1355 pressure=pressure,
+1356 pressure_std=pressure_std,
+1357 pressure_max=pressure_max,
+1358 npings=npings,
+1359 dday=dday,
+1360 yearbase=self.yearbase,
+1361 dep=self.dgrid,
+1362 editparams=self.editparams,
+1363 tgridparams=self.tgridparams,
+1364 dgridparams=self.dgridparams,
+1365 magdec=self.magdec,
+1366 lon=self.lon,
+1367 lat=self.lon,
+1368 )
+1369
+1370 self._ave2nc()
+1371 self._add_meta_data_to_ds()
+1372 self._log_processing_params()
+1373
+1374 def _safely_add_attribute_from_params(self, key, d):
+1375 """Add attribute from dictionary and throw an exception if the key is missing.
+1376
+1377 Parameters
+1378 ----------
+1379 key : str
+1380 d : dict
+1381 """
+1382 try:
+1383 setattr(self, key, d[key])
+1384 except KeyError as error:
+1385 print(f"{error} is missing in input parameters")
+1386 raise
+1387
+1388 def _set_up_logger(self):
+1389 """Set up logging to both a file and to screen.
+1390
+1391 Default for screen logging is to show only warnings unless `verbose` is
+1392 set to True.
+1393
+1394 """
+1395
+1396 # Parse metadata for naming the log.
+1397 if self.meta_data is not None:
+1398 if "sn" in self.meta_data and "mooring" in self.meta_data:
+1399 sn = self.meta_data["sn"]
+1400 mooring = self.meta_data["mooring"]
+1401 has_mooring_id = True
+1402 else:
+1403 sn = None
+1404 mooring = None
+1405 has_mooring_id = False
+1406 else:
+1407 has_mooring_id = False
+1408
+1409 # Set up a logger that will collect log info from this module and the
+1410 # other pycurrent methods as well.
+1411 logdir = Path(self.logdir)
+1412 logdir.mkdir(exist_ok=True)
+1413 if has_mooring_id:
+1414 filename = f"{mooring}_{sn}.log"
+1415 else:
+1416 filename = "adcp_proc.log"
+1417 # Delete any existing handlers. This may be bad style, but I kept adding handlers when developing this.
+1418 logger.handlers = []
+1419 logging.basicConfig(
+1420 filename=logdir.joinpath(filename),
+1421 filemode="w",
+1422 format="%(asctime)s %(name)s %(levelname)s %(message)s",
+1423 datefmt="%H:%M:%S",
+1424 level=logging.INFO,
+1425 force=True,
+1426 )
+1427 ConsoleOutputHandler = logging.StreamHandler()
+1428 ConsoleOutputHandler.setLevel(logging.WARNING)
+1429 if self.verbose:
+1430 ConsoleOutputHandler.setLevel(logging.INFO)
+1431 logger.addHandler(ConsoleOutputHandler)
+1432
+1433 # current date
+1434 datestr = np.datetime64(datetime.datetime.now()).astype(datetime.datetime)
+1435 strformat = "%Y-%m-%d"
+1436 datestr = datestr.strftime(strformat)
+1437
+1438 if has_mooring_id:
+1439 logger.info(f"Processing {mooring} SN {sn} on {datestr}")
+1440 else:
+1441 logger.warning(
+1442 "No meta data provided, logging to generic filename 'adcp_proc.log'"
+1443 )
+1444 logger.info(f"Processing on {datestr}")
+1445
+1446 def _log_processing_params(self):
+1447 """Write processing parameters to log file."""
+1448
+1449 logger.info("processing settings")
+1450 logger.info("-------------------")
+1451 self._log_params(self.dgridparams)
+1452 self._log_params(self.tgridparams)
+1453 self._log_params(self.editparams)
+1454 logger.info("-------------------")
+1455
+1456 def _log_params(self, pd):
+1457 """Print ADCP processing parameter dict to logs.
+1458
+1459 Parameters
+1460 ----------
+1461 pd : dict
+1462 Parameter dict.
+1463 """
+1464
+1465 for k, v in pd.items():
+1466 if k == "maskbins" and v is not None:
+1467 logstr = (k, ":", np.flatnonzero(v))
+1468 logger.info(logstr)
+1469 else:
+1470 logstr = (k, ":", v)
+1471 logger.info(logstr)
+1472
+1473 def _ave2nc(self):
+1474 """Convert data structure from ave to xarray.Dataset format."""
+1475 # load npz file
+1476 dat = self.ave.copy()
+1477 dat = Bunch(dat)
+1478 # identify variables
+1479 k = dat.keys()
+1480 varsint = []
+1481 vars1d = []
+1482 vars2d = []
+1483 for ki in k:
+1484 try:
+1485 tmp = dat[ki].shape
+1486 if len(tmp) == 1:
+1487 vars1d.append(ki)
+1488 elif len(tmp) == 2:
+1489 vars2d.append(ki)
+1490 except:
+1491 tmp = None
+1492 varsint.append(ki)
+1493 # print(ki, tmp)
+1494
+1495 # generate time vector
+1496 base = datetime.datetime(dat.yearbase, 1, 1, 0, 0, 0)
+1497 time = [base + datetime.timedelta(days=ti) for ti in dat.dday]
+1498 adcptime = [np.datetime64(ti, "ns") for ti in time]
+1499 # generate Dataset
+1500 out = xr.Dataset(
+1501 {"pg": (["depth", "time"], dat.pg.T)},
+1502 coords={"time": (["time"], adcptime), "depth": (["depth"], dat.dep)},
+1503 )
+1504 for vari in vars2d:
+1505 out[vari] = (["depth", "time"], dat[vari].T)
+1506 for vari in vars1d:
+1507 if vari not in ["dep", "dday"]:
+1508 out[vari] = (["time"], dat[vari])
+1509
+1510 # Percent good is currently defined everywhere. Set to NaN where we
+1511 # don't have any amplitude data (i.e. no data).
+1512 out["pg"] = out.pg.where(~np.isnan(out.amp), other=np.nan)
+1513
+1514 # Drop depth levels with all nan
+1515 out = out.dropna(how="all", dim="depth")
+1516
+1517 # Drop pressure_std, pressure_max, and e_std
+1518 dropvars = ["pressure_std", "pressure_max", "e_std"]
+1519 for var in dropvars:
+1520 if var in out:
+1521 out = out.drop(var)
+1522
+1523 # Change u/v/w std to standard error by dividing by sqrt(npings)
+1524 for var in ["u", "v", "w"]:
+1525 if f"{var}_std" in out:
+1526 out = out.rename({f"{var}_std": f"{var}_error"})
+1527 out[f"{var}_error"] = out[f"{var}_error"] / np.sqrt(out["npings"])
+1528
+1529 # Calculate transducer depth from pressure
+1530 out["xducer_depth"] = -gsw.z_from_p(out.pressure, self.lat)
+1531
+1532 # add variable names and units for plotting
+1533 out = self._add_names_and_units(out)
+1534
+1535 self.ds = out
+1536
+1537 def _add_names_and_units(self, ds):
+1538 """Add variable attributes based on CF conventions.
+1539
+1540 Parameters
+1541 ----------
+1542 ds : xarray.Dataset
+1543
+1544 Returns
+1545 -------
+1546 ds : xarray.Dataset
+1547
+1548 """
+1549 CF = io.cf_conventions()
+1550 for v in ds.variables:
+1551 if v in CF:
+1552 ds[v].attrs = CF[v]
+1553 return ds
+1554
+1555 def _add_meta_data_to_ds(self):
+1556 # Add some more info.
+1557 self.ds.attrs["orientation"] = self.orientation
+1558 self.ds.attrs["magdec"] = self.magdec
+1559 for att in ["max_e", "max_e_deviation", "min_correlation"]:
+1560 self.ds.attrs[att] = self.editparams[att]
+1561
+1562 # Add meta data if provided.
+1563 if self.meta_data is not None:
+1564 for k, v in self.meta_data.items():
+1565 self.ds.attrs[k] = v
+1566 self.ds.attrs["proc time"] = np.datetime64("now").astype("str")
+1567
+1568 def plot_echo_stats(self):
+1569 """Plot beam statistics (correlation and amplitude) from raw ADCP data."""
+1570 r = self.raw
+1571
+1572 fig, ax = plt.subplots(
+1573 nrows=1,
+1574 ncols=2,
+1575 figsize=(5, r.bin.max().data * 0.3),
+1576 constrained_layout=True,
+1577 sharey=True,
+1578 )
+1579 r.cor.mean(dim="time").plot(
+1580 hue="beam", y="bin", marker="o", linestyle="", ax=ax[0]
+1581 )
+1582 r.amp.mean(dim="time").plot(
+1583 hue="beam", y="bin", marker="o", linestyle="", ax=ax[1]
+1584 )
+1585 ax[0].invert_yaxis()
+1586 ax[1].set(ylabel="")
+1587 ax[0].set_yticks(r.bin.data)
+1588 for axi in ax:
+1589 axi.grid(True)
+1590
+1591 def plot_pressure(self):
+1592 """Plot pressure time series and mark time at depth."""
+1593 fig, ax = plt.subplots(
+1594 nrows=1,
+1595 ncols=1,
+1596 figsize=(6, 2.5),
+1597 constrained_layout=True,
+1598 )
+1599 self.raw.pressure.plot(ax=ax, label="all")
+1600 self.raw.pressure.where(self.raw.pressure > 50).plot(ax=ax, label="subsurface")
+1601 ax.invert_yaxis()
+1602 ax.set(xlabel="", ylabel="pressure [dbar]")
+1603 ax.legend()
+1604
+1605 def generate_binmask(self, indices):
+1606 binmask = self.raw.bin.data < 0
+1607 binmask[indices] = True
+1608 return binmask
+1609
+1610
+1611class ProcessADCPyml(ProcessADCP):
+1612 """Moored ADCP processing with parameters provided via .yml file.
+1613
+1614 An example parameter file is included at
+1615 [`notebooks/parameters.yml`](https://github.com/modscripps/velosearaptor/tree/main/notebooks/parameters.yml)
+1616 """
+1617
+1618 def __init__(self, parameter_file, mooring, sn, **kwargs):
+1619 p = io.parse_yaml_input(parameter_file, mooring, sn)
+1620 super().__init__(
+1621 p["data_dir"],
+1622 meta_data=p["meta_data"],
+1623 driftparams=p["driftparams"],
+1624 tgridparams=p["tgridparams"],
+1625 dgridparams=p["dgridparams"],
+1626 editparams=p["editparams"],
+1627 **kwargs,
+1628 )
+