Fix GMAO

CHANGELOG.md

@@ -9,6 +9,7 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 ## [Unreleased]
 ### Changed
 * [743](https://github.com/dbekaert/RAiDER/pull/743) - Switched from HTTPS to DAP4 for retrieving MERRA2 data, and suppressed a warning for using DAP4 for GMAO data where doing so is not possible.
+* []() - Changed GMAO and MERRA2 fetch implementations to use xarray, eliminating pydap as a dependency.
 
 ### Fixed
 * [741](https://github.com/dbekaert/RAiDER/pull/741) - Updated mamba setup commands in CircleCI for mamba 2.0.0.

docs/WeatherModels.md

@@ -71,7 +71,7 @@ ECMWF requires a license agreement to be able to access, download, and use their
 
 ## 4. NASA weather models (GMAO, MERRA2)
 
-1. The Global Modeling and Assimilation Office (__[GMAO](https://www.nccs.nasa.gov/services/data-collections/coupled-products/geos5-forecast#:~:text=The%20Global%20Modeling%20and%20Assimilation,near%2Dreal%2Dtime%20production.)__) at NASA generates reanalysis weather models. GMAO datasets can also be accessed without a license agreement through the pyDAP interface implemented in RAiDER. GMAO has a horizontal grid spacing of approximately 33 km, and its projection is EPSG code 4326 (WGS-84). 
+1. The Global Modeling and Assimilation Office (__[GMAO](https://www.nccs.nasa.gov/services/data-collections/coupled-products/geos5-forecast#:~:text=The%20Global%20Modeling%20and%20Assimilation,near%2Dreal%2Dtime%20production.)__) at NASA generates reanalysis weather models. GMAO datasets can be accessed without a license agreement through the OPeNDAP interface implemented in RAiDER. GMAO has a horizontal grid spacing of approximately 33 km, and its projection is EPSG code 4326 (WGS-84). 
 
 2. The Modern-Era Retrospective analysis for Research and Applications, Version 2 (__[MERRA-2](https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/#:~:text=MERRA%2D2%20is%20the%20first,(say)%20Greenland%20and%20Antarctica.)__) provides data beginning in 1980. MERRA-2 is also produced by NASA and has a spatial resolution of about 50 km and a global projection (EPSG 4326, WGS-84).  
 
@@ -88,12 +88,4 @@ Reference: __[The Modern-Era Retrospective Analysis for Research and Application
 
     **Note**: the username and password represent the user's username and password.
 
-4. Add the application `NASA GESDISC DATA ARCHIVE` by clicking on the `Applications->Authorized Apps` on the menu after logging into your Earthdata profile, and then scrolling down to the application `NASA GESDISC DATA ARCHIVE` to approve it. _This seems not required for GMAO for now, but recommended to do so for all OpenDAP-based weather models._
-5. Install the OpenDAP using pip: 
-
-         pip install pydap==3.2.1
-
-
-    **Note**: this step has been included in the conda install of RAiDER, thus can be omitted if one uses the recommended conda install of RAiDER
-
-    **Note**: PyDAP v3.2.1 is required for now (thus specified in the above pip install command) because the latest v3.2.2 (as of now) has a known [bug](https://colab.research.google.com/drive/1f_ss1Oa3VzgAOd_p8sgekdnLVE5NW6s5) in accessing and slicing the GMAO data. This bug is expected to be fixed in newer versions of PyDAP.
+4. Add the application `NASA GESDISC DATA ARCHIVE` by clicking on the `Applications->Authorized Apps` on the menu after logging into your Earthdata profile, and then scrolling down to the application `NASA GESDISC DATA ARCHIVE` to approve it. _This seems not required for GMAO for now, but recommended to do so for all OPeNDAP-based weather models._

test/test_HRRR_ztd.py

@@ -4,8 +4,9 @@
 import xarray as xr
 from RAiDER.cli.raider import calcDelays
 
-def test_scenario_1(tmp_path, data_for_hrrr_ztd, mocker):
-    SCENARIO_DIR = TEST_DIR / "scenario_1"
+
+def test_hrrr_ztd(tmp_path, data_for_hrrr_ztd, mocker) -> None:
+    SCENARIO_DIR = TEST_DIR / 'scenario_1'
     test_path = SCENARIO_DIR / 'raider_example_1.yaml'
     mocker.patch('RAiDER.processWM.prepareWeatherModel',
                     side_effect=[str(data_for_hrrr_ztd)])

test/test_dem.py

@@ -7,7 +7,7 @@
 def test_download_dem_1():
     SCENARIO_1 = TEST_DIR / "scenario_4"
     hts, meta = download_dem(
-        dem_path=SCENARIO_1 / 'warpedDEM.dem', 
+        dem_path=SCENARIO_1 / 'warpedDEM.rdr', 
         overwrite=False
     )
     assert hts.shape == (45,226)

test/test_downloaders.py

@@ -13,26 +13,29 @@
 from test import random_string
 
 
-DATETIME = dt.datetime(2020, 1, 1, 0, 0, 0).replace(tzinfo=dt.timezone(offset=dt.timedelta()))
 BOUNDS = np.array([10, 10.2, -72, -72])
+DATETIME = dt.datetime(2020, 1, 1, 0, 0, 0).replace(tzinfo=dt.timezone(offset=dt.timedelta()))
+DATETIME_GMAO_OLD = dt.datetime(2017, 1, 1, 0, 0, 0).replace(tzinfo=dt.timezone(offset=dt.timedelta()))
 
 
 @pytest.mark.long
 @pytest.mark.parametrize(
-    'name,Model',
+    'Model,time',
     [
-        ('ERA5', ERA5),
-        ('ERA5T', ERA5T),
-        pytest.param('HRES', HRES, marks=pytest.mark.skip),  # Paid access
-        ('GMAO', GMAO),
-        ('MERRA2', MERRA2),
+        pytest.param(ERA5, DATETIME, id='ERA5'),
+        pytest.param(ERA5T, DATETIME, id='ERA5T'),
+        pytest.param(HRES, DATETIME, id='HRES', marks=pytest.mark.skip),  # Paid access
+        # HRRR: see test_weather_model.py
+        pytest.param(GMAO, DATETIME, id='GMAO new'),
+        pytest.param(GMAO, DATETIME_GMAO_OLD, id='GMAO old'),
+        pytest.param(MERRA2, DATETIME, id='MERRA2'),
     ],
 )
-def test_downloader(tmp_path: Path, name: str, Model: Type[WeatherModel]) -> None:
-    out_path = tmp_path / f'test_{name}.nc'
+def test_downloader(tmp_path: Path, Model: Type[WeatherModel], time: dt.datetime) -> None:
     wm = Model()
+    out_path = tmp_path / f'test_{wm._Name}.nc'
     wm.set_latlon_bounds(BOUNDS)
-    wm.fetch(out_path, DATETIME)
+    wm.fetch(out_path, time)
 
 
 @pytest.mark.long

test/test_scenario_1.py

@@ -0,0 +1,62 @@
+import datetime as dt
+from pathlib import Path
+from typing import Type
+
+import numpy as np
+import pytest
+
+from RAiDER.delay import tropo_delay
+from RAiDER.llreader import RasterRDR
+from RAiDER.losreader import Zenith
+from RAiDER.models import ERA5, ERA5T, GMAO, MERRA2
+from RAiDER.models.weatherModel import WeatherModel
+from RAiDER.processWM import prepareWeatherModel
+from RAiDER.utilFcns import rio_open
+from test import TEST_DIR
+
+
+SCENARIO_DIR = TEST_DIR / 'scenario_1'
+WM_LOC = SCENARIO_DIR / 'weather_files'
+
+_RTOL = 1e-2
+DATETIME = dt.datetime(2018, 7, 1, 0, 0)
+DATETIME_GMAO_OLD = dt.datetime(2017, 1, 1, 0, 0, 0)
+
+
+@pytest.mark.long
+@pytest.mark.parametrize(
+    'Model,time',
+    (
+        pytest.param(ERA5, DATETIME, id='ERA5'),
+        pytest.param(ERA5T, DATETIME, id='ERA5T'),
+        pytest.param(GMAO, DATETIME, id='GMAO new'),
+        pytest.param(GMAO, DATETIME_GMAO_OLD, id='GMAO old'),
+        pytest.param(MERRA2, DATETIME, id='MERRA2'),
+    ),
+)
+def test_tropo_delay(tmp_path: Path, Model: Type[WeatherModel], time: dt.datetime) -> None:
+    """
+    Scenario:
+    1: Small area, Zenith delay.
+    """
+    WM_LOC.mkdir(exist_ok=True)
+
+    los = Zenith()
+    lat_path = str(SCENARIO_DIR / 'geom/lat.dat')
+    lon_path = str(SCENARIO_DIR / 'geom/lon.dat')
+    hgt_file = str(TEST_DIR / 'test_geom/warpedDEM.rdr')
+    aoi = RasterRDR(lat_path, lon_path, hgt_file=hgt_file, output_directory=str(tmp_path))
+
+    wm = Model()
+    wm.set_latlon_bounds(aoi.bounds())
+    wm.setTime(time)
+    wm.set_wmLoc(str(WM_LOC))
+    wm_file_path = prepareWeatherModel(wm, time, aoi.bounds())
+    wet, hydro = tropo_delay(time, wm_file_path, aoi, los, zref=20_000)
+
+    # load the true delay
+    true_wet, _ = rio_open(SCENARIO_DIR / f'{wm._Name}/wet.envi', userNDV=0.0)
+    true_hydro, _ = rio_open(SCENARIO_DIR / f'{wm._Name}/hydro.envi', userNDV=0.0)
+
+    assert np.nanmax(np.abs((wet - true_wet) / true_wet)) < _RTOL
+    assert np.nanmax(np.abs((hydro - true_hydro) / true_hydro)) < _RTOL