Add automatic datatime column as needed when running raiderCombine.py

CHANGELOG.md

@@ -10,6 +10,9 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 * [811](https://github.com/dbekaert/RAiDER/pull/811) - Fixed a north-south flip in the ECMWF/ERA-5 model-level reader. `_load_model_level` reversed only the level axis of the geopotential cube (`z[::-1]`) while `t`, `q`, and `lnsp` had their latitude axis reversed, leaving heights mirrored in latitude relative to the meteorology. Surface height and surface pressure were anti-correlated (r = -0.54 where physics requires ~+1); sea-level ZTD over a 3-degree box spanned 1621-3233 mm instead of 2239-2395 mm. The error is antisymmetric in latitude, so it largely cancels in a domain average while being severe at individual points.
 
 ### Changed
+* [809](https://github.com/dbekaert/RAiDER/pull/809) - Update `raiderCombine.py` to add a `Datetime` column on-the-fly. Previously running `raiderDownloadGNSS.py` and then `raiderCombine.py` directly after could fail. `readZTDFile` now builds the column from a `Date` column (optionally with a `times` seconds-of-day column), an existing `Datetime` column, or a `YYYYMMDDTHHMMSS` timestamp in the filename, in that order, and raises an explanatory error naming the columns it did find when none of those are available. The filename fallback stamps every row with a single timestamp and now logs a warning saying so.
+  * `Lat`/`Lon`/`Hgt_m` are backfilled into the GNSS frame from the RAiDER delay file, keyed on station ID, for GNSS ZTD files that carry no station coordinates.
+  * The workflow now aborts with an explanatory error when the RAiDER and GNSS files share no common observations, rather than passing empty frames into the per-station variance analysis and failing later with a `KeyError` on a column that was never created.
 * [808](https://github.com/dbekaert/RAiDER/pull/808) - Use `xarray.open_dataset` instead of `load_dataset` for memory efficiency, and guard weather-model file reads with explicit error handling so that datasets are always closed. Dropped the redundant full-cube `shutil.copy` from the ECMWF model-level download, and added regression tests covering the layout contract between `ECMWF._get_from_cds` and `ECMWF._makeDataCubes`.
 
 ## [0.6.0]

test/test_gnss.py

@@ -6,7 +6,13 @@
 import pytest
 
 from RAiDER.gnss.downloadGNSSDelays import download_tropo_delays, filterToBBox, get_station_list, get_stats_by_llh
-from RAiDER.gnss.processDelayFiles import addDateTimeToFiles, concatDelayFiles, getDateTime
+from RAiDER.gnss.processDelayFiles import (
+    addDateTimeToFiles,
+    concatDelayFiles,
+    getDateTime,
+    readZTDFile,
+)
+from RAiDER.gnss.processDelayFiles import main as process_main
 from RAiDER.models.customExceptions import NoStationDataFoundError
 from test import TEST_DIR, pushd
 
@@ -150,3 +156,169 @@ def test_filterByBBox2():
         assert stat not in new_data['ID'].to_list()
     for stat in ['FGNW', 'JPLT', 'NVTP', 'WLHG', 'WORG']:
         assert stat in new_data['ID'].to_list()
+
+
+def write_ztd_csv(path: Path, **columns) -> Path:
+    """Write a GNSS ZTD file with whatever columns a test needs."""
+    pd.DataFrame(columns).to_csv(path, index=False)
+    return path
+
+
+def test_readZTDFile_date_and_times(tmp_path):
+    """'Date' plus a 'times' seconds-of-day column combine into 'Datetime'."""
+    f = write_ztd_csv(
+        tmp_path / 'UNRcombinedGPS_ztd.csv',
+        ID=['STAT1', 'STAT2'],
+        Date=['2020-01-30', '2020-01-30'],
+        times=[0, 45296],  # midnight and 12:34:56
+        ZTD=[2.3, 2.4],
+    )
+    data = readZTDFile(f)
+    assert data['Datetime'].to_list() == [
+        datetime.datetime(2020, 1, 30, 0, 0, 0),
+        datetime.datetime(2020, 1, 30, 12, 34, 56),
+    ]
+
+
+def test_readZTDFile_date_without_times(tmp_path):
+    """Without a 'times' column every observation lands at midnight."""
+    f = write_ztd_csv(
+        tmp_path / 'UNRcombinedGPS_ztd.csv',
+        ID=['STAT1', 'STAT2'],
+        Date=['2020-01-30', '2020-01-31'],
+        ZTD=[2.3, 2.4],
+    )
+    data = readZTDFile(f)
+    assert data['Datetime'].to_list() == [
+        datetime.datetime(2020, 1, 30),
+        datetime.datetime(2020, 1, 31),
+    ]
+
+
+def test_readZTDFile_unparseable_times_are_treated_as_midnight(tmp_path):
+    """A junk 'times' value is coerced to zero rather than failing the read."""
+    f = write_ztd_csv(
+        tmp_path / 'UNRcombinedGPS_ztd.csv',
+        ID=['STAT1', 'STAT2'],
+        Date=['2020-01-30', '2020-01-30'],
+        times=['not-a-number', 3600],
+        ZTD=[2.3, 2.4],
+    )
+    data = readZTDFile(f)
+    assert data['Datetime'].to_list() == [
+        datetime.datetime(2020, 1, 30, 0, 0, 0),
+        datetime.datetime(2020, 1, 30, 1, 0, 0),
+    ]
+
+
+def test_readZTDFile_existing_datetime_column(tmp_path):
+    """An existing 'Datetime' column is used as-is when there is no 'Date'."""
+    f = write_ztd_csv(
+        tmp_path / 'UNRcombinedGPS_ztd.csv',
+        ID=['STAT1'],
+        Datetime=['2020-01-30 13:52:45'],
+        ZTD=[2.3],
+    )
+    data = readZTDFile(f)
+    assert data['Datetime'].to_list() == [datetime.datetime(2020, 1, 30, 13, 52, 45)]
+
+
+def test_readZTDFile_datetime_from_filename(tmp_path, caplog):
+    """With no time column at all, the filename timestamp is used, with a warning."""
+    f = write_ztd_csv(
+        tmp_path / 'ERA5_Delay_20210308T120000_ztd.csv',
+        ID=['STAT1', 'STAT2'],
+        ZTD=[2.3, 2.4],
+    )
+    data = readZTDFile(f)
+
+    assert data['Datetime'].to_list() == [datetime.datetime(2021, 3, 8, 12, 0, 0)] * 2
+    # Every row gets the same stamp, so the user has to be told.
+    assert 'no longer distinguishable in time' in caplog.text
+
+
+def test_readZTDFile_raises_without_any_time_information(tmp_path):
+    """No time column and no timestamp in the filename is an error, not a guess."""
+    f = write_ztd_csv(
+        tmp_path / 'stations.csv',
+        ID=['STAT1'],
+        Lat=[15.0],
+        Lon=[-100.0],
+        ZTD=[2.3],
+    )
+    with pytest.raises(ValueError, match='cannot be added automatically'):
+        readZTDFile(f)
+
+
+def test_readZTDFile_renames_the_delay_column(tmp_path):
+    """The delay column named by col_name is renamed to 'ZTD'."""
+    f = write_ztd_csv(
+        tmp_path / 'UNRcombinedGPS_ztd.csv',
+        ID=['STAT1'],
+        Date=['2020-01-30'],
+        totalDelay=[2.3],
+    )
+    data = readZTDFile(f, col_name='totalDelay')
+    assert 'totalDelay' not in data.columns
+    assert data['ZTD'].to_list() == [2.3]
+
+
+def write_delay_pair(directory: Path, ztd_days, raider_days, with_coords=True):
+    """Write a matched RAiDER/GNSS delay file pair for the combine workflow.
+
+    The GNSS file optionally omits Lat/Lon/Hgt_m, as UNR per-station delay
+    files do, so that the backfill from the RAiDER file can be exercised.
+    """
+    stations = ['STA1', 'STA2']
+    raider_rows, ztd_rows = [], []
+    for station in stations:
+        for k, day in enumerate(raider_days):
+            raider_rows.append({
+                'ID': station, 'Lat': 34.0, 'Lon': -118.0, 'Hgt_m': 100.0,
+                'Datetime': day + pd.Timedelta(hours=12),
+                'wetDelay': 0.1, 'hydroDelay': 2.2, 'totalDelay': 2.3 + 0.01 * k,
+            })
+        for k, day in enumerate(ztd_days):
+            # Residuals alternate in sign so the per-station variance analysis
+            # has something to work with rather than a constant offset.
+            ztd_rows.append({
+                'ID': station, 'Datetime': day + pd.Timedelta(hours=12),
+                'ZTD': 2.3 + 0.01 * k + (0.02 if k % 2 else -0.02),
+                'sigZTD': 0.005, 'times': 43200,
+            })
+
+    ztd_frame = pd.DataFrame(ztd_rows)
+    if with_coords:
+        ztd_frame['Lat'] = 34.0
+        ztd_frame['Lon'] = -118.0
+        ztd_frame['Hgt_m'] = 100.0
+
+    raider_file = directory / 'raider_delays.csv'
+    ztd_file = directory / 'gnss_delays.csv'
+    pd.DataFrame(raider_rows).to_csv(raider_file, index=False)
+    ztd_frame.to_csv(ztd_file, index=False)
+    return raider_file, ztd_file
+
+
+def test_main_backfills_station_coordinates_from_raider_file(tmp_path):
+    """A GNSS file with no station coordinates inherits them from the RAiDER file."""
+    days = pd.date_range('2020-01-01', periods=12, freq='D')
+    raider_file, ztd_file = write_delay_pair(tmp_path, days, days, with_coords=False)
+    assert 'Lat' not in pd.read_csv(ztd_file).columns
+
+    out_path = tmp_path / 'combined.csv'
+    process_main(raider_file, ztd_file, out_path=out_path)
+
+    combined = pd.read_csv(out_path)
+    assert combined['Lat'].dropna().unique().tolist() == [34.0]
+    assert combined['Lon'].dropna().unique().tolist() == [-118.0]
+
+
+def test_main_raises_when_files_share_no_observations(tmp_path):
+    """Non-overlapping dates fail with an explanation, not a later KeyError."""
+    raider_days = pd.date_range('2020-01-01', periods=12, freq='D')
+    ztd_days = pd.date_range('2021-06-01', periods=12, freq='D')
+    raider_file, ztd_file = write_delay_pair(tmp_path, ztd_days, raider_days)
+
+    with pytest.raises(ValueError, match='No common observations'):
+        process_main(raider_file, ztd_file, out_path=tmp_path / 'combined.csv')

tools/RAiDER/gnss/processDelayFiles.py

@@ -240,28 +240,69 @@ def local_time_filter(raiderFile, ztdFile, dfr, dfz, localTime):
 
 
 def readZTDFile(filename, col_name='ZTD'):
-    """Read and parse a GPS zenith delay file."""
-    try:
-        data = pd.read_csv(filename, parse_dates=['Date'])
+    """Read and parse a GPS zenith delay file.
+
+    A 'Datetime' column is constructed from whatever time information the
+    file provides, in this order of preference:
+      1. a 'Date' column (optionally combined with a 'times' column of
+         seconds-of-day),
+      2. an existing 'Datetime' column,
+      3. a YYYYMMDDTHHMMSS timestamp parsed from the filename.
+    If none of these are available, a clear error is raised instead of a
+    cryptic pandas parsing error.
+    """
+    data = pd.read_csv(filename)
+
+    if 'Date' in data.columns:
         date0 = pd.to_datetime(data['Date'],
             errors='raise',
             format='%Y-%m-%d')
 
-        # If present, convert seconds → pandas Timedelta; otherwise zero
+        # If present, convert seconds → pandas Timedelta; otherwise zero.
+        # sec has to stay a Series either way: pd.to_timedelta on a bare 0
+        # returns a scalar Timedelta, which has no .values for the sum below.
         if 'times' in data.columns:
             sec = pd.to_numeric(data['times'], errors='coerce').fillna(0)
-            td = pd.to_timedelta(sec, unit='s')
         else:
-            td = pd.to_timedelta(0, unit='s')
+            sec = pd.Series(0, index=data.index)
+        td = pd.to_timedelta(sec, unit='s')
 
         # Combine using numpy/pandas arrays
         # (stays in datetime64[ns], never Python objects)
-        dt_vals = date0.values + td.values 
+        dt_vals = date0.values + td.values
 
         # Assign back
         data['Datetime'] = pd.to_datetime(dt_vals)
-    except (KeyError, ValueError):
-        data = pd.read_csv(filename, parse_dates=['Datetime'])
+    elif 'Datetime' in data.columns:
+        data['Datetime'] = pd.to_datetime(data['Datetime'], errors='raise')
+    else:
+        # Neither a 'Date' nor a 'Datetime' column is present. Try to add a
+        # 'Datetime' column automatically from a timestamp in the filename
+        # (e.g. ..._20210308T000000.csv); if that is not possible, fail with
+        # an informative message rather than a cryptic pandas error.
+        try:
+            file_datetime = getDateTime(Path(filename))
+        except (AttributeError, ValueError):
+            raise ValueError(
+                f"File '{filename}' has no 'Date' or 'Datetime' column and no "
+                "parseable YYYYMMDDTHHMMSS timestamp in its filename, so a "
+                "'Datetime' column cannot be added automatically. Columns "
+                f"found: {list(data.columns)}. If you passed a directory, "
+                "check that it points at the GNSS/model delay files and not, "
+                "e.g., a station-list file."
+            )
+        # The filename carries one timestamp, so every row in the file gets it.
+        # That is correct for a single-epoch delay file and wrong for anything
+        # holding several epochs, so say so rather than silently collapsing them.
+        logger.warning(
+            'File %s has no "Date" or "Datetime" column, so all %d rows are being stamped '
+            'with %s, parsed from the filename. If this file covers more than one epoch, '
+            'those observations are no longer distinguishable in time.',
+            filename,
+            len(data),
+            file_datetime,
+        )
+        data['Datetime'] = file_datetime
 
     data.rename(columns={col_name: 'ZTD'}, inplace=True)
     return data
@@ -724,6 +765,16 @@ def main(
     expected_data_columns = ['ID', 'Lat', 'Lon', 'Hgt_m', 'Datetime', 'wetDelay', 'hydroDelay', raider_delay]
     dfr = dfr.drop(columns=[col for col in dfr if col not in expected_data_columns])
 
+    # Some GNSS ZTD files (e.g. UNR per-station delay files) carry only
+    # ID/Date/ZTD/... and no station coordinates. Backfill Lat/Lon/Hgt_m from
+    # the RAiDER delay file, which is keyed on the same station IDs, so the
+    # GNSS frame is self-consistent for the lat/lon mapping and local-time
+    # estimation below. IDs absent from the RAiDER file get NaNs and are
+    # dropped downstream (they could not be matched anyway).
+    for _coord in ('Lat', 'Lon', 'Hgt_m'):
+        if _coord not in dfz.columns and _coord in dfr.columns:
+            dfz[_coord] = dfz['ID'].map(dict(zip(dfr['ID'], dfr[_coord])))
+
     # Create dictionaries mapping ID → Lat and ID → Lon from GNSS file
     lat_map = dict(zip(dfz["ID"], dfz["Lat"]))
     lon_map = dict(zip(dfz["ID"], dfz["Lon"]))
@@ -758,6 +809,20 @@ def main(
     dfz = pass_common_obs(dfr, dfz)
     dfr = pass_common_obs(dfz, dfr)
 
+    # Bail out early with an informative message if the GNSS and RAiDER files
+    # share no common observations. Otherwise the empty frames propagate to the
+    # per-station variance analysis and surface as a cryptic KeyError on a
+    # column that was never created (e.g. 'sigma_model_neg').
+    if dfz.empty or dfr.empty:
+        raise ValueError(
+            f"No common observations between RAiDER file '{raider_file}' and "
+            f"GNSS file '{ztd_file}'. The two datasets do not overlap in "
+            "station ID and/or calendar date (matching is by date and ID). "
+            "Check that the RAiDER delays were computed for the same dates as "
+            "the GNSS observations; a systematic offset of even one day will "
+            "produce zero matches."
+        )
+
     # If specified, convert to local-time reference frame WRT 0 longitude
     common_keys = ['Datetime', 'ID']
     if local_time is not None: