"""User-facing utilities for Seaglider-to-OG1 conversion.
Sensor metadata extraction, variable-name mapping, QC-flag conversion, dtype
coercion, profile and phase numbering, time encoding, and dataset merging.
"""
import logging
import gsw
import numpy as np
import pandas as pd
import xarray as xr
from seagliderOG1 import vocabularies
_log = logging.getLogger(__name__)
__all__ = [
"add_depth",
"add_phase",
"add_profile_number",
"add_dive_number",
"convert_qc_flags",
"encode_times_og1",
"split_by_unique_dims",
]
# Variables measured directly by the CTD.
# Variables calculated using CTD measurements.
# Alternative names used for the same physical instrument.
[docs]
def add_dive_number(ds: xr.Dataset, dive_number: int | None = None) -> xr.Dataset:
"""Add dive number as a variable to the dataset. Assumes present in the basestation attributes.
Parameters
----------
ds
The dataset to which the dive number will be added.
dive_number, optional
The dive number to add. If None, extracts from dataset attributes.
Returns
-------
xarray.Dataset
The dataset with the dive number added.
"""
if dive_number is None:
dive_number = ds.attrs.get("dive_number", np.nan)
dive_var = xr.DataArray(
np.full(ds.sizes["N_MEASUREMENTS"], dive_number),
dims=["N_MEASUREMENTS"],
attrs=dict(vocabularies.vocab_attrs["DIVE_NUMBER"]),
)
return ds.assign(DIVE_NUMBER=dive_var)
[docs]
def add_profile_number(ds: xr.Dataset, ds1: xr.Dataset) -> xr.Dataset:
"""Assign profile numbers to measurements based on dive phases (down/up casts).
This function separates each dive into two profiles: descent (down cast) and
ascent (up cast) phases. The dive is split at the maximum pressure point,
with the descent phase getting the dive number and ascent getting dive + 0.5.
Profile numbers are then calculated as 2 * DIVE_NUMBER - 1.
Parameters
----------
ds
The dataset to add profile numbers to.
ds1
Dataset containing dive number information in attributes.
Returns
-------
xarray.Dataset
Dataset with 'DIVE_NUMBER' and 'PROFILE_NUMBER' variables added.
Notes
-----
- Requires pressure variable (PRES)
"""
dive_number = ds1.attrs["dive_number"]
ds = add_dive_number(ds, dive_number)
fill_value = -9999
profile_numbers = np.full(ds.sizes["N_MEASUREMENTS"], fill_value, dtype=int)
if profile_numbers.size and not np.isnan(dive_number):
pmax_index = int(np.nanargmax(ds["PRES"].values))
profile_numbers[: pmax_index + 1] = 2 * dive_number - 1
profile_numbers[pmax_index + 1 :] = 2 * dive_number
ds = ds.assign(
PROFILE_NUMBER=xr.DataArray(
profile_numbers,
dims=["N_MEASUREMENTS"],
# Built here, not mapped from a basestation variable, so the vocabulary
# attributes are applied here (copied: vocab_attrs is shared).
attrs=dict(vocabularies.vocab_attrs["PROFILE_NUMBER"]),
)
)
ds["PROFILE_NUMBER"].encoding["_FillValue"] = fill_value
return ds
[docs]
def add_phase(ds: xr.Dataset) -> xr.Dataset:
"""Add 'PHASE' and 'PHASE_QC' variables indicating the phase of each measurement.
The phase is determined from the pressure readings ('PRES') for each unique
dive number ('dive_num').
Parameters
----------
ds : xarray.Dataset
The input dataset containing 'dive_num' and 'PRES' variables.
Returns
-------
xarray.Dataset
The dataset with additional 'PHASE' and 'PHASE_QC' variables, where:
- 'PHASE' indicates the phase of each measurement:
- Phase 2 up to and including the maximum pressure point.
- Phase 1 after the maximum pressure point.
- 'PHASE_QC' is an additional variable with no QC applied.
Notes
-----
In this formulation, we are only separating into dives and climbs based on
when the glider is at the maximum depth. Future work needs to separate out
the other phases and generate a PHASE_QC:
https://github.com/OceanGlidersCommunity/OG-format-user-manual/blob/main/vocabularyCollection/phase.md
"""
# Determine the correct keystring for divenum
if "dive_number" in ds.variables:
divenum_str = "dive_number"
elif "divenum" in ds.variables:
divenum_str = "divenum"
elif "dive_num" in ds.variables:
divenum_str = "dive_num"
elif "DIVE_NUMBER" in ds.variables:
divenum_str = "DIVE_NUMBER"
else:
raise ValueError("No valid dive number variable found in the dataset.")
# Initialize the new variable with the same dimensions as dive_num
ds["PHASE"] = (
["N_MEASUREMENTS"],
np.full(ds.sizes["N_MEASUREMENTS"], np.nan),
dict(vocabularies.vocab_attrs["PHASE"]),
)
# Initialize the new variable PHASE_QC with the same dimensions as dive_num
ds["PHASE_QC"] = (
["N_MEASUREMENTS"],
np.zeros(ds.sizes["N_MEASUREMENTS"], dtype=int),
)
# Iterate over each unique dive_num
for dive in np.unique(ds[divenum_str]):
# Get the indices for the current dive
dive_indices = np.where(ds[divenum_str] == dive)[0]
# Find the start and end index for the current dive
start_index = dive_indices[0]
end_index = dive_indices[-1]
# Find the index of the maximum pressure between start_index and end_index
pmax = np.nanmax(ds["PRES"][start_index : end_index + 1].values)
# Find the index where PRES attains the value pmax between start_index and end_index
pmax_index = start_index + np.argmax(
ds["PRES"][start_index : end_index + 1].values == pmax
)
# Assign phase 2 to all values up to and including the point where pmax is reached
ds["PHASE"][start_index : pmax_index + 1] = 2
# Assign phase 1 to all values after pmax is reached
ds["PHASE"][pmax_index + 1 : end_index + 1] = 1
# Assign phase 3 to the time at the beginning of the dive, between the first valid TIME_GPS and the second valid TIME_GPS
valid_time_gps_indices = np.where(
~np.isnan(ds["TIME_GPS"][start_index : end_index + 1].values)
)[0]
if len(valid_time_gps_indices) >= 2:
first_valid_index = start_index + valid_time_gps_indices[0]
second_valid_index = start_index + valid_time_gps_indices[1]
ds["PHASE"][first_valid_index : second_valid_index + 1] = 3
return ds
##-----------------------------------------------------------------------------------------------------------
## Calculations for new variables
##-----------------------------------------------------------------------------------------------------------
[docs]
def add_depth(ds: xr.Dataset) -> xr.Dataset:
"""Calculate the depth (Z position) of the glider using the gsw library to convert pressure to depth.
Parameters
----------
ds
The input dataset containing 'PRES', 'LATITUDE', and 'LONGITUDE' variables.
Returns
-------
xarray.Dataset
The dataset with an additional 'DEPTH' variable.
"""
# Ensure the required variables are present
if "PRES" not in ds.variables or "LATITUDE" not in ds.variables:
raise ValueError("Dataset must contain 'PRES' and 'LATITUDE' variables.")
# Convert pressure to depth using gsw (pressure in dbar, latitude in degrees)
depth = gsw.z_from_p(
ds["PRES"], ds["LATITUDE"]
).compute() # Compute to handle dask arrays
# Add depth to dataset
ds["DEPTH_Z"] = (["N_MEASUREMENTS"], depth.data)
# Assign the calculated depth to a new variable in the dataset
ds["DEPTH_Z"].attrs = {
"units": "meters",
"positive": "up",
"standard_name": "depth",
"comment": "Depth calculated from pressure using gsw library, positive up.",
}
return ds
[docs]
def split_by_unique_dims(ds: xr.Dataset) -> dict:
"""Split an xarray dataset into datasets grouped by each variable's dimension set.
Parameters
----------
ds
The input xarray dataset containing various variables.
Returns
-------
dict
A dictionary mapping dimension tuples to datasets, each with variables sharing the same set of dimensions.
"""
# Group the variables first and build each dataset once: adding variables one at a
# time re-aligns the growing dataset on every assignment (quadratic in the number
# of variables, about a second per basestation dive).
groups: dict[tuple, dict[str, xr.DataArray]] = {}
for var_name, var_data in ds.data_vars.items():
groups.setdefault(tuple(var_data.sizes), {})[var_name] = var_data
return {dims: xr.Dataset(variables) for dims, variables in groups.items()}
[docs]
def convert_qc_flags(dsa: xr.Dataset, qc_name: str) -> xr.Dataset:
"""Convert QC flag variables to proper integer format and update attributes.
This function converts QC flag variables from string format to int8,
handles NaN values appropriately, removes ``QC_`` prefixes from flag meanings,
and adds proper metadata including long_name and standard_name.
Parameters
----------
dsa
The dataset containing QC flag variables.
qc_name
The name of the QC flag variable to process.
Returns
-------
xarray.Dataset
Dataset with converted QC flag variable and updated attributes.
Notes
-----
Must be called after the main variable has been assigned its OG1 long_name.
"""
# Must be called *after* var_name has OG1 long_name
var_name = qc_name[:-3]
if qc_name in list(dsa):
# For scalar log qc's that might change over the mission are stored as json, don't convert to int8
if (
qc_name in dsa.variables
and dsa[qc_name].attrs.get("serialization") == "json"
):
return dsa
# Seaglider default type was a string. Convert to int8 and take care of NaNs
# dsa[qc_name].values = dsa[qc_name].values.astype("int8")
values = dsa[qc_name].values
# Convert byte strings to regular strings (if necessary)
if values.dtype.type is np.bytes_:
values = values.astype(str)
# Use pandas to handle NaNs safely
original_shape = np.asarray(values).shape
values = pd.to_numeric(
np.asarray(values).reshape(-1),
errors="coerce",
).reshape(original_shape) # Convert strings to numbers, NaNs stay NaNs
# Assign back to dataset
dsa[qc_name].values = values
### Set the nan values to 6 (unsampled flag) and convert to int8
### Before it had just set all values to 0, which is no change flag
### Alternative could be to set to 9 (missing value)
dsa[qc_name].values = dsa[qc_name].fillna(6).astype("int8")
# A flag variable has no missing value (6 = unsampled is itself a flag), so drop
# any inherited float _FillValue that would be invalid on the int8 result.
# Clear both encoding and attrs: set_best_dtype skips QC, so nothing else removes
# a _FillValue left in attrs, and xarray rejects _FillValue in attrs on write.
dsa[qc_name].encoding.pop("_FillValue", None)
dsa[qc_name].attrs.pop("_FillValue", None)
# Seaglider default flag_meanings were prefixed with 'QC_'. Remove this prefix.
if "flag_meaning" in dsa[qc_name].attrs:
flag_meaning = dsa[qc_name].attrs["flag_meaning"]
dsa[qc_name].attrs["flag_meaning"] = flag_meaning.replace("QC_", "")
# The parent variable may be absent: an orphan QC, or a numbered name like
# 'TEMP_QC2' whose 'TEMP_' parent does not exist (qc_name[:-3] strips only the
# '_QC' suffix). Guard before dereferencing dsa[var_name].
parent_present = var_name in dsa.variables
# Add a long_name attribute to the QC variable
parent_long_name = (
dsa[var_name].attrs.get("long_name", "") if parent_present else ""
)
dsa[qc_name].attrs["long_name"] = f"{parent_long_name} quality flag".strip()
dsa[qc_name].attrs["standard_name"] = "status_flag"
# Mention the QC variable in the parent's attributes, when the parent exists.
if parent_present:
dsa[var_name].attrs["ancillary_variables"] = qc_name
return dsa
# Variables that are semantically integer and are stored as the smallest signed integer
# that holds their range. Named explicitly: a float variable that happens to be
# integer-valued in one deployment is not an integer variable.
# OG1 canonical serialisation for time variables (CF seconds since the epoch, UTC).
# Used by the time encoders here and by writers.write; keep in sync with the time
# variable units declared in config/OG1_variable_attributes.yaml.
OG1_TIME_UNITS = "seconds since 1970-01-01T00:00:00Z"
OG1_TIME_CALENDAR = "gregorian"
[docs]
def encode_times_og1(ds: xr.Dataset) -> xr.Dataset:
"""Encode time variables according to OG1 format specifications.
Parameters
----------
ds
Dataset containing time variables to encode.
Returns
-------
xarray.Dataset
Dataset with OG1-formatted time variables.
"""
for var_name in ds.variables:
if "axis" in ds[var_name].attrs:
ds[var_name].attrs.pop("axis")
# Only real datetime variables are CF-time encoded: their units and calendar
# belong in encoding and are written there, so any stray attrs are dropped. A
# numeric duration such as DELTA_TIME_S merely has "time" in its name and must
# keep its unit attribute (seconds), so the name is not used to select here.
if np.issubdtype(ds[var_name].dtype, np.datetime64):
for drop_attr in ["units", "calendar", "dtype"]:
if drop_attr in ds[var_name].attrs:
ds[var_name].attrs.pop(drop_attr)
if drop_attr in ds[var_name].encoding:
ds[var_name].encoding.pop(drop_attr)
if var_name.lower() == "time":
ds[var_name].attrs["units"] = OG1_TIME_UNITS
ds[var_name].attrs["calendar"] = OG1_TIME_CALENDAR
return ds
# ===============================================================================
# Unused functions
# ===============================================================================
# --- deprecated v0.0.2 names (removed in 0.3.0; see seagliderOG1._deprecated) ----------
from seagliderOG1 import _deprecated # noqa: E402 # shims follow the definitions above
calc_Z = _deprecated._renamed(add_depth, "tools.calc_Z", "tools.add_depth")
assign_phase = _deprecated._renamed(add_phase, "tools.assign_phase", "tools.add_phase")
assign_profile_number = _deprecated._renamed(
add_profile_number, "tools.assign_profile_number", "tools.add_profile_number"
)
# The three unit helpers now resolve through the units view (vocabularies.UNITS); the ad-hoc
# unit tables they used are gone. Kept as warning shims for one release (removed 0.3.0).
reformat_units_str = _deprecated.reformat_units_str
reformat_units_var = _deprecated.reformat_units_var
convert_units_var = _deprecated.convert_units_var
_REMOVED_TOOLS = {
"OG1_name_mapping": (
"the basestation -> OG1 name mapping is built inside seagliderOG1.convert; "
"'seagliderOG1 inspect --conversion FILE.nc' prints what was renamed."
),
"print_OG1_mapping_summary": (
"the mapping is built inside seagliderOG1.convert; "
"'seagliderOG1 inspect --conversion FILE.nc' prints what was renamed."
),
}
_REMOVED_TOOLS.update(
dict.fromkeys(
(
"og1_name_mapping",
"gather_sensor_info",
"add_sensor_to_dataset",
"merge_datasets_along_time",
"extract_scalar_parameters",
"add_scalar_parameters",
"set_best_dtype",
),
"a conversion-pipeline internal, private since 0.1.0; call seagliderOG1.convert.",
)
)
_deprecated.install_removed(globals(), _REMOVED_TOOLS)