[FEAT][TEST] Adding relevant variables for flight model

seagliderOG1/config/OG1_var_names.yaml

@@ -86,5 +86,18 @@ eng_vbdCC: VBD_CC
 buoyancy: BUOYANCY
 sound_velocity: SOUND_VELOCITY
 
+### flight model parameters ###
+sg_cal_vbd_min_cnts: VBD_MIN_CNTS
+log_VBD_MIN: VBD_MIN_CNTS
+sg_cal_vbd_cnts_per_cc: VBD_CNTS_PER_CC
+log_VBD_CNV: VBD_CC_PER_CNTS
+sg_cal_mass: MASS
+sg_cal_volmax: VOLMAX
+log_C_VBD: C_VBD
+sg_cal_hd_a: HD_A
+sg_cal_hd_b: HD_B
+sg_cal_hd_c: HD_C
+sg_cal_vbdbias: VBD_BIAS
+
 sbe41: Seabird unpumped CTD
 wlbb2f: Wetlabs BB2FL-VMT

seagliderOG1/convertOG1.py

@@ -72,6 +72,10 @@ def convert_to_OG1(list_of_datasets: list[xr.Dataset] | xr.Dataset, contrib_to_a
     for key, value in ordered_attributes.items():
         ds_og1.attrs[key] = value
 
+    ### Add information needed/used for hydrodynamic (flight) model (hdm)
+    hdm_parameters = tools.extract_hdm_parameters(list_of_datasets)
+    ds_og1 = tools.add_hdm_parameters(ds_og1, hdm_parameters)
+
     # Construct the platform serial number
     if "platform_id" in ds1_base.attrs:
         PLATFORM_SERIAL_NUMBER = ds1_base.attrs["platform_id"].lower()

seagliderOG1/tools.py

@@ -1122,6 +1122,122 @@ def combine_two_dim_of_dataset(
 
     return updated_ds
 
+standard_names = vocabularies.standard_names
+
+def extract_hdm_parameters(list_datasets):
+    """
+    Extracts HDM parameters and their attributes from a list of datasets. If the parameter has the same value across all datasets,
+    it keeps a single value; otherwise, it returns the full list.
+
+    Parameters:
+    -----------
+    list_datasets (list): List of xarray.Dataset objects.
+    standard_names (dict): Vocabulary mapping internal names to standard names.
+
+    Returns:
+    --------
+    dict: A nested dictionary where keys are standard names and values contain
+            the 'data' and 'attributes'.
+    """
+    potential_parameters_OG1 = ['VBD_MIN_CNTS','VBD_CNTS_PER_CC','VBD_CC_PER_CNTS','VBD_BIAS','MASS','VOLMAX','C_VBD','HD_A','HD_B','HD_C']
+    potential_parameters = [key for key, value in standard_names.items() if value in potential_parameters_OG1]
+    hdm_variables = {}
+
+    for param in potential_parameters:
+        # Determine the key name using standard_names mapping
+        param_key = standard_names.get(param, param)
+        if param_key == param and param not in standard_names:
+            print(f"Warning: '{param}' not found in standard names. Using original name.")
+
+        # 1. Check if the parameter exists in the datasets
+        if param not in list_datasets[0].variables:
+            continue
+
+        # 2. Extract values from all datasets
+        all_values = [ds[param].values for ds in list_datasets]
+
+        # 3. Determine if we keep a single value or the full list
+        # We flatten to handle case where .values might be arrays
+        unique_vals = np.unique(np.array(all_values))
+
+        final_value = unique_vals[0] if len(unique_vals) == 1 else all_values
+
+        if isinstance(final_value, list):
+            final_value = np.array(final_value)
+
+        # 4. Store as a dictionary to accommodate both value and attributes and add long_name attribute
+        hdm_variables[param_key] = {
+            "values": final_value,
+            "attributes": list_datasets[0][param].attrs  # Add long_name attribute
+        }
+        hdm_variables[param_key]["attributes"]["original_name"] = param_key
+
+    # 5. Print what parameters from potential_parameters were found and which couldn't not be found in the datasets
+    found_params = [param for param in potential_parameters_OG1 if param in hdm_variables]
+    not_found_params = [param for param in potential_parameters_OG1 if param not in hdm_variables]
+    print(f"The following HDM parameters were found: {found_params}")
+    if not_found_params:
+        print(f"Warning: The following potential HDM parameters were not found in the datasets: {not_found_params}")
+
+    return hdm_variables
+
+def add_hdm_parameters(ds_OG1, hdm_parameters):
+    """
+    Add HDM parameters to the OG1 dataset as new variables with their attributes.
+
+    Parameters:
+    -----------
+    ds_OG1 (xarray.Dataset): The OG1 dataset to which HDM parameters will be added.
+    hdm_parameters (dict): A dictionary containing HDM parameters and their attributes
+                            in the format {standard_name: {"value": ..., "attributes": {...}}}.
+    Returns:
+    --------
+    xarray.Dataset: Updated OG1 dataset with HDM parameters added as variables.
+    """
+    ds_updated = ds_OG1.copy()
+
+    # Get unique profile numbers (e.g., [1, 2, 3, 4...])
+    unique_profiles = np.sort(np.unique(ds_updated.PROFILE_NUMBER.values))
+    num_profiles = len(unique_profiles)
+
+    for param_name, param_info in hdm_parameters.items():
+        # Using .get() because you used "value" in extract and "values" in your draft
+        values = param_info.get("value") or param_info.get("values")
+        attributes = param_info["attributes"]
+
+        # Check if it's a single value (scalar)
+        if np.size(values) == 1:
+            # item() converts a 1-element array to a native python scalar
+            ds_updated[param_name] = values.item() if hasattr(values, "item") else values
+            ds_updated[param_name].attrs = attributes
+
+        # Check if it's dive-based (1 value per 2 profiles)
+        elif np.size(values) == num_profiles // 2:
+            # Initialize an array of zeros/NaNs matching the measurement dimension
+            # Use the same dtype as our input values
+            mapped_array = np.full(ds_updated.N_MEASUREMENTS.shape, np.nan)
+
+            # Iterate through dives (each dive = 2 profiles)
+            for i, dive_val in enumerate(values):
+                # Create mask for the two profiles corresponding to one dive
+                p_idx_1, p_idx_2 = unique_profiles[2 * i], unique_profiles[2 * i + 1]
+
+                # Find all measurement indices belonging to these two profiles
+                mask = (ds_updated.PROFILE_NUMBER == p_idx_1) | (ds_updated.PROFILE_NUMBER == p_idx_2)
+
+                # Fill the array for those specific measurements
+                mapped_array[mask] = dive_val
+
+            # Add to dataset with the N_MEASUREMENTS dimension
+            ds_updated[param_name] = (("N_MEASUREMENTS",), mapped_array)
+            ds_updated[param_name].attrs = attributes
+
+        else:
+            print(f"Warning: {param_name} size ({np.size(values)}) does not match "
+                  f"scalar or dive count ({num_profiles // 2}). Skipping.")
+
+    return ds_updated
+
 
 # ===============================================================================
 # Unused functions

tests/test_tools.py

@@ -8,7 +8,7 @@
 import numpy as np
 import xarray as xr
 import gsw
-from seagliderOG1 import tools
+from seagliderOG1 import tools, readers, convertOG1
 
 
 def test_convert_units_var():
@@ -53,3 +53,16 @@ def test_calc_z():
     depth_z = tools.calc_Z(dataset)["DEPTH_Z"].values
 
     assert np.array_equal(depth, depth_z)
+
+def test_add_hdm_parameters():
+    source = str(parent_dir / "data/demo_sg005")
+    print("Testing load_basestation_files with source:", source)
+    start_profile = 1
+    end_profile = 5
+    datasets = readers.load_basestation_files(source, start_profile, end_profile)
+    hdm_parameters = tools.extract_hdm_parameters(datasets)
+    ds_OG1, vars = convertOG1.convert_to_OG1(datasets)
+    ds_OG1 = tools.add_hdm_parameters(ds_OG1, hdm_parameters)
+    ### check if the hdm parameters are added to the dataset and have the expected values
+    for param in hdm_parameters:
+        assert param in ds_OG1