Fix MERRA-2 access

CHANGELOG.md

@@ -6,6 +6,10 @@ 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).
 
+## Latest Changes
+### Fixed
+* [620](https://github.com/dbekaert/RAiDER/issues/620) - Fix MERRA-2 access because API was updated
+
 
 ## [0.4.7]
 ### Fixed

test/__init__.py

@@ -16,7 +16,7 @@
 WM_DIR   = os.path.join(TEST_DIR, 'weather_files')
 ORB_DIR  = os.path.join(TEST_DIR, 'orbit_files')
 
-WM = 'GMAO'
+WM = 'MERRA2'
 
 @contextmanager
 def pushd(dir):

test/test_GUNW.py

@@ -277,7 +277,7 @@ def test_azimuth_timing_interp_against_center_time_interp(weather_model_name: st
             assert np.nanmax(abs_diff_mm) < 1
 
 
-@pytest.mark.parametrize('weather_model_name', ['GMAO', 'HRRR', 'HRES', 'ERA5', 'ERA5'])
+@pytest.mark.parametrize('weather_model_name', ['HRRR', 'HRES', 'ERA5', 'ERA5T'])
 def test_check_weather_model_availability(test_gunw_path_factory, weather_model_name, mocker):
     # Should be True for all weather models
     # S1-GUNW-D-R-071-tops-20200130_20200124-135156-34956N_32979N-PP-913f-v2_0_4.nc
@@ -288,12 +288,12 @@ def test_check_weather_model_availability(test_gunw_path_factory, weather_model_
     mocker.patch("RAiDER.aria.prepFromGUNW.get_acq_time_from_slc_id", side_effect=[pd.Timestamp('2015-01-01'),
                                                                                    pd.Timestamp('2014-01-01')])
     cond = check_weather_model_availability(test_gunw_path, weather_model_name)
-    if weather_model_name in ['HRRR', 'GMAO']:
+    if weather_model_name in ['HRRR', 'MERRA2']:
         cond = not cond
     assert cond
 
 
-@pytest.mark.parametrize('weather_model_name', ['GMAO', 'HRRR'])
+@pytest.mark.parametrize('weather_model_name', ['HRRR'])
 def test_check_weather_model_availability_over_alaska(test_gunw_path_factory, weather_model_name, mocker):
     # Should be True for all weather models
     # S1-GUNW-D-R-059-tops-20230320_20220418-180300-00179W_00051N-PP-c92e-v2_0_6.nc
@@ -309,7 +309,7 @@ def test_check_weather_model_availability_over_alaska(test_gunw_path_factory, we
     assert cond
 
 
-@pytest.mark.parametrize('weather_model_name', ['HRRR', 'GMAO'])
+@pytest.mark.parametrize('weather_model_name', ['HRRR'])
 @pytest.mark.parametrize('location', ['california-t71', 'alaska'])
 def test_weather_model_availability_integration_using_valid_range(location,
                                                                   test_gunw_path_factory,

tools/RAiDER/aria/prepFromGUNW.py

@@ -26,7 +26,7 @@
 from RAiDER.s1_orbits import ensure_orbit_credentials
 
 ## cube spacing in degrees for each model
-DCT_POSTING = {'HRRR': 0.05, 'HRES': 0.10, 'GMAO': 0.10, 'ERA5': 0.10, 'ERA5T': 0.10}
+DCT_POSTING = {'HRRR': 0.05, 'HRES': 0.10, 'GMAO': 0.10, 'ERA5': 0.10, 'ERA5T': 0.10, 'MERRA2': 0.1}
 
 
 def _get_acq_time_from_gunw_id(gunw_id: str, reference_or_secondary: str) -> datetime:
@@ -95,7 +95,7 @@ def check_weather_model_availability(gunw_path: str,
     ----------
     gunw_path : str
     weather_model_name : str
-        Should be one of 'HRRR', 'HRES', 'ERA5', 'ERA5T', 'GMAO'.
+        Should be one of 'HRRR', 'HRES', 'ERA5', 'ERA5T', 'GMAO', 'MERRA2'.
     Returns
     -------
     bool:

tools/RAiDER/models/__init__.py

@@ -3,5 +3,7 @@
 from .gmao import GMAO
 from .hres import HRES
 from .hrrr import HRRR, HRRRAK
+from .merra2 import MERRA2
+from .ncmr import NCMR
 
-__all__ = ['HRRR', 'HRRRAK', 'GMAO', 'ERA5', 'ERA5T', 'HRES']
+__all__ = ['HRRR', 'HRRRAK', 'GMAO', 'ERA5', 'ERA5T', 'HRES', 'MERRA2']

tools/RAiDER/models/allowed.py

@@ -3,5 +3,7 @@
     'ERA5T',
     'HRRR',
     'GMAO',
-    'HRES'
+    'HRES',
+    'MERRA2',
+    'NCMR',
 ]

tools/RAiDER/models/gmao.py

@@ -9,7 +9,7 @@
 
 from RAiDER.models.weatherModel import WeatherModel, TIME_RES
 from RAiDER.logger import logger
-from RAiDER.utilFcns import writeWeatherVars2NETCDF4, round_date, requests_retry_session
+from RAiDER.utilFcns import writeWeatherVarsXarray, round_date, requests_retry_session
 from RAiDER.models.model_levels import (
     LEVELS_137_HEIGHTS,
 )
@@ -145,7 +145,7 @@ def _fetch(self, out):
 
         try:
             # Note that lat/lon gets written twice for GMAO because they are the same as y/x
-            writeWeatherVars2NETCDF4(self, lats, lons, h, qv, p, t, outName=out)
+            writeWeatherVarsXarray(lat, lon, h, q, p, t, dt, crs, outName=None, NoDataValue=None, chunk=(1, 91, 144))
         except Exception:
             logger.exception("Unable to save weathermodel to file")
 

tools/RAiDER/models/merra2.py

@@ -1,4 +1,6 @@
 import os
+import xarray
+
 import datetime as dt
 import numpy as np
 import pydap.cas.urs
@@ -8,7 +10,7 @@
 
 from RAiDER.models.weatherModel import WeatherModel
 from RAiDER.logger import logger
-from RAiDER.utilFcns import writeWeatherVars2NETCDF4, read_EarthData_loginInfo
+from RAiDER.utilFcns import writeWeatherVarsXarray, read_EarthData_loginInfo
 from RAiDER.models.model_levels import (
     LEVELS_137_HEIGHTS,
 )
@@ -94,6 +96,8 @@ def _fetch(self, out):
             self._lon_res
         )
 
+        lon, lat = np.meshgrid(lons, lats)
+
         if time.year < 1992:
             url_sub = 100
         elif time.year < 2001:
@@ -107,39 +111,22 @@ def _fetch(self, out):
         DT = time - T0
         time_ind = int(DT.total_seconds() / 3600.0 / 3.0)
 
-        ml_min = 0
-        ml_max = 71
-
         # Earthdata credentials
         earthdata_usr, earthdata_pwd = read_EarthData_loginInfo(EARTHDATA_RC)
 
         # open the dataset and pull the data
-        url = 'https://goldsmr5.gesdisc.eosdis.nasa.gov:443/opendap/MERRA2/M2I3NVASM.5.12.4/' + time.strftime('%Y/%m') + '/MERRA2_' + str(url_sub) + '.inst3_3d_asm_Nv.' + time.strftime('%Y%m%d') + '.nc4'
+        url = 'https://goldsmr5.gesdisc.eosdis.nasa.gov/opendap/MERRA2/M2T3NVASM.5.12.4/' + time.strftime('%Y/%m') + '/MERRA2_' + str(url_sub) + '.tavg3_3d_asm_Nv.' + time.strftime('%Y%m%d') + '.nc4'
+
         session = pydap.cas.urs.setup_session(earthdata_usr, earthdata_pwd, check_url=url)
-        ds = pydap.client.open_url(url, session=session)
+        stream = pydap.client.open_url(url, session=session)
 
-        ############# The MERRA-2 server changes the pydap data retrieval format frequently between these two formats; so better to retain both of them rather than only using either one of them #############
-        try:
-            q = ds['QV'].data[0][time_ind, ml_min:(ml_max + 1), lat_min_ind:(lat_max_ind + 1), lon_min_ind:(lon_max_ind + 1)][0]
-            p = ds['PL'].data[0][time_ind, ml_min:(ml_max + 1), lat_min_ind:(lat_max_ind + 1), lon_min_ind:(lon_max_ind + 1)][0]
-            t = ds['T'].data[0][time_ind, ml_min:(ml_max + 1), lat_min_ind:(lat_max_ind + 1), lon_min_ind:(lon_max_ind + 1)][0]
-            h = ds['H'].data[0][time_ind, ml_min:(ml_max + 1), lat_min_ind:(lat_max_ind + 1), lon_min_ind:(lon_max_ind + 1)][0]
-        except IndexError:
-            q = ds['QV'].data[time_ind, ml_min:(ml_max + 1), lat_min_ind:(lat_max_ind + 1), lon_min_ind:(lon_max_ind + 1)][0]
-            p = ds['PL'].data[time_ind, ml_min:(ml_max + 1), lat_min_ind:(lat_max_ind + 1), lon_min_ind:(lon_max_ind + 1)][0]
-            t = ds['T'].data[time_ind, ml_min:(ml_max + 1), lat_min_ind:(lat_max_ind + 1), lon_min_ind:(lon_max_ind + 1)][0]
-            h = ds['H'].data[time_ind, ml_min:(ml_max + 1), lat_min_ind:(lat_max_ind + 1), lon_min_ind:(lon_max_ind + 1)][0]
-        except AttributeError:
-            q = ds['QV'][time_ind, ml_min:(ml_max + 1), lat_min_ind:(lat_max_ind + 1), lon_min_ind:(lon_max_ind + 1)][0]
-            p = ds['PL'][time_ind, ml_min:(ml_max + 1), lat_min_ind:(lat_max_ind + 1), lon_min_ind:(lon_max_ind + 1)][0]
-            t = ds['T'][time_ind, ml_min:(ml_max + 1), lat_min_ind:(lat_max_ind + 1), lon_min_ind:(lon_max_ind + 1)][0]
-            h = ds['H'][time_ind, ml_min:(ml_max + 1), lat_min_ind:(lat_max_ind + 1), lon_min_ind:(lon_max_ind + 1)][0]
-        except BaseException:
-            logger.exception("MERRA-2: Unable to read weathermodel data")
-        ########################################################################################################################
+        q = stream['QV'][0,:,lat_min_ind:lat_max_ind + 1, lon_min_ind:lon_max_ind + 1].data.squeeze()
+        p = stream['PL'][0,:,lat_min_ind:lat_max_ind + 1, lon_min_ind:lon_max_ind + 1].data.squeeze()
+        t = stream['T'][0,:,lat_min_ind:lat_max_ind + 1, lon_min_ind:lon_max_ind + 1].data.squeeze()
+        h = stream['H'][0,:,lat_min_ind:lat_max_ind + 1, lon_min_ind:lon_max_ind + 1].data.squeeze()
 
         try:
-            writeWeatherVars2NETCDF4(self, lats, lons, h, q, p, t, outName=out)
+            writeWeatherVarsXarray(lat, lon, h, q, p, t, time, self._proj, outName=out)
         except Exception as e:
             logger.debug(e)
             logger.exception("MERRA-2: Unable to save weathermodel to file")
@@ -161,51 +148,39 @@ def _load_model_level(self, filename):
         '''
 
         # adding the import here should become absolute when transition to netcdf
-        from netCDF4 import Dataset
-        with Dataset(filename, mode='r') as f:
-            lons = np.array(f.variables['x'][:])
-            lats = np.array(f.variables['y'][:])
-            h = np.array(f.variables['H'][:])
-            q = np.array(f.variables['QV'][:])
-            p = np.array(f.variables['PL'][:])
-            t = np.array(f.variables['T'][:])
-
-        # restructure the 3-D lat/lon/h in regular grid
-        _lons = np.broadcast_to(lons[np.newaxis, np.newaxis, :], t.shape)
-        _lats = np.broadcast_to(lats[np.newaxis, :, np.newaxis], t.shape)
+        ds = xarray.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
 
         # Re-structure everything from (heights, lats, lons) to (lons, lats, heights)
         p = np.transpose(p)
         q = np.transpose(q)
         t = np.transpose(t)
         h = np.transpose(h)
-        _lats = np.transpose(_lats)
-        _lons = np.transpose(_lons)
 
         # check this
         # data cube format should be lats,lons,heights
         p = p.swapaxes(0, 1)
         q = q.swapaxes(0, 1)
         t = t.swapaxes(0, 1)
         h = h.swapaxes(0, 1)
-        _lats = _lats.swapaxes(0, 1)
-        _lons = _lons.swapaxes(0, 1)
 
         # For some reason z is opposite the others
         p = np.flip(p, axis=2)
         q = np.flip(q, axis=2)
         t = np.flip(t, axis=2)
         h = np.flip(h, axis=2)
-        _lats = np.flip(_lats, axis=2)
-        _lons = np.flip(_lons, axis=2)
 
         # assign the regular-grid (lat/lon/h) variables
-
         self._p = p
         self._q = q
         self._t = t
-        self._lats = _lats
-        self._lons = _lons
-        self._xs = _lons
-        self._ys = _lats
+        self._lats = lats
+        self._lons = lons
+        self._xs = lons
+        self._ys = lats
         self._zs = h

tools/RAiDER/models/ncmr.py

@@ -13,9 +13,9 @@
 from RAiDER.models.weatherModel import WeatherModel, TIME_RES
 from RAiDER.logger import logger
 from RAiDER.utilFcns import (
-    writeWeatherVars2NETCDF4,
     read_NCMR_loginInfo,
-    show_progress
+    show_progress,
+    writeWeatherVarsXarray,
 )
 from RAiDER.models.model_levels import (
     LEVELS_137_HEIGHTS,
@@ -174,7 +174,7 @@ def _download_ncmr_file(self, out, date_time, bounding_box):
         ########################################################################################################################
 
         try:
-            writeWeatherVars2NETCDF4(self, lats, lons, hgt, q, p, t, outName=out)
+            writeWeatherVarsXarray(lats, lons, hgt, q, p, t, self._time, self._proj, outName=out)
         except Exception:
             logger.exception("Unable to save weathermodel to file")
 

tools/RAiDER/models/weatherModel.py

@@ -20,7 +20,7 @@
 from RAiDER.models import plotWeather as plots, weatherModel
 from RAiDER.models.customExceptions import *
 from RAiDER.utilFcns import (
-    robmax, robmin, write2NETCDF4core, calcgeoh, transform_coords, clip_bbox
+    robmax, robmin, writeWeatherVarsXarray, calcgeoh, transform_coords, clip_bbox
 )
 
 TIME_RES = {'GMAO': 3,