diff --git a/.gitignore b/.gitignore index b9390d473..306cb0b21 100644 --- a/.gitignore +++ b/.gitignore @@ -35,3 +35,4 @@ orbits/ /GUNW*.yaml /*.tif .coverage* +TODO.md diff --git a/CHANGELOG.md b/CHANGELOG.md index 72741bf20..fbb1980f5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,12 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [PEP 440](https://www.python.org/dev/peps/pep-0440/) and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +### Changed +* [802](https://github.com/dbekaert/RAiDER/pull/802) - Update the CDS API to request all dates simultaneously when pulling model levels +- Use the `era5-pressure-levels` archive when accessing pressure level data. +- add unit tests to cover the changes +* [795](https://github.com/dbekaert/RAiDER/pull/795) - use xarray.open_dataset instead of load_dataset for memory efficiency. Also convert GNSS heights to geoid instead of native ellipsoid for increased accuracy when intersecting with the weather model cube + ## [0.6.0] ### Removed * [764](https://github.com/dbekaert/RAiDER/pull/764) - Removed Python 3.8 support. Python 3.9 is now the minimum version officially required to run RAiDER. diff --git a/pytest.ini b/pytest.ini index 42b90946c..256698e61 100644 --- a/pytest.ini +++ b/pytest.ini @@ -1,3 +1,4 @@ [pytest] markers = long: mark a test as a long + network: mark a test as requiring external network access (e.g. PROJ CDN for EGM96 grids) diff --git a/test/test_ecmwf_era5.py b/test/test_ecmwf_era5.py new file mode 100644 index 000000000..4279f9d61 --- /dev/null +++ b/test/test_ecmwf_era5.py @@ -0,0 +1,562 @@ +"""Unit tests for ECMWF/ERA5/ERA5T weather model loaders. + +These tests exercise the loading and axis-ordering logic that was broken +in the multiday_download branch (ecmwf.py commit 3b64ad3..d734fd7). The +primary failure mode was _load_pressure_level dividing geopotential heights +by g0 a second time and applying axis flips to the wrong dimensions, +producing hydrostatic ZTDs of ~1 m instead of ~2.3 m. +""" + +from __future__ import annotations + +import datetime as dt +from pathlib import Path +from unittest.mock import MagicMock, patch + +import numpy as np +import pytest +import xarray as xr + +from RAiDER.models.era5 import ERA5 +from RAiDER.models.era5t import ERA5T + + +# --------------------------------------------------------------------------- +# Helpers to build synthetic netCDF files +# --------------------------------------------------------------------------- + +G0 = 9.80665 # m/s² + +# Small but physically plausible grid +_NLAT = 4 +_NLON = 5 +_NLEV = 6 + +# Lats: N→S descending, as CDS returns them; loader must flip to S→N. +_LATS_DESC = np.array([52.0, 51.0, 50.0, 49.0]) +_LATS_ASC = _LATS_DESC[::-1] + +# Lons: ascending, centred on Europe +_LONS = np.array([10.0, 11.0, 12.0, 13.0, 14.0]) + +# Pressure levels in hPa, surface-first (1000 → 10), as the batch writer +# produces them from a standard CDS download. +_LEVS_HPA_DESC = np.array([1000.0, 500.0, 200.0, 100.0, 50.0, 10.0]) +_LEVS_PA_DESC = _LEVS_HPA_DESC * 100.0 + +# Geopotential heights (m) matching the pressure levels above, surface first. +# Values are broadly consistent with the standard atmosphere. +_Z_GHT_SURFACE_FIRST = np.array([100.0, 5_800.0, 11_800.0, 16_200.0, 20_600.0, 31_100.0]) + + +def _make_t(nlat=_NLAT, nlon=_NLON, nlev=_NLEV) -> np.ndarray: + """Return a (nlat, nlon, nlev) temperature array in K (surface first in lev).""" + sfc_t = 290.0 + lapse = np.linspace(0, 80, nlev) # decreases with altitude + return np.broadcast_to(sfc_t - lapse, (nlat, nlon, nlev)).copy().astype(np.float32) + + +def _make_q(nlat=_NLAT, nlon=_NLON, nlev=_NLEV) -> np.ndarray: + """Return a (nlat, nlon, nlev) specific humidity array (surface first in lev).""" + q_sfc = 0.01 + decay = np.exp(-np.linspace(0, 5, nlev)) + return np.broadcast_to(q_sfc * decay, (nlat, nlon, nlev)).copy().astype(np.float32) + + +def _make_z_ght(nlat=_NLAT, nlon=_NLON) -> np.ndarray: + """Return a (nlat, nlon, nlev) geopotential-height array, surface first.""" + return np.broadcast_to(_Z_GHT_SURFACE_FIRST, (nlat, nlon, _NLEV)).copy().astype(np.float64) + + +def write_pl_file_batch_format( + path: Path, + *, + lats: np.ndarray = _LATS_DESC, + lons: np.ndarray = _LONS, + levs_hpa: np.ndarray = _LEVS_HPA_DESC, + z_ght: np.ndarray | None = None, + time_dim_name: str = 'valid_time', + lev_dim_name: str = 'pressure_level', + z_as_geopotential: bool = False, +) -> Path: + """Write a synthetic pressure-level file in the format _batch_get_from_cds produces. + + * dims (time, lat, lon, lev) with z in geopotential-height metres (< 50 km). + * Optionally multiply z by g0 to produce a raw-geopotential file. + * ``time_dim_name`` and ``lev_dim_name`` let us probe the alternate coord names. + """ + nlat, nlon, nlev = len(lats), len(lons), len(levs_hpa) + if z_ght is None: + z_ght = np.broadcast_to(_Z_GHT_SURFACE_FIRST[:nlev], (nlat, nlon, nlev)).copy() + + t = _make_t(nlat, nlon, nlev) + q = _make_q(nlat, nlon, nlev) + z_out = z_ght * G0 if z_as_geopotential else z_ght.copy() + + t4d = t[np.newaxis] # (1, nlat, nlon, nlev) + q4d = q[np.newaxis] + z4d = z_out[np.newaxis] + + time_val = np.array(['2020-01-01'], dtype='datetime64[ns]') + + ds = xr.Dataset( + { + 'z': xr.Variable((time_dim_name, 'latitude', 'longitude', lev_dim_name), z4d), + 't': xr.Variable((time_dim_name, 'latitude', 'longitude', lev_dim_name), t4d), + 'q': xr.Variable((time_dim_name, 'latitude', 'longitude', lev_dim_name), q4d), + }, + coords={ + time_dim_name: time_val, + 'latitude': lats, + 'longitude': lons, + lev_dim_name: levs_hpa, + }, + ) + ds.to_netcdf(path) + return path + + +def write_ml_file(path: Path) -> Path: + """Write a minimal model-level file matching the format _makeDataCubes expects.""" + nlat, nlon, nlev = _NLAT, _NLON, _NLEV + lats = _LATS_ASC + lons = _LONS + + t = _make_t(nlat, nlon, nlev).transpose(2, 0, 1) # (nlev, nlat, nlon) + q = _make_q(nlat, nlon, nlev).transpose(2, 0, 1) + z = np.broadcast_to( + _Z_GHT_SURFACE_FIRST[:, np.newaxis, np.newaxis], + (nlev, nlat, nlon), + ).copy().astype(np.float32) + lnsp = np.full((nlat, nlon), np.log(95_000.0), dtype=np.float32) + + ds = xr.Dataset( + { + 't': xr.Variable(('valid_time', 'model_level', 'latitude', 'longitude'), t[np.newaxis]), + 'q': xr.Variable(('valid_time', 'model_level', 'latitude', 'longitude'), q[np.newaxis]), + 'z': xr.Variable(('valid_time', 'model_level', 'latitude', 'longitude'), z[np.newaxis]), + 'lnsp': xr.Variable(('valid_time', 'sfc_level', 'latitude', 'longitude'), lnsp[np.newaxis, np.newaxis]), + }, + coords={ + 'valid_time': np.array(['2020-01-01'], dtype='datetime64[ns]'), + 'model_level': np.arange(1, nlev + 1), + 'latitude': lats, + 'longitude': lons, + }, + ) + ds.to_netcdf(path) + return path + + +# --------------------------------------------------------------------------- +# Fixtures +# --------------------------------------------------------------------------- + +@pytest.fixture() +def era5(): + return ERA5() + + +@pytest.fixture() +def era5t(): + return ERA5T() + + +# --------------------------------------------------------------------------- +# Init / attribute tests +# --------------------------------------------------------------------------- + +class TestERA5Init: + def test_default_level_type_is_ml(self, era5: ERA5) -> None: + assert era5._model_level_type == 'ml' + + def test_name(self, era5: ERA5) -> None: + assert era5._Name == 'ERA-5' + + def test_humidity_type(self, era5: ERA5) -> None: + assert era5._humidityType == 'q' + + def test_projection(self, era5: ERA5) -> None: + assert era5._proj.to_epsg() == 4326 + + def test_valid_range_start(self, era5: ERA5) -> None: + assert era5._valid_range[0] == dt.datetime(1950, 1, 1, tzinfo=dt.timezone.utc) + + def test_load_weather_dispatches_to_pl(self, era5: ERA5, tmp_path: Path) -> None: + """load_weather should call _load_pressure_level for pl type.""" + f = write_pl_file_batch_format(tmp_path / 'test.nc') + era5.setLevelType('pl') + era5.set_latlon_bounds([48.0, 53.0, 9.0, 15.0]) + era5.files = [f] + + era5.load_weather(f) + + assert era5._t is not None + assert era5._q is not None + + +class TestERA5TInit: + def test_name(self, era5t: ERA5T) -> None: + assert era5t._Name == 'ERA-5T' + + def test_expver(self, era5t: ERA5T) -> None: + assert era5t._expver == '0005' + + def test_dataset(self, era5t: ERA5T) -> None: + assert era5t._dataset == 'era5t' + + def test_inherits_ml_level_type(self, era5t: ERA5T) -> None: + assert era5t._model_level_type == 'ml' + + def test_lag_time_is_one_day(self, era5t: ERA5T) -> None: + assert era5t._lag_time == dt.timedelta(days=1) + + def test_valid_range_extends_to_now(self, era5t: ERA5T) -> None: + # ERA5T covers up to "now"; check the upper bound is in the future relative to + # a known past date + past = dt.datetime(2020, 1, 1, tzinfo=dt.timezone.utc) + assert era5t._valid_range[1] > past + + +# --------------------------------------------------------------------------- +# _load_pressure_level tests +# --------------------------------------------------------------------------- + +class TestLoadPressureLevel: + """All variants of _load_pressure_level — axis ordering, unit guards, coord names.""" + + def _setup_and_load(self, path: Path, era5: ERA5) -> None: + era5.set_latlon_bounds([48.0, 53.0, 9.0, 15.0]) + era5._load_pressure_level(path) + + # --- shape correctness ------------------------------------------------- + + def test_output_shapes_are_lat_lon_lev(self, era5: ERA5, tmp_path: Path) -> None: + f = write_pl_file_batch_format(tmp_path / 'test.nc') + self._setup_and_load(f, era5) + + assert era5._t.shape == (_NLAT, _NLON, _NLEV) + assert era5._q.shape == (_NLAT, _NLON, _NLEV) + assert era5._zs.shape == (_NLAT, _NLON, _NLEV) + assert era5._p.shape == (_NLAT, _NLON, _NLEV) + + # --- lat/lon ordering -------------------------------------------------- + + def test_lats_are_ascending_south_to_north(self, era5: ERA5, tmp_path: Path) -> None: + """Loader must flip descending CDS lats to ascending.""" + f = write_pl_file_batch_format(tmp_path / 'test.nc', lats=_LATS_DESC) + self._setup_and_load(f, era5) + + # _lats is a 2D meshgrid; first row should be the southernmost + assert era5._lats[0, 0] < era5._lats[-1, 0] + + def test_lats_already_ascending_unchanged(self, era5: ERA5, tmp_path: Path) -> None: + f = write_pl_file_batch_format(tmp_path / 'test.nc', lats=_LATS_ASC) + self._setup_and_load(f, era5) + + assert era5._lats[0, 0] < era5._lats[-1, 0] + + def test_lons_gt_180_normalized(self, era5: ERA5, tmp_path: Path) -> None: + """Longitudes > 180 must be shifted to [-180, 180].""" + lons_wrapped = np.array([350.0, 351.0, 352.0, 353.0, 354.0]) + f = write_pl_file_batch_format(tmp_path / 'test.nc', lons=lons_wrapped) + # Widen bounds to include these lons + era5.set_latlon_bounds([48.0, 53.0, -15.0, 0.0]) + era5._load_pressure_level(f) + + assert np.all(era5._lons <= 180.0) + assert np.all(era5._lons >= -180.0) + + def test_descending_lons_flipped(self, era5: ERA5, tmp_path: Path) -> None: + lons_desc = np.array([14.0, 13.0, 12.0, 11.0, 10.0]) + f = write_pl_file_batch_format(tmp_path / 'test.nc', lons=lons_desc) + self._setup_and_load(f, era5) + + assert era5._lons[0, 0] < era5._lons[0, -1] + + # --- level ordering (pressure / height) -------------------------------- + + def test_pressure_surface_first(self, era5: ERA5, tmp_path: Path) -> None: + """p[...,0] must be the largest pressure (surface).""" + f = write_pl_file_batch_format(tmp_path / 'test.nc') + self._setup_and_load(f, era5) + + assert era5._p[0, 0, 0] > era5._p[0, 0, -1] + + def test_pressure_toa_last(self, era5: ERA5, tmp_path: Path) -> None: + f = write_pl_file_batch_format(tmp_path / 'test.nc') + self._setup_and_load(f, era5) + + # Top-of-atmosphere pressure should be <= 1000 Pa + assert era5._p[0, 0, -1] <= 1_000.0 + + def test_ascending_hpa_levels_are_flipped(self, era5: ERA5, tmp_path: Path) -> None: + """When CDS returns levels ascending (TOA-first), the loader must flip them.""" + levs_asc = _LEVS_HPA_DESC[::-1] # 10 → 1000 hPa (TOA first) + z_toa_first = _Z_GHT_SURFACE_FIRST[::-1] # heights matching TOA-first levels + z4d = np.broadcast_to(z_toa_first, (_NLAT, _NLON, _NLEV)).copy() + f = write_pl_file_batch_format(tmp_path / 'test.nc', levs_hpa=levs_asc, z_ght=z4d) + self._setup_and_load(f, era5) + + assert era5._p[0, 0, 0] > era5._p[0, 0, -1], 'pressure must be surface-first' + + def test_heights_monotonically_increasing(self, era5: ERA5, tmp_path: Path) -> None: + """Heights must increase from surface to TOA (index 0 → -1).""" + f = write_pl_file_batch_format(tmp_path / 'test.nc') + self._setup_and_load(f, era5) + + dz = np.diff(era5._zs, axis=2) + assert np.all(dz > 0) + + # --- z unit guard (geopotential vs geopotential-height) ---------------- + + def test_geopotential_z_divided_by_g0(self, era5: ERA5, tmp_path: Path) -> None: + """When z is stored as geopotential (m²/s²) the loader divides by g0.""" + f_ght = write_pl_file_batch_format(tmp_path / 'ght.nc', z_as_geopotential=False) + f_gp = write_pl_file_batch_format(tmp_path / 'gp.nc', z_as_geopotential=True) + + era5_ght = ERA5() + era5_ght.set_latlon_bounds([48.0, 53.0, 9.0, 15.0]) + era5_ght._load_pressure_level(f_ght) + + era5_gp = ERA5() + era5_gp.set_latlon_bounds([48.0, 53.0, 9.0, 15.0]) + era5_gp._load_pressure_level(f_gp) + + # Heights should be (approximately) the same once g0 division is applied + np.testing.assert_allclose(era5_ght._zs, era5_gp._zs, rtol=1e-4) + + def test_geopotential_height_file_not_double_divided(self, era5: ERA5, tmp_path: Path) -> None: + """A file with z already in metres must NOT be divided again by g0.""" + f = write_pl_file_batch_format(tmp_path / 'test.nc') + self._setup_and_load(f, era5) + + # The highest point on our test grid should be ~31 km, not ~3 km + assert era5._zs.max() > 25_000 + + # --- alternate coordinate names ---------------------------------------- + + def test_time_dim_named_time(self, era5: ERA5, tmp_path: Path) -> None: + """Loader must handle 'time' as the time-dimension name.""" + f = write_pl_file_batch_format(tmp_path / 'test.nc', time_dim_name='time') + self._setup_and_load(f, era5) + + assert era5._t.shape == (_NLAT, _NLON, _NLEV) + + def test_lev_dim_named_level(self, era5: ERA5, tmp_path: Path) -> None: + """Loader must handle 'level' as the level-dimension name.""" + f = write_pl_file_batch_format(tmp_path / 'test.nc', lev_dim_name='level') + self._setup_and_load(f, era5) + + assert era5._t.shape == (_NLAT, _NLON, _NLEV) + + # --- crude physics sanity check ---------------------------------------- + + def test_hydrostatic_ztd_near_sea_level_is_plausible(self, era5: ERA5, tmp_path: Path) -> None: + """Integrate k1 * p/t to get a rough hydrostatic ZTD. + + This test would have caught the original bug: the double g0 division + compressed the atmosphere from ~48 km to ~4.9 km, giving a ZTD + of ~1 m instead of the expected ~2.2 m. + """ + f = write_pl_file_batch_format(tmp_path / 'test.nc') + self._setup_and_load(f, era5) + + # Pick the column closest to sea level + i, j = np.unravel_index(np.argmin(era5._zs[..., 0]), era5._zs[..., 0].shape) + zz = era5._zs[i, j] + pp = era5._p[i, j] + tt = era5._t[i, j] + + n_hydro = 0.776 * pp / tt * 1e-6 + ztd_h = np.trapezoid(n_hydro, zz) + + assert 1.5 < ztd_h < 3.0, ( + f'Hydrostatic ZTD {ztd_h:.3f} m is outside the expected range [1.5, 3.0] m. ' + 'This may indicate an axis-ordering or unit bug.' + ) + + +# --------------------------------------------------------------------------- +# _makeDataCubes / _load_model_level tests +# --------------------------------------------------------------------------- + +class TestMakeDataCubes: + def test_returns_correct_shapes(self, era5: ERA5, tmp_path: Path) -> None: + f = write_ml_file(tmp_path / 'ml.nc') + era5.set_latlon_bounds([48.0, 53.0, 9.0, 15.0]) + era5.setLevelType('ml') + + lats, lons, _, _, t, q, lnsp, z = era5._makeDataCubes(f) + + assert t.shape == (_NLEV, _NLAT, _NLON) + assert q.shape == (_NLEV, _NLAT, _NLON) + assert lnsp.shape == (_NLAT, _NLON) + assert z.shape == (_NLEV, _NLAT, _NLON) + assert len(lats) == _NLAT + assert len(lons) == _NLON + + def test_raises_when_mask_excludes_all_data(self, era5: ERA5, tmp_path: Path) -> None: + """Bounds that exclude all data must raise an error (RuntimeError or ValueError).""" + f = write_ml_file(tmp_path / 'ml.nc') + era5.set_latlon_bounds([80.0, 85.0, 90.0, 95.0]) + era5.setLevelType('ml') + + # xarray + netCDF4 raises ValueError when trying to index a 0-sized + # dimension; _makeDataCubes raises RuntimeError for the explicit check. + # Either is acceptable — both indicate no usable data. + with pytest.raises((RuntimeError, ValueError)): + era5._makeDataCubes(f) + + +def _mock_calculategeoh(nlat, nlon, nlev): + """Return a mock for _calculategeoh that produces plausible (nlev, nlat, nlon) arrays. + + _calculategeoh requires num_levels == t.shape[0] (137 for ML mode), so a + synthetic 6-level file would crash the real implementation. We mock it to + keep the ML loader tests fast and self-contained. + + _calculategeoh returns data with model levels numbered from the HIGHEST + elevation (TOA) to the LOWEST (surface), so index 0 is TOA. The loader + then transposes and flips to produce surface-first (lat, lon, lev) arrays. + """ + # TOA-first ordering: highest height / lowest pressure at index 0 + z_toa_first = _Z_GHT_SURFACE_FIRST[:nlev][::-1] # [31100, ..., 100] + p_toa_first = _LEVS_PA_DESC[:nlev][::-1] # [1000, ..., 100000] + z_3d = np.broadcast_to( + z_toa_first[:, np.newaxis, np.newaxis], + (nlev, nlat, nlon), + ).copy() + p_3d = np.broadcast_to( + p_toa_first[:, np.newaxis, np.newaxis], + (nlev, nlat, nlon), + ).copy() + return z_3d, p_3d, z_3d + + +class TestLoadModelLevel: + """Tests for _load_model_level. + + _calculategeoh is mocked here because it requires num_levels == t.shape[0] + (137 for ERA5 ML mode) which would force an impractically large synthetic file. + The physics of _calculategeoh are tested separately in test_util.py / calcgeoh. + """ + + def _load(self, era5: ERA5, path: Path) -> None: + era5.set_latlon_bounds([48.0, 53.0, 9.0, 15.0]) + era5.setLevelType('ml') + mock_return = _mock_calculategeoh(_NLAT, _NLON, _NLEV) + with patch.object(era5, '_calculategeoh', return_value=mock_return): + era5._load_model_level(path) + + def test_output_shapes_are_lat_lon_lev(self, era5: ERA5, tmp_path: Path) -> None: + f = write_ml_file(tmp_path / 'ml.nc') + self._load(era5, f) + + assert era5._t.shape == (_NLAT, _NLON, _NLEV) + assert era5._q.shape == (_NLAT, _NLON, _NLEV) + assert era5._zs.shape == (_NLAT, _NLON, _NLEV) + assert era5._p.shape == (_NLAT, _NLON, _NLEV) + + def test_heights_monotonically_increasing(self, era5: ERA5, tmp_path: Path) -> None: + f = write_ml_file(tmp_path / 'ml.nc') + self._load(era5, f) + + dz = np.diff(era5._zs, axis=2) + assert np.all(dz > 0) + + def test_pressure_surface_first(self, era5: ERA5, tmp_path: Path) -> None: + f = write_ml_file(tmp_path / 'ml.nc') + self._load(era5, f) + + assert era5._p[0, 0, 0] > era5._p[0, 0, -1] + + +# --------------------------------------------------------------------------- +# batch_fetch tests (ERA5) +# --------------------------------------------------------------------------- + +class TestBatchFetch: + def test_skips_already_downloaded_file(self, era5: ERA5, tmp_path: Path) -> None: + """Files that already exist must not be re-queued.""" + existing = tmp_path / 'already_there.nc' + existing.touch() + + era5.set_latlon_bounds([48.0, 53.0, 9.0, 15.0]) + + with patch.object(era5, '_batch_get_from_cds') as mock_batch: + era5.batch_fetch([(dt.datetime(2020, 1, 1, 0), existing)]) + + mock_batch.assert_not_called() + + def test_queues_missing_file(self, era5: ERA5, tmp_path: Path) -> None: + """A file that does not exist must be passed to _batch_get_from_cds.""" + missing = tmp_path / 'missing.nc' + + era5.set_latlon_bounds([48.0, 53.0, 9.0, 15.0]) + + with patch.object(era5, '_batch_get_from_cds') as mock_batch: + era5.batch_fetch([(dt.datetime(2020, 1, 1, 0), missing)]) + + mock_batch.assert_called_once() + queued_paths = [p for _, p in mock_batch.call_args[0][0]] + assert missing in queued_paths + + def test_mixed_existing_and_missing(self, era5: ERA5, tmp_path: Path) -> None: + """Only missing files must be queued; existing ones silently skipped.""" + existing = tmp_path / 'exists.nc' + existing.touch() + missing = tmp_path / 'new.nc' + + era5.set_latlon_bounds([48.0, 53.0, 9.0, 15.0]) + + with patch.object(era5, '_batch_get_from_cds') as mock_batch: + era5.batch_fetch([ + (dt.datetime(2020, 1, 1, 0), existing), + (dt.datetime(2020, 1, 2, 0), missing), + ]) + + mock_batch.assert_called_once() + queued = mock_batch.call_args[0][0] + queued_paths = [p for _, p in queued] + assert missing in queued_paths + assert existing not in queued_paths + + def test_empty_input_does_not_call_batch(self, era5: ERA5, tmp_path: Path) -> None: + era5.set_latlon_bounds([48.0, 53.0, 9.0, 15.0]) + with patch.object(era5, '_batch_get_from_cds') as mock_batch: + era5.batch_fetch([]) + mock_batch.assert_not_called() + + def test_time_rounding(self, era5: ERA5, tmp_path: Path) -> None: + """Times that are not on an exact hour boundary must be rounded.""" + missing = tmp_path / 'rounded.nc' + + era5.set_latlon_bounds([48.0, 53.0, 9.0, 15.0]) + + with patch.object(era5, '_batch_get_from_cds') as mock_batch: + era5.batch_fetch([(dt.datetime(2020, 1, 1, 0, 31), missing)]) + + mock_batch.assert_called_once() + queued_dt, _ = mock_batch.call_args[0][0][0] + assert queued_dt.minute == 0 + + +# --------------------------------------------------------------------------- +# ERA5T inherits ERA5 loading unchanged +# --------------------------------------------------------------------------- + +class TestERA5TLoading: + def test_load_weather_uses_pl_loader(self, era5t: ERA5T, tmp_path: Path) -> None: + f = write_pl_file_batch_format(tmp_path / 'test.nc') + era5t.setLevelType('pl') + era5t.set_latlon_bounds([48.0, 53.0, 9.0, 15.0]) + + era5t.load_weather(f) + + assert era5t._t is not None + assert era5t._t.shape == (_NLAT, _NLON, _NLEV) + + def test_batch_fetch_inherited(self, era5t: ERA5T, tmp_path: Path) -> None: + """ERA5T inherits batch_fetch from ERA5.""" + assert callable(getattr(era5t, 'batch_fetch', None)) diff --git a/test/test_gnss.py b/test/test_gnss.py index 18bad7ed2..f0173f72b 100644 --- a/test/test_gnss.py +++ b/test/test_gnss.py @@ -14,6 +14,14 @@ SCENARIO2_DIR = os.path.join(TEST_DIR, "scenario_2") +# Tests that perform real bulk GNSS-delay downloads can take many minutes and +# hang CI. They are skipped unless RAIDER_RUN_NETWORK_TESTS is set. +requires_network = pytest.mark.skipif( + not os.environ.get("RAIDER_RUN_NETWORK_TESTS"), + reason="Real network download; set RAIDER_RUN_NETWORK_TESTS=1 to run", +) + + def file_len(path: Path) -> int: with path.open('rb') as f: return sum(1 for _ in f) @@ -91,6 +99,7 @@ def test_concatDelayFiles(tmp_path, temp_file): assert file_len(out_path) == file_length +@requires_network def test_get_stats_by_llh2(): stations = get_stats_by_llh(llhBox=[10, 18, -93, -88]) assert isinstance(stations, pd.DataFrame) @@ -114,11 +123,13 @@ def test_download_tropo_delays1(): 2022], gps_repo='dummy_repo') -def test_download_tropo_delays2(): +@requires_network +def test_download_tropo_delays_no_data(): with pytest.raises(NoStationDataFoundError): download_tropo_delays(stats=['dummy_station'], years=[2022]) +@requires_network def test_download_tropo_delays2(tmp_path): with pushd(tmp_path): stations, output_file = get_station_list( diff --git a/test/test_llreader.py b/test/test_llreader.py index b22ad3a73..24fc37b6c 100644 --- a/test/test_llreader.py +++ b/test/test_llreader.py @@ -1,9 +1,12 @@ +import logging import os from pathlib import Path import pytest +from unittest.mock import MagicMock import numpy as np import pandas as pd +import pyproj from test import GEOM_DIR, TEST_DIR from pyproj import CRS @@ -12,7 +15,8 @@ from RAiDER.utilFcns import rio_open from RAiDER.llreader import ( - StationFile, RasterRDR, BoundingBox, GeocodedFile, bounds_from_latlon_rasters, bounds_from_csv + StationFile, RasterRDR, BoundingBox, GeocodedFile, bounds_from_latlon_rasters, + bounds_from_csv, _parse_crs, _ellipsoidal_to_geometric, ) SCENARIO0_DIR = TEST_DIR / "scenario_0" @@ -147,3 +151,285 @@ def test_GeocodedFile(): x,y = aoi.readLL() assert z.shape == (569,558) assert x.shape == z.shape + + +# --------------------------------------------------------------------------- +# Fixtures shared by CRS / geoid tests +# --------------------------------------------------------------------------- + +@pytest.fixture +def tmp_station_file(tmp_path): + """Three-station CSV with known ellipsoidal heights across CONUS latitudes.""" + csv = tmp_path / 'stations.csv' + csv.write_text( + 'ID,Lat,Lon,Hgt_m\n' + 'STA1,34.0,-118.0,50.0\n' + 'STA2,36.0,-115.0,150.0\n' + 'STA3,49.0,-122.0,195.28\n' + ) + return csv + + +@pytest.fixture +def conus_lats(): + return np.array([34.0, 36.0, 49.0]) + + +@pytest.fixture +def conus_lons(): + return np.array([-118.0, -115.0, -122.0]) + + +@pytest.fixture +def conus_heights(): + return np.array([50.0, 150.0, 195.28]) + + +# --------------------------------------------------------------------------- +# _parse_crs +# --------------------------------------------------------------------------- + +def test_parse_crs_int(): + assert _parse_crs(4326).to_epsg() == 4326 + + +def test_parse_crs_string_bare(): + assert _parse_crs('4979').to_epsg() == 4979 + + +def test_parse_crs_string_with_prefix(): + assert _parse_crs('EPSG:4979').to_epsg() == 4979 + + +def test_parse_crs_crs_object_is_identity(): + """Passing a CRS object should return the same object unchanged.""" + src = CRS.from_epsg(32631) + assert _parse_crs(src) is src + + +def test_parse_crs_invalid_raises(): + with pytest.raises(Exception): + _parse_crs('definitely_not_a_crs') + + +# --------------------------------------------------------------------------- +# _ellipsoidal_to_geometric – unit tests (no network required) +# --------------------------------------------------------------------------- + +def test_ellipsoidal_to_geometric_2d_crs_is_noop(conus_lats, conus_lons, conus_heights): + """A 2D CRS (EPSG:4326) must return heights unchanged without touching PROJ.""" + result = _ellipsoidal_to_geometric( + conus_lats, conus_lons, conus_heights, CRS.from_epsg(4326) + ) + assert np.array_equal(result, conus_heights) + + +def test_ellipsoidal_to_geometric_cdn_fallback_applied( + conus_lats, conus_lons, conus_heights, monkeypatch +): + """Noop local grid → CDN enabled → real corrections applied.""" + correction = np.full_like(conus_heights, 20.0) + h_expected = conus_heights + correction + call_count = {'n': 0} + + def mock_from_crs(*args, **kwargs): + call_count['n'] += 1 + t = MagicMock() + if call_count['n'] == 1: # first build: local grid missing → noop + t.to_proj4.return_value = '+proj=noop' + t.transform.return_value = (conus_lons, conus_lats, conus_heights.copy()) + else: # second build: after CDN enabled → real grid + t.to_proj4.return_value = None + t.transform.return_value = (conus_lons, conus_lats, h_expected.copy()) + return t + + enabled_states = [] + monkeypatch.setattr(pyproj.Transformer, 'from_crs', mock_from_crs) + monkeypatch.setattr(pyproj.network, 'is_network_enabled', lambda: False) + monkeypatch.setattr(pyproj.network, 'set_network_enabled', + lambda v: enabled_states.append(v)) + + result = _ellipsoidal_to_geometric( + conus_lats, conus_lons, conus_heights, CRS.from_epsg(4979) + ) + + assert np.allclose(result, h_expected) + assert call_count['n'] == 2, 'Transformer should be built twice (local then CDN)' + assert True in enabled_states, 'CDN should have been enabled' + assert False in enabled_states, 'CDN state should have been restored' + + +def test_ellipsoidal_to_geometric_cdn_also_noop_warns( + conus_lats, conus_lons, conus_heights, monkeypatch, caplog +): + """When CDN is also unavailable (still noop), return heights with a warning.""" + noop_mock = MagicMock() + noop_mock.to_proj4.return_value = '+proj=noop' + noop_mock.transform.return_value = (conus_lons, conus_lats, conus_heights.copy()) + + monkeypatch.setattr(pyproj.Transformer, 'from_crs', lambda *a, **kw: noop_mock) + monkeypatch.setattr(pyproj.network, 'is_network_enabled', lambda: False) + monkeypatch.setattr(pyproj.network, 'set_network_enabled', lambda v: None) + + with caplog.at_level(logging.WARNING, logger='RAiDER'): + result = _ellipsoidal_to_geometric( + conus_lats, conus_lons, conus_heights, CRS.from_epsg(4979) + ) + + assert np.array_equal(result, conus_heights) + assert 'EGM96 grid unavailable' in caplog.text + + +def test_ellipsoidal_to_geometric_exception_returns_unchanged( + conus_lats, conus_lons, conus_heights, monkeypatch, caplog +): + """An unexpected PROJ exception must return heights unchanged with a warning.""" + def raising(*args, **kwargs): + raise RuntimeError('simulated PROJ error') + + monkeypatch.setattr(pyproj.Transformer, 'from_crs', raising) + + with caplog.at_level(logging.WARNING, logger='RAiDER'): + result = _ellipsoidal_to_geometric( + conus_lats, conus_lons, conus_heights, CRS.from_epsg(4979) + ) + + assert np.array_equal(result, conus_heights) + assert 'conversion failed' in caplog.text.lower() + + +# --------------------------------------------------------------------------- +# _ellipsoidal_to_geometric – integration test (requires PROJ CDN or local grids) +# --------------------------------------------------------------------------- + +@pytest.mark.network +def test_ellipsoidal_to_geometric_egm96_conus_values( + conus_lats, conus_lons, conus_heights +): + """EGM96 orthometric heights in CONUS are 10–45 m larger than WGS84 ellipsoidal.""" + result = _ellipsoidal_to_geometric( + conus_lats, conus_lons, conus_heights, CRS.from_epsg(4979) + ) + corrections = result - conus_heights + # EGM96 geoid is 15–35 m below the WGS84 ellipsoid in CONUS, so + # orthometric heights are larger than ellipsoidal heights. + assert np.all(corrections > 10), f'Expected >10 m corrections; got {corrections}' + assert np.all(corrections < 45), f'Expected <45 m corrections; got {corrections}' + + +# --------------------------------------------------------------------------- +# StationFile – crs parameter and update_crs method +# --------------------------------------------------------------------------- + +def test_station_file_default_crs_is_4326(tmp_station_file): + sf = StationFile(tmp_station_file) + assert sf._crs.to_epsg() == 4326 + + +def test_station_file_crs_int(tmp_station_file): + sf = StationFile(tmp_station_file, crs=4979) + assert sf._crs.to_epsg() == 4979 + + +def test_station_file_crs_string(tmp_station_file): + sf = StationFile(tmp_station_file, crs='4979') + assert sf._crs.to_epsg() == 4979 + + +def test_station_file_crs_epsg_string(tmp_station_file): + sf = StationFile(tmp_station_file, crs='EPSG:4979') + assert sf._crs.to_epsg() == 4979 + + +def test_station_file_crs_object(tmp_station_file): + crs = CRS.from_epsg(4979) + sf = StationFile(tmp_station_file, crs=crs) + assert sf._crs.to_epsg() == 4979 + + +def test_update_crs_changes_datum(tmp_station_file): + sf = StationFile(tmp_station_file) + assert sf._crs.to_epsg() == 4326 + sf.update_crs(4979) + assert sf._crs.to_epsg() == 4979 + + +def test_update_crs_string_input(tmp_station_file): + sf = StationFile(tmp_station_file) + sf.update_crs('EPSG:4979') + assert sf._crs.to_epsg() == 4979 + + +def test_update_crs_crs_object(tmp_station_file): + sf = StationFile(tmp_station_file) + sf.update_crs(CRS.from_epsg(4979)) + assert sf._crs.to_epsg() == 4979 + + +# --------------------------------------------------------------------------- +# StationFile.readZ – height conversion plumbing +# --------------------------------------------------------------------------- + +def test_readZ_default_crs_heights_unchanged(tmp_station_file): + """crs=4326 must return Hgt_m values from the CSV unmodified.""" + z = StationFile(tmp_station_file).readZ() + assert np.allclose(z, [50.0, 150.0, 195.28]) + + +def test_readZ_3d_crs_calls_conversion(tmp_station_file, monkeypatch): + """crs=4979 must route Hgt_m through _ellipsoidal_to_geometric.""" + import RAiDER.llreader as lr + + captured = {} + + def mock_convert(lats, lons, heights, crs): + captured['heights'] = heights.copy() + captured['crs_epsg'] = crs.to_epsg() + return heights + 25.0 + + monkeypatch.setattr(lr, '_ellipsoidal_to_geometric', mock_convert) + + z = StationFile(tmp_station_file, crs=4979).readZ() + + assert captured['crs_epsg'] == 4979 + assert np.allclose(captured['heights'], [50.0, 150.0, 195.28]) + assert np.allclose(z, [75.0, 175.0, 220.28]) + + +def test_readZ_update_crs_uses_new_datum(tmp_station_file, monkeypatch): + """update_crs before readZ must pass the updated CRS to the conversion.""" + import RAiDER.llreader as lr + + captured_epsg = [] + + def mock_convert(lats, lons, heights, crs): + captured_epsg.append(crs.to_epsg()) + return heights + + monkeypatch.setattr(lr, '_ellipsoidal_to_geometric', mock_convert) + + sf = StationFile(tmp_station_file) # default 4326 + sf.update_crs(4979) + sf.readZ() + + assert captured_epsg == [4979] + + +def test_readZ_2d_crs_conversion_is_noop(tmp_station_file, monkeypatch): + """_ellipsoidal_to_geometric is still called for 2D CRS but returns heights unchanged.""" + import RAiDER.llreader as lr + + call_count = {'n': 0} + + def mock_convert(lats, lons, heights, crs): + call_count['n'] += 1 + # The function itself handles the 2D noop, but verify it receives crs=4326 + assert crs.to_epsg() == 4326 + return heights # same as real behaviour for 2D + + monkeypatch.setattr(lr, '_ellipsoidal_to_geometric', mock_convert) + + z = StationFile(tmp_station_file).readZ() + + assert call_count['n'] == 1 + assert np.allclose(z, [50.0, 150.0, 195.28]) diff --git a/test/test_synthetic.py b/test/test_synthetic.py index 1bd4b3059..7b14bb98f 100644 --- a/test/test_synthetic.py +++ b/test/test_synthetic.py @@ -19,6 +19,17 @@ from test import ORB_DIR, TEST_DIR, WM_DIR, pushd +# These tests fall through to a real ERA5 download from the CDS API, which +# fails with an HTTP 401 in CI (a placeholder ~/.cdsapirc is often present, so +# detecting credential *presence* is not enough -- they have to be valid). +# Gate them behind an explicit opt-in flag so CI always skips and only someone +# with working CDS credentials runs them deliberately. +requires_cds = pytest.mark.skipif( + not os.environ.get("RAIDER_RUN_NETWORK_TESTS"), + reason="Real CDS download; set RAIDER_RUN_NETWORK_TESTS=1 (with valid credentials) to run", +) + + def update_model(wm_file: str, wm_eq_type: str, wm_dir: str = "weather_files_synth"): """Update weather model file by the equation to test, write it to disk. @@ -41,6 +52,7 @@ def update_model(wm_file: str, wm_eq_type: str, wm_dir: str = "weather_files_syn e = ds["e"] if wm_eq_type == "hydro": p = t + e = e * 0 # keep N_hydro = k1 constant under the k1*P/Tv formulation elif wm_eq_type == "wet_linear": e = t Obj._k3 = 0 @@ -214,6 +226,7 @@ def test_dl_real(tmp_path, region, mod="ERA5"): assert proc.returncode == 0, 'RAiDER did not complete successfully' +@requires_cds @pytest.mark.parametrize("region", "AK LA Fort".split()) def test_hydrostatic_eq(tmp_path, region, mod="ERA-5"): """Test hydrostatic equation: Hydro Refractivity = k1 * (Pressure/Temp). @@ -277,6 +290,7 @@ def test_hydrostatic_eq(tmp_path, region, mod="ERA-5"): del da +@requires_cds @pytest.mark.parametrize("region", "AK LA Fort".split()) def test_wet_eq_linear(tmp_path, region, mod="ERA-5"): """Test linear part of wet equation. @@ -353,6 +367,7 @@ def test_wet_eq_linear(tmp_path, region, mod="ERA-5"): shutil.rmtree(dir_to_del) +@requires_cds @pytest.mark.parametrize("region", "AK LA Fort".split()) def test_wet_eq_nonlinear(tmp_path, region, mod="ERA-5"): """Test the nonlinear part of the wet equation.""" diff --git a/test/test_weather_model.py b/test/test_weather_model.py index 205673f1a..25c4ea0e7 100644 --- a/test/test_weather_model.py +++ b/test/test_weather_model.py @@ -66,22 +66,25 @@ def load_weather(self, *args, **kwargs) -> None: # noqa: D102 self._xs = np.arange(-3, 4) + _LON0 self._zs = np.linspace(0, 1e5, _N_Z) self._t = np.ones((len(self._ys), len(self._xs), _N_Z)) - self._e = self._t.copy() - self._e[:, 3:, :] = 2 - _p = np.arange(31, -1, -1) + # Exponential p and e profiles with a common scale height (like the + # real atmosphere), so both refractivities are single exponentials + # and the true ZTDs have exact closed forms independent of the + # quadrature: int_z^ztop N0*exp(-z'/H) dz' = H * (N(z) - N(ztop)) + _H = 2.0e4 + _decay = np.exp(-self._zs / _H) + self._e = np.ones(self._t.shape) * _decay + self._e[:, 3:, :] *= 2 + _p = 32 * _decay self._p = np.broadcast_to(_p, self._t.shape) - self._true_hydro_refr = np.broadcast_to(_p, (self._t.shape)) - self._true_wet_ztd = 1e-6 * 2 * np.broadcast_to(np.flip(self._zs), (self._t.shape)) - self._true_wet_ztd[:, 3:] = 2 * self._true_wet_ztd[:, 3:] + # wet refractivity: k2*e/t + k3*e/t^2 = 2*e (k2 = k3 = t = 1) + self._true_wet_refr = 2 * self._e + self._true_wet_ztd = 1e-6 * 2 * _H * (self._e - self._e[..., -1:]) - self._true_hydro_ztd = np.zeros(self._t.shape) - for layer in range(len(self._zs)): - self._true_hydro_ztd[:, :, layer] = 1e-6 * 0.5 * (self._zs[-1] - self._zs[layer]) * _p[layer] - - self._true_wet_refr = 2 * np.ones(self._t.shape) - self._true_wet_refr[:, 3:] = 4 + # hydrostatic refractivity: k1*(p - (1 - Rd/Rv)*e)/t (k1 = t = 1) + self._true_hydro_refr = self._p - (1 - self._R_d / self._R_v) * self._e + self._true_hydro_ztd = 1e-6 * _H * (self._true_hydro_refr - self._true_hydro_refr[..., -1:]) def interpWet(self): # noqa: ANN201, D102 _ifWet = rgi((self._ys, self._xs, self._zs), self._true_wet_refr) diff --git a/tools/RAiDER/aria/prepFromGUNW.py b/tools/RAiDER/aria/prepFromGUNW.py index bca72baf8..9bf2d187a 100644 --- a/tools/RAiDER/aria/prepFromGUNW.py +++ b/tools/RAiDER/aria/prepFromGUNW.py @@ -360,6 +360,10 @@ def main(args: CalcDelaysArgs) -> tuple[Path, float]: raider_cfg = { 'weather_model': args.weather_model, + # Optionally override the model's default vertical level representation + # (model/native vs. pressure levels). Left out when not specified so the + # model's built-in default is used. + **({'weather_model_levels': args.model_levels} if getattr(args, 'model_levels', None) is not None else {}), 'look_dir': GUNWObj.look_dir, 'aoi_group': {'bounding_box': GUNWObj.SNWE}, 'height_group': {'height_levels': GUNWObj.heights}, diff --git a/tools/RAiDER/aria/types.py b/tools/RAiDER/aria/types.py index 56034ddd0..9a80c623e 100644 --- a/tools/RAiDER/aria/types.py +++ b/tools/RAiDER/aria/types.py @@ -15,6 +15,7 @@ class CalcDelaysArgsUnparsed(argparse.Namespace): api_key: Optional[str] interpolate_time: TimeInterpolationMethod output_directory: Path + model_levels: Optional[str] class CalcDelaysArgs(argparse.Namespace): bucket: Optional[str] @@ -26,3 +27,4 @@ class CalcDelaysArgs(argparse.Namespace): api_key: Optional[str] interpolate_time: TimeInterpolationMethod output_directory: Path + model_levels: Optional[str] diff --git a/tools/RAiDER/cli/examples/example_LA_GNSS/example_LA_GNSS.yaml b/tools/RAiDER/cli/examples/example_LA_GNSS/example_LA_GNSS.yaml index 148e1fd62..50e6d24f7 100644 --- a/tools/RAiDER/cli/examples/example_LA_GNSS/example_LA_GNSS.yaml +++ b/tools/RAiDER/cli/examples/example_LA_GNSS/example_LA_GNSS.yaml @@ -31,6 +31,12 @@ look_dir: right ## for more details and information on licensing weather_model: GMAO +## OPTIONAL: Vertical level representation used by the weather model. +## FORMATS: string, one of 'ml'/'model' (model/native levels) or 'pl'/'pressure' (pressure levels) +## If left blank, each model's built-in default is used. Note that not every +## model supports every level type (GMAO supports only model levels). +weather_model_levels: + ########## 2. Date ## REQUIRED: TRUE diff --git a/tools/RAiDER/cli/examples/example_LA_bbox/example_LA_bbox.yaml b/tools/RAiDER/cli/examples/example_LA_bbox/example_LA_bbox.yaml index e1ce8a8e9..a58a19a5f 100644 --- a/tools/RAiDER/cli/examples/example_LA_bbox/example_LA_bbox.yaml +++ b/tools/RAiDER/cli/examples/example_LA_bbox/example_LA_bbox.yaml @@ -31,6 +31,12 @@ look_dir: right ## for more details and information on licensing weather_model: GMAO +## OPTIONAL: Vertical level representation used by the weather model. +## FORMATS: string, one of 'ml'/'model' (model/native levels) or 'pl'/'pressure' (pressure levels) +## If left blank, each model's built-in default is used. Note that not every +## model supports every level type (GMAO supports only model levels). +weather_model_levels: + ########## 2. Date ## REQUIRED: TRUE diff --git a/tools/RAiDER/cli/examples/example_UK_isce/example_UK_isce.yaml b/tools/RAiDER/cli/examples/example_UK_isce/example_UK_isce.yaml index caf3127d5..c655a8e98 100644 --- a/tools/RAiDER/cli/examples/example_UK_isce/example_UK_isce.yaml +++ b/tools/RAiDER/cli/examples/example_UK_isce/example_UK_isce.yaml @@ -31,6 +31,12 @@ look_dir: right ## for more details and information on licensing weather_model: ERA5 +## OPTIONAL: Vertical level representation used by the weather model. +## FORMATS: string, one of 'ml'/'model' (model/native levels) or 'pl'/'pressure' (pressure levels) +## If left blank, each model's built-in default is used. The ECMWF family +## (ERA-5, ERA-5T, HRES) supports both model and pressure levels. +weather_model_levels: + ########## 2. Date ## REQUIRED: TRUE diff --git a/tools/RAiDER/cli/examples/template/template.yaml b/tools/RAiDER/cli/examples/template/template.yaml index 52e8be916..a80bb4997 100644 --- a/tools/RAiDER/cli/examples/template/template.yaml +++ b/tools/RAiDER/cli/examples/template/template.yaml @@ -31,6 +31,15 @@ look_dir: right ## for more details and information on licensing weather_model: +## OPTIONAL: Vertical level representation used by the weather model. +## FORMATS: string, one of 'ml'/'model' (model/native levels) or 'pl'/'pressure' (pressure levels) +## If left blank, each model's built-in default is used (model levels for the +## ECMWF family, GMAO, MERRA-2, NCMR; native levels for HRRR). +## NOTE: Not every model supports every level type. The ECMWF family +## (ERA-5, ERA-5T, HRES) supports both; HRRR supports only native/model levels; +## GMAO, MERRA-2, and NCMR support only model levels. +weather_model_levels: + ########## 2. Date ## REQUIRED: TRUE @@ -70,12 +79,17 @@ time_group: ## 2. Specify a geocoded file, e.g. ARIA GUNW product, from which the AOI will be determined ## 3/4. lat/lon raster files (such as those produced by the ISCE software) ## 5. A comma-delimited file (station_file) containing at least the columns Lat and Lon, and optionally Hgt_m +## station_file_crs declares the height datum of the Hgt_m column: +## 4326 (default) = heights are geoid-referenced (MSL), matching the weather model z-axis +## 4979 = heights are WGS84 ellipsoidal (e.g. GNSS/UNR station lists); they will be +## converted to geoid heights before the weather model is sampled aoi_group: bounding_box: geocoded_file: lat_file: lon_file: station_file: + station_file_crs: ########## 5. Height info diff --git a/tools/RAiDER/cli/raider.py b/tools/RAiDER/cli/raider.py index 0622c4b4d..741d32fd6 100644 --- a/tools/RAiDER/cli/raider.py +++ b/tools/RAiDER/cli/raider.py @@ -126,7 +126,11 @@ def read_run_config_file(path: Path) -> RunConfig: return RunConfig( look_dir=yaml_data['look_dir'].lower(), - weather_model=parse_weather_model(yaml_data['weather_model'], aoi_group.aoi), + weather_model=parse_weather_model( + yaml_data['weather_model'], + aoi_group.aoi, + level_type=yaml_data.get('weather_model_levels'), + ), date_group=parse_dates(DateGroupUnparsed(**yaml_data['date_group'])), time_group=TimeGroup(**yaml_data['time_group']), aoi_group=aoi_group, @@ -267,6 +271,25 @@ def calcDelays(iargs: Optional[Sequence[str]]=None) -> list[Path]: model.set_latlon_bounds(wm_bounds, output_spacing=aoi.get_output_spacing()) + # Batch pre-download for CDS-backed models (ERA5, ERA5T): collect all + # datetimes across the full date_list and issue a single API request. + if hasattr(model, 'batch_fetch'): + _interp = run_config.time_group.interpolate_time or 'none' + _step = model.dtime() if model.dtime() is not None else 6 + _all_times: list[dt.datetime] = [] + for _t in run_config.date_group.date_list: + if _interp == 'center_time': + _all_times.extend(get_nearest_wmtimes(_t, _step)) + elif _interp == 'azimuth_time_grid': + _all_times.extend(get_times_for_azimuth_interpolation(_t, _step)) + else: + _all_times.append(_t) + RAiDER.processWM.batch_download_weather_model( + model, + list(dict.fromkeys(_all_times)), # deduplicate, preserve order + wm_bounds, + ) + wet_paths: list[Path] = [] t: dt.datetime w: str @@ -335,7 +358,8 @@ def calcDelays(iargs: Optional[Sequence[str]]=None) -> list[Path]: if len(wfiles) == 0: logger.error('No weather model data was successfully processed.') - raise NoWeatherModelData('Weather model processing failed for all times') + # raise NoWeatherModelData('Weather model processing failed for all times') + continue # Get the weather model file weather_model_file = getWeatherFile(wfiles, times, t, model._Name, interp_method) @@ -570,6 +594,19 @@ def calcDelaysGUNW(iargs: Optional[list[str]] = None) -> Optional[xr.Dataset]: help='Weather model API KEY [key, password], depending on model.' ) + p.add_argument( + '-l', + '--model-levels', + default=None, + choices=['ml', 'model', 'pl', 'pressure'], + help=( + 'Vertical level representation used by the weather model: model/native ' + "levels ('ml'/'model') or pressure levels ('pl'/'pressure'). If not " + "specified, the model's built-in default is used. Note that not every " + 'model supports every level type.' + ), + ) + p.add_argument( '-interp', '--interpolate-time', @@ -814,41 +851,45 @@ def combine_weather_files(wfiles: list[Path], time: dt.datetime, model: str, int # read the individual datetime datasets datasets = [xr.open_dataset(f) for f in wfiles] - # Pull the datetimes from the datasets - times: list[dt.datetime] = [] - for ds in datasets: - times.append(dt.datetime.strptime(ds.attrs['datetime'], '%Y_%m_%dT%H_%M_%S')) + try: + # Pull the datetimes from the datasets + times: list[dt.datetime] = [] + for ds in datasets: + times.append(dt.datetime.strptime(ds.attrs['datetime'], '%Y_%m_%dT%H_%M_%S')) - if len(times) == 0: - raise NoWeatherModelData() + if len(times) == 0: + raise NoWeatherModelData() - # calculate relative weights of each dataset - if interp_method == 'center_time': - wgts = get_weights_time_interp(times, time) - elif interp_method == 'azimuth_time_grid': - time_grid = get_time_grid_for_aztime_interp(datasets, time, model) - wgts = get_inverse_weights_for_dates(time_grid, times) - else: # interp_method == 'none' - raise ValueError('Interpolating weather files is not available with interpolation method "none"') - - # combine datasets - ds_out = datasets[0] - for var in ['wet', 'hydro', 'wet_total', 'hydro_total']: - ds_out[var] = sum([wgt * ds[var] for (wgt, ds) in zip(wgts, datasets)]) - ds_out.attrs['Date1'] = 0 - ds_out.attrs['Date2'] = 0 - - # Give the weighted combination a new file name - weather_model_file = wfiles[0].parent / ( - wfiles[0].name.split('_')[0] - + '_' - + time.strftime('%Y_%m_%dT%H_%M_%S') - + STYLE[interp_method] - + '_'.join(wfiles[0].name.split('_')[-4:]) - ) + # calculate relative weights of each dataset + if interp_method == 'center_time': + wgts = get_weights_time_interp(times, time) + elif interp_method == 'azimuth_time_grid': + time_grid = get_time_grid_for_aztime_interp(datasets, time, model) + wgts = get_inverse_weights_for_dates(time_grid, times) + else: # interp_method == 'none' + raise ValueError('Interpolating weather files is not available with interpolation method "none"') + + # combine datasets + ds_out = datasets[0] + for var in ['wet', 'hydro', 'wet_total', 'hydro_total']: + ds_out[var] = sum([wgt * ds[var] for (wgt, ds) in zip(wgts, datasets)]) + ds_out.attrs['Date1'] = 0 + ds_out.attrs['Date2'] = 0 + + # Give the weighted combination a new file name + weather_model_file = wfiles[0].parent / ( + wfiles[0].name.split('_')[0] + + '_' + + time.strftime('%Y_%m_%dT%H_%M_%S') + + STYLE[interp_method] + + '_'.join(wfiles[0].name.split('_')[-4:]) + ) - # write the combined results to disk - ds_out.to_netcdf(weather_model_file) + # write the combined results to disk (must happen before closing datasets) + ds_out.to_netcdf(weather_model_file) + finally: + for ds in datasets: + ds.close() return weather_model_file @@ -858,36 +899,40 @@ def combine_files_using_azimuth_time(wfiles, time: dt.datetime, times: list[dt.d # read the individual datetime datasets datasets = [xr.open_dataset(f) for f in wfiles] - # Pull the datetimes from the datasets - times: list[dt.datetime] = [] - for ds in datasets: - times.append(dt.datetime.strptime(ds.attrs['datetime'], '%Y_%m_%dT%H_%M_%S')) - - model = datasets[0].attrs['model_name'] + try: + # Pull the datetimes from the datasets + times: list[dt.datetime] = [] + for ds in datasets: + times.append(dt.datetime.strptime(ds.attrs['datetime'], '%Y_%m_%dT%H_%M_%S')) - time_grid = get_time_grid_for_aztime_interp(datasets, times, time, model) + model = datasets[0].attrs['model_name'] - wgts = get_inverse_weights_for_dates(time_grid, times) + time_grid = get_time_grid_for_aztime_interp(datasets, times, time, model) - # combine datasets - ds_out = datasets[0] - for var in ['wet', 'hydro', 'wet_total', 'hydro_total']: - ds_out[var] = sum([wgt * ds[var] for (wgt, ds) in zip(wgts, datasets)]) - ds_out.attrs['Date1'] = 0 - ds_out.attrs['Date2'] = 0 + wgts = get_inverse_weights_for_dates(time_grid, times) - # Give the weighted combination a new file name - weather_model_file = os.path.join( - os.path.dirname(wfiles[0]), - os.path.basename(wfiles[0]).split('_')[0] - + '_' - + time.strftime('%Y_%m_%dT%H_%M_%S') - + '_timeInterpAziGrid_' - + '_'.join(wfiles[0].split('_')[-4:]), - ) + # combine datasets + ds_out = datasets[0] + for var in ['wet', 'hydro', 'wet_total', 'hydro_total']: + ds_out[var] = sum([wgt * ds[var] for (wgt, ds) in zip(wgts, datasets)]) + ds_out.attrs['Date1'] = 0 + ds_out.attrs['Date2'] = 0 + + # Give the weighted combination a new file name + weather_model_file = os.path.join( + os.path.dirname(wfiles[0]), + os.path.basename(wfiles[0]).split('_')[0] + + '_' + + time.strftime('%Y_%m_%dT%H_%M_%S') + + '_timeInterpAziGrid_' + + '_'.join(wfiles[0].split('_')[-4:]), + ) - # write the combined results to disk - ds_out.to_netcdf(weather_model_file) + # write the combined results to disk (must happen before closing datasets) + ds_out.to_netcdf(weather_model_file) + finally: + for ds in datasets: + ds.close() return weather_model_file diff --git a/tools/RAiDER/cli/types.py b/tools/RAiDER/cli/types.py index 246fa10f5..4d78a60c5 100644 --- a/tools/RAiDER/cli/types.py +++ b/tools/RAiDER/cli/types.py @@ -112,6 +112,11 @@ class AOIGroupUnparsed: lat_file: Optional[str] = None lon_file: Optional[str] = None station_file: Optional[str] = None + # Height datum of the station file's Hgt_m column. 4326 (default) means + # heights are already geoid-referenced (MSL); 4979 means WGS84 ellipsoidal + # heights (e.g. GNSS/UNR station lists), which will be converted to geoid + # heights before sampling the weather model. + station_file_crs: Optional[Union[int, str]] = None geo_cube: Optional[str] = None @dataclasses.dataclass diff --git a/tools/RAiDER/cli/validators.py b/tools/RAiDER/cli/validators.py index 350de1e08..77928fcfa 100755 --- a/tools/RAiDER/cli/validators.py +++ b/tools/RAiDER/cli/validators.py @@ -35,7 +35,7 @@ _BUFFER_SIZE = 0.2 # default buffer size in lat/lon degrees -def parse_weather_model(weather_model_name: str, aoi: AOI) -> WeatherModel: +def parse_weather_model(weather_model_name: str, aoi: AOI, level_type: Optional[str] = None) -> WeatherModel: weather_model_name = weather_model_name.upper().replace('-', '') try: _, Model = get_wm_by_name(weather_model_name) @@ -48,6 +48,16 @@ def parse_weather_model(weather_model_name: str, aoi: AOI) -> WeatherModel: model: WeatherModel = Model() model.checkValidBounds(aoi.bounds()) + # Optionally override the default vertical level representation (pressure vs. + # model levels). If unspecified, the model's built-in default is used. + if level_type is not None: + try: + model.setLevelType(level_type) + except (RuntimeError, NotImplementedError) as e: + raise ValueError( + f'weather_model_levels="{level_type}" is not valid for model {weather_model_name}: {e}' + ) + return model @@ -155,7 +165,11 @@ def get_query_region(aoi_group: AOIGroupUnparsed, height_group: HeightGroupUnpar ) elif aoi_group.station_file is not None: - query = StationFile(aoi_group.station_file, cube_spacing_in_m=cube_spacing_in_m) + query = StationFile( + aoi_group.station_file, + cube_spacing_in_m=cube_spacing_in_m, + crs=aoi_group.station_file_crs if aoi_group.station_file_crs is not None else 4326, + ) elif aoi_group.bounding_box is not None: bbox = parse_bbox(aoi_group.bounding_box) diff --git a/tools/RAiDER/delay.py b/tools/RAiDER/delay.py index 32bae9d45..820fea023 100755 --- a/tools/RAiDER/delay.py +++ b/tools/RAiDER/delay.py @@ -62,8 +62,8 @@ def tropo_delay( """ crs = CRS(out_proj) - # Load CRS from weather model file - with xr.load_dataset(weather_model_file) as ds: + # Load CRS and heights from weather model file (single open) + with xr.open_dataset(weather_model_file) as ds: try: wm_proj = CRS.from_wkt(ds['proj'].attrs['crs_wkt']) except KeyError: @@ -71,9 +71,6 @@ def tropo_delay( "WARNING: I can't find a CRS in the weather model file, so I will assume you are using WGS84" ) wm_proj = CRS.from_epsg(4326) - - # get heights - with xr.load_dataset(weather_model_file) as ds: wm_levels = ds.z.values toa = wm_levels.max() - 1 @@ -137,7 +134,7 @@ def _get_delays_on_cube(datetime: dt.datetime, weather_model_file, wm_proj, aoi, try: aoi.xpts except AttributeError: - with xr.load_dataset(weather_model_file) as ds: + with xr.open_dataset(weather_model_file) as ds: x_spacing = ds.x.diff(dim='x').values.mean() y_spacing = ds.y.diff(dim='y').values.mean() aoi.set_output_spacing(ll_res=np.min([x_spacing, y_spacing])) diff --git a/tools/RAiDER/delayFcns.py b/tools/RAiDER/delayFcns.py index 3bcc3f9e9..1fe0a27d9 100755 --- a/tools/RAiDER/delayFcns.py +++ b/tools/RAiDER/delayFcns.py @@ -28,17 +28,22 @@ def getInterpolators(wm_file: Union[xr.Dataset, Path, str], kind: str='pointwise The interpolator grid is (y, x, z) """ # Get the weather model data - ds = wm_file if isinstance(wm_file, xr.Dataset) else xr.load_dataset(wm_file) - - xs_wm = np.array(ds.variables['x'][:]) - ys_wm = np.array(ds.variables['y'][:]) - zs_wm = np.array(ds.variables['z'][:]) - - wet = ds.variables['wet_total' if kind == 'total' else 'wet'][:] - hydro = ds.variables['hydro_total' if kind == 'total' else 'hydro'][:] - - wet = np.array(wet).transpose(1, 2, 0) - hydro = np.array(hydro).transpose(1, 2, 0) + _close_ds = not isinstance(wm_file, xr.Dataset) + ds = wm_file if isinstance(wm_file, xr.Dataset) else xr.open_dataset(wm_file) + + try: + xs_wm = np.array(ds.variables['x'][:]) + ys_wm = np.array(ds.variables['y'][:]) + zs_wm = np.array(ds.variables['z'][:]) + + wet = ds.variables['wet_total' if kind == 'total' else 'wet'][:] + hydro = ds.variables['hydro_total' if kind == 'total' else 'hydro'][:] + + wet = np.array(wet).transpose(1, 2, 0) + hydro = np.array(hydro).transpose(1, 2, 0) + finally: + if _close_ds: + ds.close() if np.any(np.isnan(wet)) or np.any(np.isnan(hydro)): logger.critical('Weather model contains NaNs!') diff --git a/tools/RAiDER/gnss/downloadGNSSDelays.py b/tools/RAiDER/gnss/downloadGNSSDelays.py index 76b82dc1a..22530da92 100755 --- a/tools/RAiDER/gnss/downloadGNSSDelays.py +++ b/tools/RAiDER/gnss/downloadGNSSDelays.py @@ -79,6 +79,15 @@ def get_stats_by_llh(llhBox=None, baseURL=_UNR_URL): url_igs20 = f'{baseURL}gps_timeseries/IGS20/llh/llh.out' col_names = ['ID', 'Lat', 'Lon', 'Hgt_m'] + # Validate the bounding box convention before any network access so bad + # input fails fast rather than after downloading the station holdings. + W, E = llhBox[2], llhBox[3] + if (W > 180.) or (E > 180.): + raise ValueError( + f"Check input -b:{llhBox} longitudes appear to be in the " + "[0, 360] convention. Expected [-180, 180] convention." + ) + # 1. Fetch IGS20 list try: stat_igs = pd.read_csv(url_igs20, sep=r'\s+', names=col_names) diff --git a/tools/RAiDER/gnss/processDelayFiles.py b/tools/RAiDER/gnss/processDelayFiles.py index bae56ddba..cb5598bc8 100644 --- a/tools/RAiDER/gnss/processDelayFiles.py +++ b/tools/RAiDER/gnss/processDelayFiles.py @@ -240,9 +240,20 @@ 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') @@ -256,12 +267,28 @@ def readZTDFile(filename, col_name='ZTD'): # 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: + data['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." + ) data.rename(columns={col_name: 'ZTD'}, inplace=True) return data @@ -724,6 +751,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 +795,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: diff --git a/tools/RAiDER/llreader.py b/tools/RAiDER/llreader.py index bd1e2d5d4..379dd35a5 100644 --- a/tools/RAiDER/llreader.py +++ b/tools/RAiDER/llreader.py @@ -191,15 +191,138 @@ def set_output_xygrid(self, dst_crs: Union[int, str]=4326) -> None: self.ypts = np.arange(out_snwe[1], out_snwe[0] - out_spacing, -out_spacing) +def _parse_crs(crs: Union[int, str, CRS]) -> CRS: + """Parse an EPSG code (int or str), CRS string, or CRS object into a pyproj CRS.""" + if isinstance(crs, CRS): + return crs + try: + return CRS.from_epsg(crs) + except pyproj.exceptions.CRSError: + return CRS(crs) + + +def _ellipsoidal_to_geometric( + lats: np.ndarray, + lons: np.ndarray, + heights: np.ndarray, + crs: CRS, +) -> np.ndarray: + """Convert heights from ellipsoidal to geometric (geoid-referenced) if needed. + + ERA5 geopotential heights are referenced to the geoid (mean sea level). + GNSS-derived station heights (e.g. from UNR MAGNET in IGS20) are ellipsoidal + heights above the WGS84 ellipsoid. In CONUS the geoid sits ~15–35 m below + the ellipsoid, so sampling the RAiDER cube with uncorrected GNSS heights + places the integration surface too low by that amount, causing a ~5–9 mm + positive ZTD bias. + + Conversion targets EPSG:9707 (WGS 84 + EGM96 height). PROJ requires the + EGM96 geoid grid, which ships with the ``proj-data`` conda package. If the + grid is absent locally the function transparently enables the PROJ CDN for a + one-time download and then restores the previous network state. If the + download also fails, input heights are returned unchanged with a warning. + + Args: + lats: station latitudes in degrees + lons: station longitudes in degrees + heights: station heights in the coordinate system described by ``crs`` + crs: pyproj CRS describing the input height datum + + Returns: + heights above the EGM96 geoid (~MSL / orthometric), + or the input heights unchanged if the CRS has no 3D ellipsoidal vertical. + """ + # A 2D CRS (e.g. EPSG:4326) carries no vertical datum — return unchanged. + if len(crs.axis_info) < 3: + return heights + + _TARGET = CRS.from_epsg(9707) # WGS 84 + EGM96 height + + def _build_and_apply() -> tuple: + """Return (transformer, converted_heights).""" + from pyproj import Transformer + t = Transformer.from_crs(crs, _TARGET, always_xy=True) + _, _, h = t.transform(lons, lats, heights) + return t, h + + try: + t, h_geoid = _build_and_apply() + + # PROJ silently falls back to a noop when the EGM96 grid is missing. + # Detect this by inspecting the WKT pipeline string. + # to_proj4() returns '+proj=noop' when PROJ fell back to a no-op + # because the EGM96 grid is absent. A real geoid transform returns + # None (it is too complex to express as a PROJ4 string). + if t.to_proj4() == '+proj=noop': + logger.info( + 'EGM96 geoid grid not found locally; attempting download from ' + 'the PROJ CDN. Install the proj-data conda package for fully ' + 'offline use.' + ) + _was_enabled = pyproj.network.is_network_enabled() + pyproj.network.set_network_enabled(True) + try: + t, h_geoid = _build_and_apply() + if t.to_proj4() == '+proj=noop': + logger.warning( + 'EGM96 grid unavailable (no local data and CDN ' + 'unreachable); ellipsoidal heights used unchanged.' + ) + return heights + finally: + pyproj.network.set_network_enabled(_was_enabled) + + logger.debug( + 'Converted ellipsoidal heights to EGM96 geoid heights; ' + 'mean geoid undulation applied: %.2f m', + float(np.nanmean(h_geoid - heights)), + ) + return h_geoid + + except Exception as exc: + logger.warning( + 'Ellipsoidal-to-geoid height conversion failed (%s); ' + 'using input heights unchanged. ' + 'Ensure proj-data grids are installed for accurate results.', + exc, + ) + return heights + + class StationFile(AOI): - """Use a .csv file containing at least Lat, Lon, and optionally Hgt_m columns.""" + """Use a .csv file containing at least Lat, Lon, and optionally Hgt_m columns. + + By default heights are assumed to be already geoid-referenced (MSL), matching + the ERA5 z-axis convention (``crs=4326``). Pass ``crs=4979`` when the file + contains WGS84 ellipsoidal heights (e.g. from UNR MAGNET / IGS20) so that + they are automatically converted to geoid heights before the weather model + cube is sampled. + """ - def __init__(self, station_file, demFile=None, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None: + def __init__( + self, + station_file: Union[str, Path], + demFile: Optional[Union[str, Path]] = None, + cube_spacing_in_m: Optional[float] = None, + output_directory: Union[str, Path] = Path.cwd(), + crs: Union[int, str, CRS] = 4326, + ) -> None: super().__init__(cube_spacing_in_m, output_directory) self._filename = station_file self._demfile = demFile self._bounding_box = bounds_from_csv(station_file) self._type = 'station_file' + self._crs = _parse_crs(crs) + + def update_crs(self, new_crs: Union[int, str, CRS]) -> None: + """Update the CRS describing the height datum of the station file. + + Args: + new_crs: EPSG code (int or str), CRS string, or pyproj CRS object. + ``4326`` (default) — heights are geoid-referenced (MSL). + ``4979`` — heights are WGS84 ellipsoidal (e.g. from GNSS/UNR). + """ + self._crs = _parse_crs(new_crs) def readLL(self) -> tuple[np.ndarray, np.ndarray]: """Read the station lat/lons from the csv file.""" @@ -207,12 +330,21 @@ def readLL(self) -> tuple[np.ndarray, np.ndarray]: return df['Lat'].to_numpy(), df['Lon'].to_numpy() def readZ(self): - """Read the station heights from the file, or download a DEM if not present.""" + """Read station heights, converting to geoid heights if the CRS is 3D ellipsoidal. + + When constructed with ``crs=4979`` (or any 3D CRS with an ellipsoidal + vertical), heights are converted from WGS84 ellipsoidal to EGM96 geoid + heights so they are consistent with the ERA5 z-axis reference. + DEM-derived heights are already MSL and are never converted. + """ df = pd.read_csv(self._filename).drop_duplicates(subset=['Lat', 'Lon']) if 'Hgt_m' in df.columns: - return df['Hgt_m'].values + heights = df['Hgt_m'].values + return _ellipsoidal_to_geometric( + df['Lat'].values, df['Lon'].values, heights, self._crs + ) else: - # Download the DEM + # Download the DEM (DEM heights are MSL; no conversion needed) from RAiDER.dem import download_dem from RAiDER.interpolator import interpolateDEM diff --git a/tools/RAiDER/models/ecmwf.py b/tools/RAiDER/models/ecmwf.py index 144c61952..401ca1f65 100755 --- a/tools/RAiDER/models/ecmwf.py +++ b/tools/RAiDER/models/ecmwf.py @@ -1,5 +1,4 @@ import datetime as dt -import shutil import tempfile from pathlib import Path @@ -132,12 +131,22 @@ def _get_from_cds( if not corrected_DT == acqTime: logger.warning('Rounded given datetime from %s to %s', acqTime, corrected_DT) + if self._model_level_type == 'pl': + # Pressure-level downloads share the batch writer so that all + # pl files have the same layout and units (z in geopotential m) + self._batch_get_from_cds([(corrected_DT, out_path)], lat_min, lat_max, lon_min, lon_max) + return + + param = ['lnsp', 'z', 'q', 't'] + dataset = 'reanalysis-era5-complete' + with tempfile.TemporaryDirectory() as temp_dir_str: temp_dir = Path(temp_dir_str) - out_path_lnsp_z = temp_dir / f'{out_path.stem}_lnsp_z' - out_path_t_q = temp_dir / f'{out_path.stem}_t_q' + out_path_combined = temp_dir / f'{out_path.stem}_combined' # Developed from https://confluence.ecmwf.int/display/CKB/How+to+download+ERA5 + # All four variables are fetched in a single CDS request to halve queue wait time. + # lnsp and z are surface-only fields; t and q span all model levels. params = { 'class': 'ea', 'expver': '1', @@ -153,54 +162,188 @@ def _get_from_cds( 'area': [lat_max, lon_min, lat_min, lon_max], 'grid': [0.25, 0.25], 'format': 'netcdf', + 'param': param, } - # Make two separate requests: one for lnsp and z, and the other for t and q. - params['param'] = ['lnsp', 'z'] - c.retrieve('reanalysis-era5-complete', params, out_path_lnsp_z) - params['param'] = ['q', 't'] - c.retrieve('reanalysis-era5-complete', params, out_path_t_q) - - # RAiDER requires z data for all levels, but ERA-5 only provides it for - # the surface level. z can be computed for all levels using lnsp, t, and - # q, so that is what we will do. - # We will use the t/q dataset as a base to make a full dataset with - # lnsp, t, and q, and full z data. - shutil.copy(out_path_t_q, out_path) - - with xr.open_dataset(out_path_lnsp_z) as ds_lnsp_z, xr.open_dataset(out_path_t_q) as ds_t_q: - # Compute full z + c.retrieve(dataset, params, out_path_combined) + + with xr.open_dataset(out_path_combined) as ds: + # ERA-5 only provides z at the surface level; compute it at all + # model levels from lnsp, t, and q via the hypsometric equation. + # .squeeze() removes the size-1 time (and level, for lnsp/z) dims, + # since calcgeoh expects (level, lat, lon) or (lat, lon) arrays. z_full, _, _ = util.calcgeoh( - # .squeeze(): data comes in with dimensions: - # (valid time, model level, latitude, longitude), - # and we need it in: - # (model level, latitude, longitude), - # and there is always exactly one valid time - # (i.e., shape is always (1, ..., ..., ...)). - lnsp=ds_lnsp_z['lnsp'].values.squeeze(), - z_surface=ds_lnsp_z['z'].values.squeeze(), - t=ds_t_q['t'].values.squeeze(), - q=ds_t_q['q'].values.squeeze(), + lnsp=ds['lnsp'].values.squeeze(), + z_surface=ds['z'].values.squeeze(), + t=ds['t'].values.squeeze(), + q=ds['q'].values.squeeze(), a=self._a, b=self._b, num_levels=self._levels, R_d=self._R_d, ) - # Add the full z cube to the output dataset. - ds_out = ds_t_q.assign( + # Replace the surface-only z with the full model-level z cube, + # broadcast to match t's (time, level, lat, lon) dimensions. + ds_out = ds.assign( z=xr.Variable( - # Copy over t's dimensions as z's. - # Could also have used q; all three are the same. - dims=ds_t_q['t'].dims, - # Present z_full as though it were wrapped in an array to - # match the shape of the rest of the data. + dims=ds['t'].dims, data=np.broadcast_to(z_full, (1, *z_full.shape)), ), - # To fit the shape of the rest of the dataset, this is NaN on - # every level but the first (the surface). - lnsp=ds_lnsp_z['lnsp'], ) ds_out.to_netcdf(out_path) + def _batch_get_from_cds( + self, + times_and_paths: list[tuple[dt.datetime, Path]], + lat_min: float, + lat_max: float, + lon_min: float, + lon_max: float, + ) -> None: + """Download multiple ERA5 time steps in a single CDS request and split into per-datetime files.""" + import cdsapi + + c = cdsapi.Client(verify=1) + + if c.url == 'https://cds.climate.copernicus.eu/api/v2': + logger.warning( + 'Old CDS API configuration detected: ECMWF released a breaking change in late 2024 that expired all ' + 'existing credentials. This run may fail with a 404 HTTP error, in which case you may have to ' + 'regenerate your CDS API credentials at https://cds.climate.copernicus.eu/how-to-api.' + ) + + # Build unique date and time lists (CDS takes the Cartesian product) + seen_dates: dict[str, None] = {} + seen_times: dict[str, None] = {} + for corrected_dt, _ in times_and_paths: + seen_dates[corrected_dt.strftime('%Y-%m-%d')] = None + seen_times[corrected_dt.strftime('%H:%M')] = None + date_str = '/'.join(seen_dates) + time_str = '/'.join(seen_times) + + base_params = { + 'class': 'ea', + 'expver': '1', + 'levelist': 'all', + 'levtype': self._model_level_type, + 'stream': 'oper', + 'type': 'an', + 'date': date_str, + 'time': time_str, + 'step': '0', + 'area': [lat_max, lon_min, lat_min, lon_max], + 'grid': [0.25, 0.25], + 'format': 'netcdf', + } + + if self._model_level_type == 'ml': + param = ['lnsp', 'z', 'q', 't'] + dataset = 'reanalysis-era5-complete' + else: + param = ['z', 't', 'q'] + dataset = 'reanalysis-era5-pressure-levels' + + with tempfile.TemporaryDirectory() as temp_dir_str: + temp_dir = Path(temp_dir_str) + surface_file = temp_dir / 'batch_surface.nc' + ml_file = temp_dir / 'batch_ml.nc' + + if self._model_level_type == 'pl': + # MARS for reanalysis-era5-pressure-levels allows only one date per month + # per request; multiple same-month dates cause "Duplicate value for month". + # The CDS API v3 also auto-splits a list by month before sending to MARS, + # so batching is not possible. Issue one request per datetime. + for corrected_dt, out_path in times_and_paths: + pl_params = { + 'levelist': 'all', + 'levtype': 'pl', + 'date': corrected_dt.strftime('%Y-%m-%d'), + 'time': corrected_dt.strftime('%H:%M'), + 'area': [lat_max, lon_min, lat_min, lon_max], + 'data_format': 'netcdf', + 'param': param, + } + c.retrieve(dataset, pl_params, ml_file) + + with xr.open_dataset(ml_file) as ds: + tc = 'valid_time' if 'valid_time' in ds.coords else 'time' + target = corrected_dt.replace(tzinfo=None) + ds_slice = ds.isel({tc: 0}) + + z_v = (ds_slice['z'].values.squeeze() / self._g0).transpose(1, 2, 0)[np.newaxis] + t_v = ds_slice['t'].values.squeeze().transpose(1, 2, 0)[np.newaxis] + q_v = ds_slice['q'].values.squeeze().transpose(1, 2, 0)[np.newaxis] + + ds_out = xr.Dataset( + { + 'z': xr.Variable((tc, 'latitude', 'longitude', 'pressure_level'), z_v), + 't': xr.Variable((tc, 'latitude', 'longitude', 'pressure_level'), t_v), + 'q': xr.Variable((tc, 'latitude', 'longitude', 'pressure_level'), q_v), + }, + coords={ + tc: np.array([target], dtype='datetime64[ns]'), + 'pressure_level': ds_slice['pressure_level'].values, + 'latitude': ds_slice['latitude'].values, + 'longitude': ds_slice['longitude'].values, + }, + ) + ds_out.to_netcdf(out_path) + else: + # Model level requests only include lnsp and z at the surface, so we have to make two requests to get + # the full model-level z cube. This is a bit less efficient, but still much better than making + # separate requests for each datetime. All four variables are fetched in a single CDS request to + # halve queue wait time. + # lnsp and z are surface fields; t and q span all model levels. + # For multi-datetime batch requests the CDS server returns mixed-level + # requests as separate files, so request the two groups independently. + c.retrieve('reanalysis-era5-complete', {**base_params, 'param': ['lnsp', 'z']}, surface_file) + c.retrieve('reanalysis-era5-complete', {**base_params, 'param': ['t', 'q']}, ml_file) + + with xr.open_dataset(surface_file) as ds_surface, xr.open_dataset(ml_file) as ds_ml: + # CDS API uses 'valid_time' in newer versions, 'time' in older ones + tc_s = 'valid_time' if 'valid_time' in ds_surface.coords else 'time' + tc_m = 'valid_time' if 'valid_time' in ds_ml.coords else 'time' + + for corrected_dt, out_path in times_and_paths: + # Strip timezone: numpy datetime64 coordinates are timezone-naive + target = corrected_dt.replace(tzinfo=None) + surf_slice = ds_surface.sel({tc_s: target}) + ml_slice = ds_ml.sel({tc_m: target}) + + lnsp_v = surf_slice['lnsp'].values.squeeze() # (nlat, nlon) + z_sfc_v = surf_slice['z'].values.squeeze() # (nlat, nlon) + t_v = ml_slice['t'].values.squeeze() # (nlev, nlat, nlon) + q_v = ml_slice['q'].values.squeeze() # (nlev, nlat, nlon) + + z_full, _, _ = util.calcgeoh( + lnsp=lnsp_v, + z_surface=z_sfc_v, + t=t_v, + q=q_v, + a=self._a, + b=self._b, + num_levels=self._levels, + R_d=self._R_d, + ) + + # Write a single-time file matching the structure _makeDataCubes + # expects: t/q/z as (time, model_level, lat, lon) and lnsp as + # (time, sfc_level, lat, lon) so that [0,0] indexing gives (lat,lon). + ds_out = xr.Dataset( + { + 't': xr.Variable((tc_m, 'model_level', 'latitude', 'longitude'), t_v[np.newaxis]), + 'q': xr.Variable((tc_m, 'model_level', 'latitude', 'longitude'), q_v[np.newaxis]), + 'z': xr.Variable((tc_m, 'model_level', 'latitude', 'longitude'), z_full[np.newaxis]), + 'lnsp': xr.Variable((tc_m, 'sfc_level', 'latitude', 'longitude'), lnsp_v[np.newaxis, np.newaxis]), + }, + coords={ + tc_m: np.array([target], dtype='datetime64[ns]'), + 'model_level': ml_slice['model_level'].values, + 'latitude': ml_slice['latitude'].values, + 'longitude': ml_slice['longitude'].values, + }, + ) + ds_out.to_netcdf(out_path) + def _download_ecmwf(self, lat_min, lat_max, lat_step, lon_min, lon_max, lon_step, time, out: Path) -> None: """Used for HRES.""" from ecmwfapi import ECMWFService @@ -242,7 +385,7 @@ def _load_model_level(self, filename, *args, **kwargs) -> None: # read data from netcdf file lats, lons, _, _, t, q, lnsp, z = self._makeDataCubes(Path(filename)) - # ECMWF appears to give me this backwards + # data ordering if lats[0] > lats[1]: z: FloatArray3D = z[::-1] lnsp: FloatArray2D = lnsp[::-1] @@ -289,29 +432,48 @@ def _load_model_level(self, filename, *args, **kwargs) -> None: self._zs = np.flip(self._zs, axis=2) def _load_pressure_level(self, filename) -> None: - with xr.open_dataset(filename) as block: - # Pull the data - z = np.squeeze(block['z'].values) - t = np.squeeze(block['t'].values) - q = np.squeeze(block['q'].values) - lats = np.squeeze(block['latitude'].values) - lons = np.squeeze(block['longitude'].values) - levels = np.squeeze(block['level'].values) * 100 - - z = np.flip(z, axis=1) - - # ECMWF appears to give me this backwards + with xr.open_dataset(filename) as ds: + # Drop the singleton time dimension (name varies by CDS API version) + for time_dim in ('valid_time', 'time'): + if time_dim in ds.dims: + ds = ds.isel({time_dim: 0}) + break + + lev_dim = 'pressure_level' if 'pressure_level' in ds.dims else 'level' + # Normalize dimension order by name so the file's on-disk layout + # doesn't matter + ds = ds.transpose('latitude', 'longitude', lev_dim) + + z = ds['z'].values.astype(np.float64) + t = ds['t'].values + q = ds['q'].values + lats = ds['latitude'].values + lons = ds['longitude'].values + levels = ds[lev_dim].values * 100 # hPa -> Pa + + # Files written by _batch_get_from_cds store z as geopotential height + # (already divided by g0); raw CDS files store geopotential (m^2/s^2). + # Geopotential height tops out below ~50 km, so use that to distinguish. + if np.nanmax(z) > 100_000: + z = z / self._g0 + + # Reorder axes (consistently across all cubes) so lats and lons are + # ascending and levels go surface -> TOA if lats[0] > lats[1]: z = z[::-1] - t = t[:, ::-1] - q = q[:, ::-1] + t = t[::-1] + q = q[::-1] lats = lats[::-1] - # Lons is usually ok, but we'll throw in a check to be safe if lons[0] > lons[1]: + z = z[:, ::-1] + t = t[:, ::-1] + q = q[:, ::-1] + lons = lons[::-1] + if levels[0] < levels[-1]: z = z[..., ::-1] t = t[..., ::-1] q = q[..., ::-1] - lons = lons[::-1] + levels = levels[::-1] # pyproj gets fussy if the latitude is wrong, plus our # interpolator isn't clever enough to pick up on the fact that # they are the same @@ -320,27 +482,17 @@ def _load_pressure_level(self, filename) -> None: self._t = t self._q = q - geo_hgt = (z / self._g0).transpose(1, 2, 0) - # re-assign lons, lats to match heights self._lons, self._lats = np.meshgrid(lons, lats) # correct heights for latitude - self._get_heights(self._lats, geo_hgt) + self._get_heights(self._lats, z) self._p = np.broadcast_to(levels[np.newaxis, np.newaxis, :], self._zs.shape) - # Re-structure from (heights, lats, lons) to (lons, lats, heights) - self._t = self._t.transpose(1, 2, 0) - self._q = self._q.transpose(1, 2, 0) self._ys = self._lats.copy() self._xs = self._lons.copy() - # flip z to go from surface to toa - self._p = np.flip(self._p, axis=2) - self._t = np.flip(self._t, axis=2) - self._q = np.flip(self._q, axis=2) - def _makeDataCubes( self, path: Path, diff --git a/tools/RAiDER/models/era5.py b/tools/RAiDER/models/era5.py index 610367ec5..e3934436b 100755 --- a/tools/RAiDER/models/era5.py +++ b/tools/RAiDER/models/era5.py @@ -4,6 +4,8 @@ from dateutil.relativedelta import relativedelta from pyproj import CRS +from RAiDER import utilFcns as util +from RAiDER.logger import logger from RAiDER.models.ecmwf import ECMWF @@ -31,7 +33,6 @@ def __init__(self) -> None: # Availability lag time in days self._lag_time = relativedelta(months=lag_time) - # Default, need to change to ml self.setLevelType('ml') def _fetch(self, out: Path) -> None: @@ -43,6 +44,22 @@ def _fetch(self, out: Path) -> None: # execute the search at ECMWF self._get_from_cds(lat_min, lat_max, lon_min, lon_max, time, out) + def batch_fetch(self, times_and_paths: list[tuple[dt.datetime, Path]]) -> None: + """Download multiple ERA5 datetimes in a single CDS API call.""" + rounded: list[tuple[dt.datetime, Path]] = [] + for acqTime, out_path in times_and_paths: + corrected_DT = util.round_date(acqTime, dt.timedelta(hours=self._time_res)) + if corrected_DT != acqTime: + logger.warning('Rounded given datetime from %s to %s', acqTime, corrected_DT) + if not out_path.exists(): + rounded.append((corrected_DT, out_path)) + + if not rounded: + return + + lat_min, lat_max, lon_min, lon_max = self._ll_bounds + self._batch_get_from_cds(rounded, lat_min, lat_max, lon_min, lon_max) + def load_weather(self, f=None, *args, **kwargs) -> None: """Load either pressure or model level data.""" f = self.files[0] if f is None else f diff --git a/tools/RAiDER/models/gmao.py b/tools/RAiDER/models/gmao.py index f1176ffe1..6ad44fac0 100755 --- a/tools/RAiDER/models/gmao.py +++ b/tools/RAiDER/models/gmao.py @@ -58,6 +58,9 @@ def __init__(self) -> None: # Projection self._proj = CRS.from_epsg(4326) + def __model_levels__(self) -> None: + self._zlevels = np.flipud(LEVELS_137_HEIGHTS) + def _fetch(self, out: Path) -> None: """Fetch weather model data from GMAO.""" # calculate the array indices for slicing the GMAO variable arrays diff --git a/tools/RAiDER/models/hres.py b/tools/RAiDER/models/hres.py index 9a7a5bdf2..c8ac96d33 100755 --- a/tools/RAiDER/models/hres.py +++ b/tools/RAiDER/models/hres.py @@ -72,7 +72,7 @@ def load_weather(self, f=None) -> None: self.update_a_b() self._load_model_level(f) elif self._model_level_type == 'pl': - self._load_pressure_levels(f) + self._load_pressure_level(f) def _fetch(self, out: Path) -> None: """Fetch a weather model from ECMWF.""" diff --git a/tools/RAiDER/models/merra2.py b/tools/RAiDER/models/merra2.py index 8416e95e9..357fdfe48 100755 --- a/tools/RAiDER/models/merra2.py +++ b/tools/RAiDER/models/merra2.py @@ -68,6 +68,9 @@ def __init__(self) -> None: # Projection self._proj = CRS.from_epsg(4326) + def __model_levels__(self) -> None: + self._zlevels = np.flipud(LEVELS_137_HEIGHTS) + def _fetch(self, out: Path) -> None: """Fetch weather model data from GMAO: note we only extract the lat/lon bounds for this weather model; fetching data is not needed here as we don't actually download any data using OpenDAP.""" time = self._time @@ -136,13 +139,13 @@ def load_weather(self, f=None, *args, **kwargs) -> None: def _load_model_level(self, filename) -> None: """Get the variables from the GMAO link using OpenDAP.""" # adding the import here should become absolute when transition to netcdf - ds = xr.load_dataset(filename) - lons = ds['longitude'].values - lats = ds['latitude'].values - h = ds['h'].values - q = ds['q'].values - p = ds['p'].values - t = ds['t'].values + with xr.open_dataset(filename) as ds: + lons = ds['longitude'].values + lats = ds['latitude'].values + h = ds['h'].values + q = ds['q'].values + p = ds['p'].values + t = ds['t'].values # Re-structure everything from (heights, lats, lons) to (lons, lats, heights) p = np.transpose(p) diff --git a/tools/RAiDER/models/ncmr.py b/tools/RAiDER/models/ncmr.py index 58338f398..b2e78616d 100755 --- a/tools/RAiDER/models/ncmr.py +++ b/tools/RAiDER/models/ncmr.py @@ -63,6 +63,9 @@ def __init__(self) -> None: # Projection self._proj = CRS.from_epsg(4326) + def __model_levels__(self) -> None: + self._zlevels = np.flipud(LEVELS_137_HEIGHTS) + def _fetch(self, out: Path) -> None: """ Fetch weather model data from NCMR: note we only extract the lat/lon bounds for this weather model; diff --git a/tools/RAiDER/models/weatherModel.py b/tools/RAiDER/models/weatherModel.py index e9412e9f5..69696ef36 100755 --- a/tools/RAiDER/models/weatherModel.py +++ b/tools/RAiDER/models/weatherModel.py @@ -323,8 +323,33 @@ def checkTime(self, time: dt.datetime) -> None: if time > dt.datetime.now(dt.timezone.utc) - self._lag_time: raise DatetimeOutsideRange(self.Model(), time) + def __model_levels__(self) -> None: + """Configure the model to use native/model levels. + + Models that support model levels should override this to populate the + relevant level arrays (e.g. ``self._levels``, ``self._zlevels``). + """ + raise NotImplementedError(f'Weather model {self.Model()} does not support model levels') + + def __pressure_levels__(self) -> None: + """Configure the model to use pressure levels. + + Models that support pressure levels should override this to populate the + relevant level arrays (e.g. ``self._levels``, ``self._zlevels``). + """ + raise NotImplementedError(f'Weather model {self.Model()} does not support pressure levels') + def setLevelType(self, levelType: str) -> None: - """Set the level type to model levels or pressure levels.""" + """Set the level type to model levels or pressure levels. + + Accepts the internal codes ('ml', 'nat' for model/native levels; + 'pl', 'prs' for pressure levels) as well as the user-friendly aliases + 'model' and 'pressure'. + """ + # Map user-friendly aliases onto the internal codes. + aliases = {'model': 'ml', 'pressure': 'pl'} + levelType = aliases.get(levelType.lower(), levelType) + if levelType in 'ml pl nat prs'.split(): self._model_level_type = levelType else: @@ -373,8 +398,19 @@ def _get_wet_refractivity(self) -> None: self._wet_refractivity = self._k2 * self._e / self._t + self._k3 * self._e / self._t**2 def _get_hydro_refractivity(self) -> None: - """Calculate the hydrostatic delay from pressure and temperature.""" - self._hydrostatic_refractivity = self._k1 * self._p / self._t + """Calculate the hydrostatic refractivity from pressure, temperature, and e. + + Hydrostatic refractivity is k1 * P / Tv = k1 * Rd * rho (rho = total air + density, Davis et al. 1985), written here via the exact identity + k1 * P / Tv = k1 * (P - (1 - Rd/Rv) * e) / T. Using virtual temperature + (not T) is required for consistency with the k2' = k2 - k1*Rd/Rv wet + coefficient (0.233 K/Pa) used by all models; pairing k1*P/T with k2' + double-counts part of the water-vapor contribution (~2% of the wet + delay). This split matches the GNSS ZHD/ZWD convention. + """ + self._hydrostatic_refractivity = ( + self._k1 * (self._p - (1 - self._R_d / self._R_v) * self._e) / self._t + ) def getWetRefractivity(self) -> np.ndarray: """Returns the data cube of refractivity.""" @@ -410,11 +446,11 @@ def _getZTD(self) -> None: wet = self.getWetRefractivity() hydro = self.getHydroRefractivity() - # Get the integrated ZTD - wet_total, hydro_total = np.zeros(wet.shape), np.zeros(hydro.shape) - for level in range(wet.shape[2]): - wet_total[..., level] = 1e-6 * np_trapezoid(wet[..., level:], x=self._zs[level:], axis=2) - hydro_total[..., level] = 1e-6 * np_trapezoid(hydro[..., level:], x=self._zs[level:], axis=2) + # Get the integrated ZTD. Layers are integrated assuming exponential + # variation of refractivity with height; plain trapezoid integration + # overestimates ZTD by ~1 cm on the coarse fixed z-levels (convex N). + wet_total = 1e-6 * util.cumulative_integral_from_top(wet, self._zs) + hydro_total = 1e-6 * util.cumulative_integral_from_top(hydro, self._zs) self._hydrostatic_ztd = hydro_total self._wet_ztd = wet_total @@ -451,7 +487,7 @@ def bbox(self) -> Union[list, tuple, np.ndarray]: if not Path.exists(Path(path_weather_model)): raise ValueError('Need to save cropped weather model as netcdf') - with xr.load_dataset(path_weather_model) as ds: + with xr.open_dataset(path_weather_model) as ds: try: xmin, xmax = ds.x.min(), ds.x.max() ymin, ymax = ds.y.min(), ds.y.max() @@ -872,8 +908,5 @@ def checkContainment_raw( return weather_model_box.contains(input_box) - elif weather_model_box.contains(world_box): - return True - else: - return False + return weather_model_box.contains(input_box) \ No newline at end of file diff --git a/tools/RAiDER/processWM.py b/tools/RAiDER/processWM.py index d8ab4b80d..6528d9e9b 100755 --- a/tools/RAiDER/processWM.py +++ b/tools/RAiDER/processWM.py @@ -52,7 +52,6 @@ def prepareWeatherModel( # get the path to the less processed weather model file path_wm_raw = make_raw_weather_data_filename(wmLoc, weather_model.Model(), time) - # get the path to the more processed (cropped) weather model file path_wm_crop = weather_model.out_file(wmLoc) @@ -135,6 +134,30 @@ def prepareWeatherModel( return f +def batch_download_weather_model( + weather_model, + times: list, + ll_bounds, + force_download: bool = False, +) -> None: + """Pre-download all ERA5/ERA5T datetimes in a single CDS API call. + + Writes per-datetime raw files to disk so that subsequent prepareWeatherModel + calls find them already present and skip the download step. + """ + wmLoc = weather_model.get_wmLoc() + missing: list[tuple] = [] + for t in times: + path_wm_raw = Path(make_raw_weather_data_filename(wmLoc, weather_model.Model(), t)) + if not force_download and path_wm_raw.exists() and checkContainment_raw(path_wm_raw, ll_bounds): + continue + os.makedirs(path_wm_raw.parent, exist_ok=True) + missing.append((t, path_wm_raw)) + + if missing: + weather_model.batch_fetch(missing) + + def _weather_model_debug(los, lats, lons, ll_bounds, weather_model, wmLoc, time, out, download_only) -> None: """RaiderWeatherModelDebug main function.""" logger.debug('Starting to run the weather model calculation with debugging plots') diff --git a/tools/RAiDER/utilFcns.py b/tools/RAiDER/utilFcns.py index 30f22d081..5ee2cab5f 100644 --- a/tools/RAiDER/utilFcns.py +++ b/tools/RAiDER/utilFcns.py @@ -984,3 +984,41 @@ def parse_crs(proj: CRSLike) -> CRS: np_trapezoid = np.trapezoid else: np_trapezoid = np.trapz + + +def cumulative_integral_from_top(ns: np.ndarray, zs: np.ndarray) -> np.ndarray: + """Cumulatively integrate refractivity from each height level to the column top. + + Refractivity decays quasi-exponentially with height, so the trapezoid rule + systematically overestimates each layer integral (the chord lies above a + convex curve); on the coarse fixed z-levels used for pressure-level models + this accumulates to ~1 cm of zenith delay. Instead, each layer is + integrated assuming exponential variation between its endpoints + (the logarithmic-mean rule): + + int_z0^z1 N dz = (z1 - z0) * (N0 - N1) / ln(N0 / N1) + + which is exact for N(z) = N0 * exp(-(z - z0)/H). Layers where the + exponential model is undefined (non-positive or nearly equal endpoint + values) fall back to the trapezoid rule. + + Args: + ns: refractivity, shape (..., nz), levels ascending in height along the last axis + zs: 1-D array of heights (m), length nz, ascending + + Returns: + ndarray of shape (..., nz): integral from each level to the top level + (the top level is 0 by construction). + """ + n0 = ns[..., :-1].astype(np.float64) + n1 = ns[..., 1:].astype(np.float64) + dz = np.diff(np.asarray(zs, dtype=np.float64)) + trap = 0.5 * (n0 + n1) * dz + with np.errstate(divide='ignore', invalid='ignore'): + lnr = np.log(n0 / n1) + expo = dz * (n0 - n1) / lnr + use_expo = (n0 > 0) & (n1 > 0) & (np.abs(lnr) > 1e-6) & np.isfinite(expo) + seg = np.where(use_expo, expo, trap) + out = np.zeros(ns.shape, dtype=np.float64) + out[..., :-1] = np.cumsum(seg[..., ::-1], axis=-1)[..., ::-1] + return out