[FIX] Added functions to scan and repair NetCDF files before loading

seagliderOG1/convertOG1.py

@@ -18,7 +18,10 @@
 _log = logging.getLogger(__name__)
 
 
-def convert_to_OG1(list_of_datasets: list[xr.Dataset] | xr.Dataset, contrib_to_append: dict[str, str] | None = None,) -> tuple[xr.Dataset, list[str]]:
+def convert_to_OG1(
+    list_of_datasets: list[xr.Dataset] | xr.Dataset,
+    contrib_to_append: dict[str, str] | None = None,
+) -> tuple[xr.Dataset, list[str]]:
     """Convert Seaglider basestation datasets to OG1 format.
     Processes a list of xarray datasets or a single xarray dataset, converts them to OG1 format,
     concatenates the datasets, sorts by time, and applies attributes. Main conversion function that
@@ -64,6 +67,8 @@ def convert_to_OG1(list_of_datasets: list[xr.Dataset] | xr.Dataset, contrib_to_a
 
     ds_og1 = xr.concat(processed_datasets, dim="N_MEASUREMENTS")
     ds_og1 = ds_og1.sortby("TIME")
+    # Change format of time into datetime64[ns] to avoid problems with attributes and writing to netcdf
+    # ds_og1["TIME"] = (ds_og1["TIME"].astype("float64") * 1e9).astype("datetime64[ns]")
 
     # Apply attributes
     ordered_attributes = update_dataset_attributes(
@@ -104,12 +109,18 @@ def convert_to_OG1(list_of_datasets: list[xr.Dataset] | xr.Dataset, contrib_to_a
     ds_og1["TRAJECTORY"].attrs["long_name"] = "trajectory name"
     ds_og1["TRAJECTORY"].attrs["cf_role"] = "trajectory_id"
 
-    ds_og1["DEPLOYMENT_LATITUDE"] = xr.DataArray(ds_og1.LATITUDE.values[~np.isnan(ds_og1.LATITUDE)][0],
-                                              attrs = {"long_name": "latitude of deployment"})
-    ds_og1["DEPLOYMENT_LONGITUDE"] = xr.DataArray(ds_og1.LONGITUDE.values[~np.isnan(ds_og1.LONGITUDE)][0],
-                                               attrs = {"long_name": "longitude of deployment"})
-    ds_og1["DEPLOYMENT_TIME"] = xr.DataArray(ds_og1.TIME.values[~np.isnan(ds_og1.TIME)][0],
-                                             attrs = {"long_name": "time of deployment"})
+    ds_og1["DEPLOYMENT_LATITUDE"] = xr.DataArray(
+        ds_og1.LATITUDE.values[~np.isnan(ds_og1.LATITUDE)][0],
+        attrs={"long_name": "latitude of deployment"},
+    )
+    ds_og1["DEPLOYMENT_LONGITUDE"] = xr.DataArray(
+        ds_og1.LONGITUDE.values[~np.isnan(ds_og1.LONGITUDE)][0],
+        attrs={"long_name": "longitude of deployment"},
+    )
+    ds_og1["DEPLOYMENT_TIME"] = xr.DataArray(
+        ds_og1.TIME.values[~np.isnan(ds_og1.TIME)][0],
+        attrs={"long_name": "time of deployment"},
+    )
 
     # Remove attributes from TIME_GPS
     if "TIME_GPS" in ds_og1.variables:
@@ -274,8 +285,11 @@ def process_dataset(ds1_base: xr.Dataset, firstrun: bool = False) -> tuple[
     ds_new = tools.add_sensor_to_dataset(ds_new, ds_sensor, ds_sgcal, firstrun)
 
     # To avoid problems, reset the dtype of TIME_GPS
-    ds_new['TIME_GPS'] = ds_new['TIME_GPS'].astype('datetime64[ns]')
-    vars_to_remove = vocabularies.vars_to_remove #+ ["TIME_GPS"]
+    # ds_new['TIME_GPS'] = ds_new['TIME_GPS'].astype('datetime64[ns]')
+    ds_new["TIME_GPS"] = (ds_new["TIME_GPS"].astype("float64") * 1e9).astype(
+        "datetime64[ns]"
+    )
+    vars_to_remove = vocabularies.vars_to_remove  # + ["TIME_GPS"]
     vars_present_to_remove = [var for var in vars_to_remove if var in ds_new.variables]
 
     # Drop them
@@ -288,6 +302,7 @@ def process_dataset(ds1_base: xr.Dataset, firstrun: bool = False) -> tuple[
     attr_warnings = ""
     return ds_new, attr_warnings, ds_sgcal, ds_other, ds_log
 
+
 def standardise_OG10(
     ds: xr.Dataset,
     firstrun: bool = False,
@@ -367,7 +382,10 @@ def standardise_OG10(
             )
             ### Only log a warning for variables that aren't in the vocabularies and aren't in the list of variables to keep or remove
             ### Removed varaiables will be printed in the log as being removed, so no need to log a warning for them here.
-            if orig_varname not in (*vocabularies.vars_as_is, *vocabularies.vars_to_remove):
+            if orig_varname not in (
+                *vocabularies.vars_as_is,
+                *vocabularies.vars_to_remove,
+            ):
                 vars_not_in_vocab.append(orig_varname)
 
     if firstrun and vars_not_in_vocab:
@@ -505,7 +523,9 @@ def add_gps_info_to_dataset(ds: xr.Dataset, gps_ds: xr.Dataset) -> xr.Dataset:
 ##-----------------------------------------------------------------------------------------
 ## Editing attributes
 ##-----------------------------------------------------------------------------------------
-def update_dataset_attributes(ds: xr.Dataset, contrib_to_append: dict[str, str] | None) -> dict[str, str]:
+def update_dataset_attributes(
+    ds: xr.Dataset, contrib_to_append: dict[str, str] | None
+) -> dict[str, str]:
     """Update the attributes of the dataset based on the provided attribute input.
 
     Processes contributor information, time attributes, and applies OG1
@@ -570,7 +590,9 @@ def update_dataset_attributes(ds: xr.Dataset, contrib_to_append: dict[str, str]
     return ordered_attributes
 
 
-def get_contributors(ds: xr.Dataset, values_to_append: dict[str, str] | None = None) -> dict[str, str]:
+def get_contributors(
+    ds: xr.Dataset, values_to_append: dict[str, str] | None = None
+) -> dict[str, str]:
     """Extract and format contributor information for OG1 attributes.
 
     Processes creator and contributor information from dataset attributes,
@@ -589,6 +611,7 @@ def get_contributors(ds: xr.Dataset, values_to_append: dict[str, str] | None = N
         Dictionary with formatted contributor attribute strings.
 
     """
+
     # Function to create or append to a list
     def create_or_append_list(existing_list, new_item):
         if new_item not in existing_list:
@@ -808,7 +831,9 @@ def get_time_attributes(ds: xr.Dataset) -> dict[str, str]:
     return time_attrs
 
 
-def extract_attr_to_keep(ds1: xr.Dataset, attr_as_is: list[str] = vocabularies.global_attrs["attr_as_is"]) -> dict[str, str]:
+def extract_attr_to_keep(
+    ds1: xr.Dataset, attr_as_is: list[str] = vocabularies.global_attrs["attr_as_is"]
+) -> dict[str, str]:
     """Extract attributes to retain unchanged.
 
     Parameters
@@ -835,7 +860,8 @@ def extract_attr_to_keep(ds1: xr.Dataset, attr_as_is: list[str] = vocabularies.g
 
 
 def extract_attr_to_rename(
-    ds1: xr.Dataset, attr_to_rename: dict[str, str] = vocabularies.global_attrs["attr_to_rename"]
+    ds1: xr.Dataset,
+    attr_to_rename: dict[str, str] = vocabularies.global_attrs["attr_to_rename"],
 ) -> dict[str, str]:
     """Extract and rename attributes according to OG1 vocabulary.
 
@@ -861,7 +887,12 @@ def extract_attr_to_rename(
     return renamed_attrs
 
 
-def process_and_save_data(input_location: str, save: bool = False, output_dir: str = ".", run_quietly: bool = True) -> xr.Dataset:
+def process_and_save_data(
+    input_location: str,
+    save: bool = False,
+    output_dir: str = ".",
+    run_quietly: bool = True,
+) -> xr.Dataset:
     """Process and save data from the specified input location.
 
     This function loads and concatenates datasets from the server, converts them to OG1 format,

seagliderOG1/readers.py

@@ -3,12 +3,14 @@
 import re
 import sys
 
-import pooch
 import requests
 import xarray as xr
 from bs4 import BeautifulSoup
 from tqdm import tqdm
 
+from datetime import datetime
+import shutil
+
 script_dir = pathlib.Path(__file__).parent.absolute()
 parent_dir = script_dir.parents[0]
 sys.path.append(str(parent_dir))
@@ -51,10 +53,13 @@
     },
 )"""
 
+
 # Instead of loading from the server, we will load from a local directory for testing and development purposes.
 # The local directory will contain the same files as the server, but we will not use pooch to manage them.
 # Instead, we will just read them directly from the filesystem.
-def load_sample_dataset(file_path: str = str(parent_dir / "data/demo_sg005/p0050001_20080606.nc")) -> xr.Dataset:
+def load_sample_dataset(
+    file_path: str = str(parent_dir / "data/demo_sg005/p0050001_20080606.nc"),
+) -> xr.Dataset:
     """Download sample datasets for use with seagliderOG1.
 
     Parameters
@@ -80,6 +85,7 @@ def load_sample_dataset(file_path: str = str(parent_dir / "data/demo_sg005/p0050
         msg = f"Requested sample dataset {file_path} not known. Available datasets are: {os.listdir(str(parent_dir / 'data/demo_sg005'))}"
         raise KeyError(msg)
 
+
 def _validate_filename(filename: str) -> bool:
     """Validate if filename matches expected Seaglider basestation patterns.
 
@@ -155,7 +161,11 @@ def _glider_sn_from_filename(filename: str) -> int:
     return int(filename[1:4])
 
 
-def filter_files_by_profile(file_list: list[str], start_profile: int | None = None, end_profile: int | None = None) -> list[str]:
+def filter_files_by_profile(
+    file_list: list[str],
+    start_profile: int | None = None,
+    end_profile: int | None = None,
+) -> list[str]:
     """Filter files by profile/dive number range.
 
     Filters Seaglider basestation files based on profile number range.
@@ -230,7 +240,9 @@ def load_first_basestation_file(source: str) -> xr.Dataset:
     return datasets[0]
 
 
-def load_basestation_files(source: str, start_profile: int | None = None, end_profile: int | None = None) -> list[xr.Dataset]:
+def load_basestation_files(
+    source: str, start_profile: int | None = None, end_profile: int | None = None
+) -> list[xr.Dataset]:
     """Load multiple Seaglider basestation files with optional profile filtering.
 
     Main function for loading Seaglider data from either online repositories
@@ -251,22 +263,27 @@ def load_basestation_files(source: str, start_profile: int | None = None, end_pr
         List of loaded basestation datasets, ordered by filename.
 
     """
+    ### Scan all basestation files and repair any with inconsistent time metadata before loading
+    scan_and_repair_files(source, start_profile, end_profile)
+
     file_list = list_files(source)
     filtered_files = filter_files_by_profile(file_list, start_profile, end_profile)
 
     datasets = []
 
     ### Include a tqdm progress bar
     for file in tqdm(filtered_files, desc="Loading datasets", unit="file"):
-        ds = xr.open_dataset(os.path.join(source, file), decode_timedelta=False)
+        ds = xr.open_dataset(os.path.join(source, file))
 
         datasets.append(ds)
 
     return datasets
 
 
 def list_files(
-    source: str, registry_loc: str = "seagliderOG1", registry_name: str = "seaglider_registry.txt"
+    source: str,
+    registry_loc: str = "seagliderOG1",
+    registry_name: str = "seaglider_registry.txt",
 ) -> list[str]:
     """List NetCDF files from a source (URL or local directory).
 
@@ -315,3 +332,135 @@ def list_files(
     file_list.sort()
 
     return file_list
+
+
+def scan_and_repair_files(
+    source: str,
+    start_profile: int | None = None,
+    end_profile: int | None = None,
+) -> None:
+    """
+    Scan NetCDF files and repair inconsistent time metadata only when needed.
+
+    Parameters
+    ----------
+    source : str
+        Directory path containing NetCDF files to scan and repair.
+    start_profile : int, optional
+        Minimum profile number to include (inclusive).
+    end_profile : int, optional
+        Maximum profile number to include (inclusive).
+
+    Returns
+    -------
+    None
+    """
+    file_list = list_files(source)
+    filtered_files = filter_files_by_profile(file_list, start_profile, end_profile)
+
+    fixes_dir = os.path.join(source, "metadata_fixes")
+    os.makedirs(fixes_dir, exist_ok=True)
+
+    log_path = os.path.join(fixes_dir, "repair_log.txt")
+
+    for file in tqdm(filtered_files, desc="Scanning files", unit="file"):
+        full_path = os.path.join(source, file)
+
+        try:
+            # Good files are left untouched
+            ds = xr.open_dataset(full_path, decode_timedelta=False)
+            ds.close()
+
+        except Exception as err:
+            print(f"Need repair: {file}")
+
+            fixed_vars = repair_netcdf_time_metadata_inplace(
+                full_path,
+                repair_dir=fixes_dir,
+                backup=True,
+            )
+
+            if fixed_vars:
+                log_repair(log_path, file, fixed_vars, err)
+                print(f"Repaired {file}: {fixed_vars}")
+
+
+def _repair_folder(source: str | pathlib.Path) -> pathlib.Path:
+    """Return the folder where backups and logs for repaired files are stored."""
+    repair_dir = pathlib.Path(source) / "metadata_fixes"
+    repair_dir.mkdir(exist_ok=True)
+    return repair_dir
+
+
+def _backup_path(
+    path: str | pathlib.Path, repair_dir: str | pathlib.Path
+) -> pathlib.Path:
+    """Return a backup filename stored in metadata_fixes/."""
+    path = pathlib.Path(path)
+    repair_dir = pathlib.Path(repair_dir)
+    return repair_dir / f"{path.stem}_original{path.suffix}"
+
+
+def log_repair(
+    log_path: str | pathlib.Path,
+    filename: str,
+    fixed_vars: list[str],
+    error: Exception | str,
+):
+    """Append one repair entry to the log file."""
+    with open(log_path, "a") as f:
+        f.write(f"{datetime.now().isoformat()} | {filename}\n")
+        f.write(f"  Fixed variables: {fixed_vars}\n")
+        f.write(f"  Original error: {str(error)}\n")
+        f.write("\n")
+
+
+def repair_netcdf_time_metadata_inplace(
+    path: str | pathlib.Path,
+    repair_dir: str | pathlib.Path,
+    backup: bool = True,
+) -> list[str]:
+    """Repair variables stored as text but mislabeled as CF time variables.
+
+    The original file is backed up into repair_dir before the repaired version
+    replaces the original file.
+    """
+    path = pathlib.Path(path)
+    repair_dir = pathlib.Path(repair_dir)
+    fixed_vars: list[str] = []
+
+    ds = xr.open_dataset(path, decode_times=False, decode_timedelta=False)
+    try:
+        for var_name in ds.variables:
+            var = ds[var_name]
+            units = var.attrs.get("units")
+
+            if not isinstance(units, str):
+                continue
+
+            if "since" not in units.lower():
+                continue
+
+            if getattr(var.dtype, "kind", None) in {"S", "U", "O"}:
+                var.attrs.pop("units", None)
+                var.attrs.pop("calendar", None)
+                fixed_vars.append(var_name)
+
+        if not fixed_vars:
+            return []
+
+        ds.load()
+
+        if backup:
+            backup_path = _backup_path(path, repair_dir)
+            if not backup_path.exists():
+                shutil.copy2(path, backup_path)
+
+        tmp_path = path.with_name(f"{path.stem}_tmp_repaired{path.suffix}")
+        ds.to_netcdf(tmp_path)
+
+    finally:
+        ds.close()
+
+    os.replace(tmp_path, path)
+    return fixed_vars