diff --git a/CHANGELOG.md b/CHANGELOG.md index bed54c979..f276ab1af 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/docs/WeatherModels.md b/docs/WeatherModels.md index 8a48feb64..b21e110ec 100644 --- a/docs/WeatherModels.md +++ b/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._ diff --git a/test/_scenario_1.py b/test/_scenario_1.py deleted file mode 100755 index 4f296f1e3..000000000 --- a/test/_scenario_1.py +++ /dev/null @@ -1,144 +0,0 @@ -import datetime -import os -import pytest - -import numpy as np - -from test import TEST_DIR, pushd - -from RAiDER.losreader import Zenith -from RAiDER.delay import main -from RAiDER.utilFcns import rio_open -from RAiDER.checkArgs import makeDelayFileNames -from RAiDER.cli.validators import get_wm_by_name - -SCENARIO_DIR = os.path.join(TEST_DIR, "scenario_1") -_RTOL = 1e-2 - - -@pytest.mark.long -def test_tropo_delay_ERAI(tmp_path): - ''' - Scenario: - 1: Small area, ERAI, Zenith delay - ''' - core_test_tropo_delay(tmp_path, modelName="ERAI") - - -@pytest.mark.long -def test_tropo_delay_ERA5(tmp_path): - ''' - Scenario: - 1: Small area, ERA5, Zenith delay - ''' - core_test_tropo_delay(tmp_path, modelName="ERA5") - - -@pytest.mark.long -def test_tropo_delay_ERA5T(tmp_path): - ''' - Scenario: - 1: Small area, ERA5T, Zenith delay - ''' - core_test_tropo_delay(tmp_path, modelName="ERA5T") - - -@pytest.mark.long -def test_tropo_delay_MERRA2(tmp_path): - ''' - Scenario: - 1: Small area, MERRA2, Zenith delay - ''' - core_test_tropo_delay(tmp_path, modelName="MERRA2") - - -@pytest.mark.skip(reason="NCMR keeps hanging") -def test_tropo_delay_NCMR(tmp_path): - ''' - Scenario: - 1: Small area, NCMR, Zenith delay - ''' - core_test_tropo_delay(tmp_path, modelName="NCMR") - - -@pytest.mark.long -def test_tropo_delay_GMAO(tmp_path): - ''' - Scenario: - 1: Small area, GMAO, Zenith delay - ''' - core_test_tropo_delay(tmp_path, modelName="GMAO") - - -def core_test_tropo_delay(tmp_path, modelName): - ''' - Scenario: - 1: Small area, Zenith delay - ''' - lats = rio_open(os.path.join( - SCENARIO_DIR, 'geom', 'lat.dat' - )) - lons = rio_open(os.path.join( - SCENARIO_DIR, 'geom', 'lon.dat' - )) - - if modelName == 'ERAI': - time = datetime.datetime(2018, 1, 3, 23, 0) - elif modelName == 'NCMR': - time = datetime.datetime(2018, 7, 1, 0, 0) - else: - time = datetime.datetime(2020, 1, 3, 23, 0) - - wmLoc = os.path.join(SCENARIO_DIR, 'weather_files') - if not os.path.exists(wmLoc): - os.mkdir(wmLoc) - - _, model_obj = get_wm_by_name(modelName) - wet_file, hydro_file = makeDelayFileNames( - time, Zenith, "envi", modelName, tmp_path - ) - - with pushd(tmp_path): - # packing the dictionairy - args = {} - args['los'] = Zenith - args['lats'] = lats - args['lons'] = lons - args['ll_bounds'] = (15.75, 18.25, -103.24, -99.75) - args['heights'] = ("dem", os.path.join(TEST_DIR, "test_geom", "warpedDEM.dem")) - args['pnts_file'] = 'lat_query_points.h5' - args['flag'] = "files" - args['weather_model'] = {"type": model_obj(), "files": None, "name": modelName} - args['wmLoc'] = wmLoc - args['zref'] = 20000. - args['outformat'] = "envi" - args['times'] = time - args['out'] = tmp_path - args['download_only'] = False - args['wetFilenames'] = wet_file - args['hydroFilenames'] = hydro_file - args['verbose'] = True - - (_, _) = main(args) - - # get the results - wet = rio_open(wet_file) - hydro = rio_open(hydro_file) - true_wet = rio_open( - os.path.join( - SCENARIO_DIR, - modelName + "/wet.envi" - ), - userNDV=0. - ) - true_hydro = rio_open( - os.path.join( - SCENARIO_DIR, - modelName + "/hydro.envi" - ), - userNDV=0. - ) - - # get the true delay from the weather model - assert np.nanmax(np.abs((wet - true_wet) / true_wet)) < _RTOL - assert np.nanmax(np.abs((hydro - true_hydro) / true_hydro)) < _RTOL diff --git a/test/scenario_0/geom/warpedDEM.dem b/test/scenario_0/geom/warpedDEM.rdr similarity index 100% rename from test/scenario_0/geom/warpedDEM.dem rename to test/scenario_0/geom/warpedDEM.rdr diff --git a/test/scenario_4/warpedDEM.dem b/test/scenario_4/warpedDEM.rdr similarity index 100% rename from test/scenario_4/warpedDEM.dem rename to test/scenario_4/warpedDEM.rdr diff --git a/test/scenario_4/warpedDEM.dem.aux.xml b/test/scenario_4/warpedDEM.rdr.aux.xml similarity index 100% rename from test/scenario_4/warpedDEM.dem.aux.xml rename to test/scenario_4/warpedDEM.rdr.aux.xml diff --git a/test/scenario_7/warpedDEM.dem b/test/scenario_7/warpedDEM.rdr similarity index 100% rename from test/scenario_7/warpedDEM.dem rename to test/scenario_7/warpedDEM.rdr diff --git a/test/scenario_7/warpedDEM.dem.aux.xml b/test/scenario_7/warpedDEM.rdr.aux.xml similarity index 100% rename from test/scenario_7/warpedDEM.dem.aux.xml rename to test/scenario_7/warpedDEM.rdr.aux.xml diff --git a/test/scenario_8/warpedDEM.dem b/test/scenario_8/warpedDEM.rdr similarity index 100% rename from test/scenario_8/warpedDEM.dem rename to test/scenario_8/warpedDEM.rdr diff --git a/test/scenario_8/warpedDEM.dem.aux.xml b/test/scenario_8/warpedDEM.rdr.aux.xml similarity index 100% rename from test/scenario_8/warpedDEM.dem.aux.xml rename to test/scenario_8/warpedDEM.rdr.aux.xml diff --git a/test/test_HRRR_ztd.py b/test/test_HRRR_ztd.py index 4777de4ad..cd7b38674 100644 --- a/test/test_HRRR_ztd.py +++ b/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)]) diff --git a/test/test_dem.py b/test/test_dem.py index 97fc06e44..c1575a737 100644 --- a/test/test_dem.py +++ b/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) diff --git a/test/test_downloaders.py b/test/test_downloaders.py index 5eaa2d6f2..debfc53b1 100644 --- a/test/test_downloaders.py +++ b/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 diff --git a/test/test_geom/warpedDEM.dem b/test/test_geom/warpedDEM.rdr similarity index 100% rename from test/test_geom/warpedDEM.dem rename to test/test_geom/warpedDEM.rdr diff --git a/test/test_scenario_1.py b/test/test_scenario_1.py new file mode 100644 index 000000000..ea34563f9 --- /dev/null +++ b/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 diff --git a/test/test_scenario_4.py b/test/test_scenario_4.py index f1ed25efd..e388df88d 100644 --- a/test/test_scenario_4.py +++ b/test/test_scenario_4.py @@ -1,67 +1,77 @@ -import datetime -import os -import pytest -import xarray +import datetime as dt +from pathlib import Path +from typing import Type import numpy as np +import pytest +import xarray from pyproj import CRS -from RAiDER.delay import tropo_delay, _get_delays_on_cube +from RAiDER.delay import _get_delays_on_cube, tropo_delay from RAiDER.llreader import RasterRDR from RAiDER.losreader import Zenith -from RAiDER.processWM import prepareWeatherModel +from RAiDER.models.era5 import ERA5 +from RAiDER.models.era5t import ERA5T +from RAiDER.models.gmao import GMAO from RAiDER.models.merra2 import MERRA2 - +from RAiDER.models.ncmr import NCMR +from RAiDER.models.weatherModel import WeatherModel +from RAiDER.processWM import prepareWeatherModel from test import TEST_DIR, pushd -SCENARIO_DIR = os.path.join(TEST_DIR, "scenario_4") + + +SCENARIO_DIR = TEST_DIR / 'scenario_4' + +DATETIME = dt.datetime(2020, 1, 1) + @pytest.mark.long -def test_aoi_without_xpts(tmp_path): +@pytest.mark.parametrize('Model', (GMAO, NCMR, MERRA2, ERA5T, ERA5)) +def test_aoi_without_xpts(tmp_path: Path, Model: Type[WeatherModel]) -> None: with pushd(tmp_path): los = Zenith() - latfile = os.path.join(SCENARIO_DIR, 'lat.rdr') - lonfile = os.path.join(SCENARIO_DIR, 'lon.rdr') - hgtfile = os.path.join(SCENARIO_DIR, 'hgt.rdr') + latfile = SCENARIO_DIR / 'lat.rdr' + lonfile = SCENARIO_DIR / 'lon.rdr' + hgtfile = SCENARIO_DIR / 'hgt.rdr' aoi = RasterRDR(latfile, lonfile, hgtfile) - dt = datetime.datetime(2020,1,1) - wm = MERRA2() + wm = Model() wm.set_latlon_bounds(aoi.bounds()) - wm.setTime(dt) - f = prepareWeatherModel(wm, dt, aoi.bounds()) - zen_wet, zen_hydro = tropo_delay(dt, f, aoi, los) + wm.setTime(DATETIME) + f = prepareWeatherModel(wm, DATETIME, aoi.bounds()) + zen_wet, zen_hydro = tropo_delay(DATETIME, f, aoi, los) assert len(zen_wet.shape) == 2 assert np.sum(np.isnan(zen_wet)) < np.prod(zen_wet.shape) assert np.nanmean(zen_wet) > 0 assert np.nanmean(zen_hydro) > 0 + @pytest.mark.long -def test_get_delays_on_cube(tmp_path): +@pytest.mark.parametrize('Model', (GMAO, NCMR, MERRA2, ERA5T, ERA5)) +def test_get_delays_on_cube(tmp_path: Path, Model: Type[WeatherModel]) -> None: with pushd(tmp_path): los = Zenith() - latfile = os.path.join(SCENARIO_DIR, 'lat.rdr') - lonfile = os.path.join(SCENARIO_DIR, 'lon.rdr') - hgtfile = os.path.join(SCENARIO_DIR, 'hgt.rdr') + latfile = SCENARIO_DIR / 'lat.rdr' + lonfile = SCENARIO_DIR / 'lon.rdr' + hgtfile = SCENARIO_DIR / 'hgt.rdr' aoi = RasterRDR(latfile, lonfile, hgtfile) - dt = datetime.datetime(2020,1,1) - wm = MERRA2() + wm = Model() wm.set_latlon_bounds(aoi.bounds()) - wm.setTime(dt) - f = prepareWeatherModel(wm, dt, aoi.bounds()) + wm.setTime(DATETIME) + f = prepareWeatherModel(wm, DATETIME, aoi.bounds()) with xarray.load_dataset(f) as ds: - wm_levels = ds.z.values + wm_levels = ds['z'].values wm_proj = CRS.from_wkt(ds['proj'].attrs['crs_wkt']) - + zref = 10000 with pytest.raises(AttributeError): aoi.xpts - - ds = _get_delays_on_cube(dt, f, wm_proj, aoi, wm_levels, los, wm_proj, zref) + + ds = _get_delays_on_cube(DATETIME, f, wm_proj, aoi, wm_levels, los, wm_proj, zref) assert len(ds.x) > 0 - assert ds.hydro.mean() > 0 - + assert ds['hydro'].mean() > 0 diff --git a/test/test_util.py b/test/test_util.py index a0a0e5084..230333bf9 100644 --- a/test/test_util.py +++ b/test/test_util.py @@ -33,7 +33,7 @@ transform_bbox, unproject, writeArrayToRaster, - writeWeatherVarsXarray, + write_weather_vars_to_ds, ) from test import TEST_DIR, pushd @@ -916,8 +916,8 @@ def test_UTM_to_WGS84_empty_input(): assert lat.size == 0 -# Test writeWeatherVarsXarray -def test_writeWeatherVarsXarray(tmp_path): +# Test write_weather_vars_to_ds +def test_write_weather_vars_to_ds(tmp_path): """Test writing weather variables to an xarray dataset and NetCDF file.""" # Mock inputs lat = np.random.rand(91, 144) * 180 - 90 # Random latitudes between -90 and 90 @@ -938,7 +938,7 @@ def test_writeWeatherVarsXarray(tmp_path): out_path = tmp_path / "test_output.nc" # Call the function - writeWeatherVarsXarray(lat, lon, h, q, p, t, datetime_value, crs, out_path) + write_weather_vars_to_ds(lat, lon, h, q, p, t, datetime_value, crs, out_path) # Check that the file was created assert out_path.exists() diff --git a/test/test_weather_model.py b/test/test_weather_model.py index 276a6a07d..a992b3790 100644 --- a/test/test_weather_model.py +++ b/test/test_weather_model.py @@ -403,7 +403,7 @@ def test_ztd(model: MockWeatherModel) -> None: def test_get_bounds_indices() -> None: """Test bounds indices.""" - snwe = [-10, 10, -10, 10] + snwe = (-10, 10, -10, 10) ll = np.arange(-20, 20) lats, lons = np.meshgrid(ll, ll) xmin, xmax, ymin, ymax = get_bounds_indices(snwe, lats, lons) @@ -415,7 +415,7 @@ def test_get_bounds_indices() -> None: def test_get_bounds_indices_2() -> None: """Test bounds indices.""" - snwe = [-10, 10, 170, -170] + snwe = (-10, 10, 170, -170) l = np.arange(-20, 20) l2 = (((np.arange(160, 200) + 180) % 360) - 180) lats, lons = np.meshgrid(l, l2) @@ -425,7 +425,7 @@ def test_get_bounds_indices_2() -> None: def test_get_bounds_indices_2b() -> None: """Test bounds indices.""" - snwe = [-10, 10, 170, 190] + snwe = (-10, 10, 170, 190) l = np.arange(-20, 20) l2 = np.arange(160, 200) lats, lons = np.meshgrid(l, l2) @@ -438,7 +438,7 @@ def test_get_bounds_indices_2b() -> None: def test_get_bounds_indices_3() -> None: """Test bounds indices""" - snwe = [-10, 10, -10, 10] + snwe = (-10, 10, -10, 10) l = np.arange(-20, 20) l2 = (((np.arange(160, 200) + 180) % 360) - 180) lats, lons = np.meshgrid(l, l2) @@ -448,7 +448,7 @@ def test_get_bounds_indices_3() -> None: def test_get_bounds_indices_4() -> None: """Test bounds_indices.""" - snwe = [55, 60, 175, 185] + snwe = (55, 60, 175, 185) l = np.arange(55, 60, 1) l2 = np.arange(175, 185, 1) lats, lons = np.meshgrid(l, l2) diff --git a/tools/RAiDER/cli/validators.py b/tools/RAiDER/cli/validators.py index dda62ec54..86e21164e 100755 --- a/tools/RAiDER/cli/validators.py +++ b/tools/RAiDER/cli/validators.py @@ -150,7 +150,7 @@ def get_query_region(aoi_group: AOIGroupUnparsed, height_group: HeightGroupUnpar raise ValueError('A lon_file must be specified if a lat_file is specified') query = RasterRDR( aoi_group.lat_file, aoi_group.lon_file, - height_group.height_file_rdr, height_group.dem, + hgt_file=height_group.height_file_rdr, dem_file=height_group.dem, cube_spacing_in_m=cube_spacing_in_m ) diff --git a/tools/RAiDER/delay.py b/tools/RAiDER/delay.py index 16ef04d33..fc34564a5 100755 --- a/tools/RAiDER/delay.py +++ b/tools/RAiDER/delay.py @@ -39,7 +39,7 @@ def tropo_delay( los: LOS, height_levels: Optional[list[float]] = None, out_proj: Union[int, str] = 4326, - zref: Optional[np.float64] = None, + zref: Optional[float] = None, ): """Calculate integrated delays on query points. diff --git a/tools/RAiDER/dem.py b/tools/RAiDER/dem.py index effdc48f6..4d9fbfc04 100644 --- a/tools/RAiDER/dem.py +++ b/tools/RAiDER/dem.py @@ -19,7 +19,7 @@ def download_dem( ll_bounds: Union[tuple, List, np.ndarray]=None, - dem_path: Path=Path('warpedDEM.dem'), + dem_path: Path=Path('warpedDEM.rdr'), overwrite: bool=False, writeDEM: bool=False, buf: float=0.02, diff --git a/tools/RAiDER/llreader.py b/tools/RAiDER/llreader.py index e3068d260..15205f613 100644 --- a/tools/RAiDER/llreader.py +++ b/tools/RAiDER/llreader.py @@ -36,8 +36,8 @@ class AOI: _type - Type of AOI """ - def __init__(self, cube_spacing_in_m: Optional[float]=None) -> None: - self._output_directory = os.getcwd() + def __init__(self, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None: + self._output_directory = output_directory self._bounding_box = None self._proj = CRS.from_epsg(4326) self._geotransform = None @@ -194,8 +194,8 @@ def set_output_xygrid(self, dst_crs: Union[int, str]=4326) -> None: class StationFile(AOI): """Use a .csv file containing at least Lat, Lon, and optionally Hgt_m columns.""" - def __init__(self, station_file, demFile=None, cube_spacing_in_m: Optional[float]=None) -> None: - super().__init__(cube_spacing_in_m) + def __init__(self, station_file, demFile=None, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None: + super().__init__(cube_spacing_in_m, output_directory) self._filename = station_file self._demfile = demFile self._bounding_box = bounds_from_csv(station_file) @@ -237,15 +237,15 @@ def readZ(self): # write the elevations to the file df['Hgt_m'] = z_out df.to_csv(self._filename, index=False) - self.__init__(self._filename) + self._demfile = None return z_out class RasterRDR(AOI): """Use a 2-band raster file containing lat/lon coordinates.""" - def __init__(self, lat_file, lon_file=None, hgt_file=None, dem_file=None, convention='isce', cube_spacing_in_m: Optional[float]=None) -> None: - super().__init__(cube_spacing_in_m) + def __init__(self, lat_file, lon_file=None, *, hgt_file=None, dem_file=None, convention='isce', cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None: + super().__init__(cube_spacing_in_m, output_directory) self._type = 'radar_rasters' self._latfile = lat_file self._lonfile = lon_file @@ -310,8 +310,8 @@ def readZ(self) -> np.ndarray: class BoundingBox(AOI): """Parse a bounding box AOI.""" - def __init__(self, bbox, cube_spacing_in_m: Optional[float]=None) -> None: - super().__init__(cube_spacing_in_m) + def __init__(self, bbox, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None: + super().__init__(cube_spacing_in_m, output_directory) self._bounding_box = bbox self._type = 'bounding_box' @@ -323,8 +323,8 @@ class GeocodedFile(AOI): _bounding_box: BB.SNWE _is_dem: bool - def __init__(self, path: Path, is_dem=False, cube_spacing_in_m: Optional[float]=None) -> None: - super().__init__(cube_spacing_in_m) + def __init__(self, path: Path, is_dem=False, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None: + super().__init__(cube_spacing_in_m, output_directory) from RAiDER.utilFcns import rio_extents, rio_profile, rio_stats @@ -366,9 +366,9 @@ def readZ(self): class Geocube(AOI): """Pull lat/lon/height from a georeferenced data cube.""" - def __init__(self, path_cube, cube_spacing_in_m: Optional[float]=None) -> None: + def __init__(self, path_cube, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None: from RAiDER.utilFcns import rio_stats - super().__init__(cube_spacing_in_m) + super().__init__(cube_spacing_in_m, output_directory) self.path = path_cube self._type = 'Geocube' self._bounding_box = self.get_extent() diff --git a/tools/RAiDER/models/gmao.py b/tools/RAiDER/models/gmao.py old mode 100755 new mode 100644 index f1176ffe1..6f9a367a3 --- a/tools/RAiDER/models/gmao.py +++ b/tools/RAiDER/models/gmao.py @@ -1,19 +1,15 @@ import datetime as dt import shutil -import warnings from pathlib import Path -import h5py import numpy as np -import pydap.client +import xarray as xr from pyproj import CRS from RAiDER.logger import logger -from RAiDER.models.model_levels import ( - LEVELS_137_HEIGHTS, -) +from RAiDER.models.model_levels import LEVELS_137_HEIGHTS from RAiDER.models.weatherModel import TIME_RES, WeatherModel -from RAiDER.utilFcns import requests_retry_session, round_date, writeWeatherVarsXarray +from RAiDER.utilFcns import requests_retry_session, round_date, write_weather_vars_to_ds class GMAO(WeatherModel): @@ -80,80 +76,83 @@ def _fetch(self, out: Path) -> None: if corrected_DT >= T0: # open the dataset and pull the data url = 'https://opendap.nccs.nasa.gov/dods/GEOS-5/fp/0.25_deg/assim/inst3_3d_asm_Nv' - # For warning from pydap when using HTTPS instead of DAP2 or DAP4: - # pydap is incompatible with DAP data from opendap.nccs.nasa.gov. See issue #736 - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - ds = pydap.client.open_url(url) - - q = ( - ds['qv'] - .array[ + with xr.open_dataset(url, decode_times=False) as ds: + q = ds['qv'][ time_ind, ml_min : (ml_max + 1), lat_min_ind : (lat_max_ind + 1), lon_min_ind : (lon_max_ind + 1), ] - .data[0] - ) - p = ( - ds['pl'] - .array[ - time_ind, ml_min : (ml_max + 1), lat_min_ind : (lat_max_ind + 1), lon_min_ind : (lon_max_ind + 1) + p = ds['pl'][ + time_ind, + ml_min : (ml_max + 1), + lat_min_ind : (lat_max_ind + 1), + lon_min_ind : (lon_max_ind + 1), ] - .data[0] - ) - t = ( - ds['t'] - .array[ - time_ind, ml_min : (ml_max + 1), lat_min_ind : (lat_max_ind + 1), lon_min_ind : (lon_max_ind + 1) + t = ds['t'][ + time_ind, + ml_min : (ml_max + 1), + lat_min_ind : (lat_max_ind + 1), + lon_min_ind : (lon_max_ind + 1), ] - .data[0] - ) - h = ( - ds['h'] - .array[ - time_ind, ml_min : (ml_max + 1), lat_min_ind : (lat_max_ind + 1), lon_min_ind : (lon_max_ind + 1) + h = ds['h'][ + time_ind, + ml_min : (ml_max + 1), + lat_min_ind : (lat_max_ind + 1), + lon_min_ind : (lon_max_ind + 1), ] - .data[0] - ) else: root = 'https://portal.nccs.nasa.gov/datashare/gmao/geos-fp/das/Y{}/M{:02d}/D{:02d}' - base = f'GEOS.fp.asm.inst3_3d_asm_Nv.{corrected_DT.strftime("%Y%m%d")}_{corrected_DT.hour:02}00.V01.nc4' - URL = f'{root.format(corrected_DT.year, corrected_DT.month, corrected_DT.day)}/{base}' - path = Path(f'{out.stem}_raw{out.suffix}') - if not path.exists(): - logger.info('Fetching URL: %s', URL) - session = requests_retry_session() - resp = session.get(URL, stream=True) - assert resp.ok, f'Could not access url for datetime: {corrected_DT}' - with path.open('wb') as fout: - shutil.copyfileobj(resp.raw, fout) - else: - logger.warning('Weather model already exists, skipping download') - - with h5py.File(path, 'r') as ds: - q = ds['QV'][0, :, lat_min_ind : (lat_max_ind + 1), lon_min_ind : (lon_max_ind + 1)] - p = ds['PL'][0, :, lat_min_ind : (lat_max_ind + 1), lon_min_ind : (lon_max_ind + 1)] - t = ds['T'][0, :, lat_min_ind : (lat_max_ind + 1), lon_min_ind : (lon_max_ind + 1)] - h = ds['H'][0, :, lat_min_ind : (lat_max_ind + 1), lon_min_ind : (lon_max_ind + 1)] - path.unlink() - - lats = np.arange((-90 + lat_min_ind * self._lat_res), (-90 + (lat_max_ind + 1) * self._lat_res), self._lat_res) + filename = f'GEOS.fp.asm.inst3_3d_asm_Nv.{corrected_DT.strftime("%Y%m%d")}_{corrected_DT.hour:02}00.V01.nc4' + url = f'{root.format(corrected_DT.year, corrected_DT.month, corrected_DT.day)}/{filename}' + url += '#mode=bytes' # https://github.com/pydata/xarray/issues/3653#issuecomment-832712426 + with xr.open_dataset(url) as ds: + q = ds['QV'][ + 0, # time (always just 1) + ml_min : (ml_max + 1), # lev + lat_min_ind : (lat_max_ind + 1), # lat + lon_min_ind : (lon_max_ind + 1), # lon + ] + p = ds['PL'][ + 0, + ml_min : (ml_max + 1), + lat_min_ind : (lat_max_ind + 1), + lon_min_ind : (lon_max_ind + 1), + ] + t = ds['T'][ + 0, + ml_min : (ml_max + 1), + lat_min_ind : (lat_max_ind + 1), + lon_min_ind : (lon_max_ind + 1), + ] + h = ds['H'][ + 0, + ml_min : (ml_max + 1), + lat_min_ind : (lat_max_ind + 1), + lon_min_ind : (lon_max_ind + 1), + ] + + lats = np.arange( + -90 + lat_min_ind * self._lat_res, + -90 + (lat_max_ind + 1) * self._lat_res, + self._lat_res, + ) lons = np.arange( - (-180 + lon_min_ind * self._lon_res), (-180 + (lon_max_ind + 1) * self._lon_res), self._lon_res + -180 + lon_min_ind * self._lon_res, + -180 + (lon_max_ind + 1) * self._lon_res, + self._lon_res, ) lon, lat = np.meshgrid(lons, lats) try: # Note that lat/lon gets written twice for GMAO because they are the same as y/x - writeWeatherVarsXarray(lat, lon, h, q, p, t, self._time, self._proj, out) + write_weather_vars_to_ds(lat, lon, h, q, p, t, self._time, self._proj, out) except: logger.exception('Unable to save weathermodel to file:') raise - def load_weather(self, f=None) -> None: + def load_weather(self, f=None, *args, **kwargs) -> None: """ Consistent class method to be implemented across all weather model types. As a result of calling this method, all of the variables (x, y, z, p, q, @@ -163,21 +162,18 @@ def load_weather(self, f=None) -> None: f = self.files[0] if f is None else f self._load_model_level(f) - 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 - 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'][:]) + def _load_model_level(self, filename: Path) -> None: + """Get the variables from the GMAO link using OPeNDAP.""" + with xr.open_dataset(filename) as ds: + lons = ds['x'] + lats = ds['y'] + h = ds['h'].data + q = ds['q'].data + p = ds['p'].data + t = ds['t'].data # restructure the 1-D lat/lon in regular 2D grid - _lons, _lats = np.meshgrid(lons, lats) + lons, lats = np.meshgrid(lons, lats) # Re-structure everything from (heights, lats, lons) to (lons, lats, heights) p = np.transpose(p) @@ -199,12 +195,11 @@ def _load_model_level(self, filename) -> None: h = np.flip(h, 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 diff --git a/tools/RAiDER/models/hrrr.py b/tools/RAiDER/models/hrrr.py index 40f2b5ff5..50619a616 100644 --- a/tools/RAiDER/models/hrrr.py +++ b/tools/RAiDER/models/hrrr.py @@ -1,14 +1,15 @@ import datetime as dt -import os from pathlib import Path +from typing import List, Tuple, Union, cast +from RAiDER.types import BB, Array2D, FloatArray2D import geopandas as gpd +import herbie +import herbie.accessors import numpy as np import xarray as xr -from herbie import Herbie from pyproj import CRS, Transformer from shapely.geometry import Polygon, box -from typing import Optional, Union, List, Tuple from RAiDER.logger import logger from RAiDER.models.customExceptions import NoWeatherModelData @@ -29,13 +30,13 @@ def check_hrrr_dataset_availability(datetime: dt.datetime, model='hrrr') -> bool: """Note a file could still be missing within the models valid range.""" - herbie = Herbie( + h = herbie.Herbie( datetime, model=model, product='nat', fxx=0, ) - return herbie.grib_source is not None + return h.grib_source is not None def download_hrrr_file(ll_bounds, DATE, out: Path, model='hrrr', product='nat', fxx=0, verbose=False) -> None: @@ -54,7 +55,7 @@ def download_hrrr_file(ll_bounds, DATE, out: Path, model='hrrr', product='nat', Returns: None, writes data to a netcdf file """ - herbie = Herbie( + h = herbie.Herbie( DATE.strftime('%Y-%m-%d %H:%M'), model=model, product=product, @@ -66,7 +67,12 @@ def download_hrrr_file(ll_bounds, DATE, out: Path, model='hrrr', product='nat', # Iterate through the list of datasets try: - ds_list = herbie.xarray(':(SPFH|PRES|TMP|HGT):', verbose=verbose) + # cast: Herbie.xarray can return one Dataset or a list + ds_list = cast( + Union[xr.Dataset, List[xr.Dataset]], + h.xarray(':(SPFH|PRES|TMP|HGT):', verbose=verbose), + ) + assert isinstance(ds_list, list) except ValueError as e: logger.error(e) raise @@ -74,14 +80,14 @@ def download_hrrr_file(ll_bounds, DATE, out: Path, model='hrrr', product='nat', # Note order coord names are request for `test_HRRR_ztd` matters # when both coord names are retreived by Herbie is ds_list possibly in # Different orders on different machines; `hybrid` is what is expected for the test. - ds_list_filt_0 = [ds for ds in ds_list if 'hybrid' in ds._coord_names] - ds_list_filt_1 = [ds for ds in ds_list if 'isobaricInhPa' in ds._coord_names] - if ds_list_filt_0: + ds_list_filt_0 = [ds for ds in ds_list if 'hybrid' in ds.coords] + ds_list_filt_1 = [ds for ds in ds_list if 'isobaricInhPa' in ds.coords] + if len(ds_list_filt_0) > 0: ds_out = ds_list_filt_0[0] coord = 'hybrid' # I do not think that this coord name will result in successful processing nominally as variables are # gh,gribfile_projection for test_HRRR_ztd - elif ds_list_filt_1: + elif len(ds_list_filt_1) > 0: ds_out = ds_list_filt_1[0] coord = 'isobaricInhPa' else: @@ -91,22 +97,25 @@ def download_hrrr_file(ll_bounds, DATE, out: Path, model='hrrr', product='nat', try: x_min, x_max, y_min, y_max = get_bounds_indices( ll_bounds, - ds_out.latitude.to_numpy(), - ds_out.longitude.to_numpy(), + ds_out['latitude'].to_numpy(), + ds_out['longitude'].to_numpy(), ) except NoWeatherModelData as e: logger.error(e) - logger.error('lat/lon bounds: %s', ll_bounds) - logger.error('Weather model is {}'.format(model)) + logger.error(f'lat/lon bounds: {ll_bounds}') + logger.error(f'Weather model is {model}') raise # bookkeepping ds_out = ds_out.rename({'gh': 'z', coord: 'levels'}) + # Herbie injects brainworms into all its xarray Datasets + crs = cast(herbie.accessors.HerbieAccessor, ds_out.herbie).crs + # projection information ds_out['proj'] = 0 - for k, v in CRS.from_user_input(ds_out.herbie.crs).to_cf().items(): - ds_out.proj.attrs[k] = v + for k, v in CRS.from_user_input(crs).to_cf().items(): + ds_out['proj'].attrs[k] = v for var in ds_out.data_vars: ds_out[var].attrs['grid_mapping'] = 'proj' @@ -121,53 +130,50 @@ def download_hrrr_file(ll_bounds, DATE, out: Path, model='hrrr', product='nat', t = Transformer.from_crs(4326, proj, always_xy=True) xl, yl = t.transform(ds_out['longitude'].values, ds_out['latitude'].values) - W, E, S, N = np.nanmin(xl), np.nanmax(xl), np.nanmin(yl), np.nanmax(yl) + s, n, w, e = np.nanmin(yl), np.nanmax(yl), np.nanmin(xl), np.nanmax(xl) grid_x = 3000 # meters grid_y = 3000 # meters - xs = np.arange(W, E + grid_x / 2, grid_x) - ys = np.arange(S, N + grid_y / 2, grid_y) + xs = np.arange(w, e + grid_x / 2, grid_x) + ys = np.arange(s, n + grid_y / 2, grid_y) ds_out['x'] = xs ds_out['y'] = ys ds_sub = ds_out.isel(x=slice(x_min, x_max), y=slice(y_min, y_max)) - ds_sub.to_netcdf(out, engine='netcdf4') + ds_sub.to_netcdf(out) -def get_bounds_indices(SNWE, lats, lons): +def get_bounds_indices(snwe: BB.SNWE, lats: FloatArray2D, lons: FloatArray2D) -> BB.WSEN: """Convert SNWE lat/lon bounds to index bounds.""" # Unpack the bounds and find the relevent indices - S, N, W, E = SNWE + s, n, w, e = snwe # Need to account for crossing the international date line - if W < E: - m1 = (S <= lats) & (N >= lats) & (W <= lons) & (E >= lons) - else: + if w >= e: raise ValueError( - 'Longitude is either flipped or you are crossing the international date line;' - + 'if the latter please give me longitudes from 0-360' + 'Longitude is either flipped or you are crossing the international date line; ' + 'if the latter please give me longitudes from 0-360' ) - if np.sum(m1) == 0: + m1: Array2D[np.bool] = (s <= lats) & (n >= lats) & (w <= lons) & (e >= lons) + if np.sum(m1) == 0: # All False lons = np.mod(lons, 360) - W, E = np.mod([W, E], 360) - m1 = (S <= lats) & (N >= lats) & (W <= lons) & (E >= lons) - if np.sum(m1) == 0: + w, e = np.mod([w, e], 360) + # try again + m1 = (s <= lats) & (n >= lats) & (w <= lons) & (e >= lons) + if np.sum(m1) == 0: # All False raise NoWeatherModelData('Area of Interest has no overlap with the HRRR model available extent') # Y extent - shp = lats.shape - m1_y = np.argwhere(np.sum(m1, axis=1) != 0) - y_min = max(m1_y[0][0], 0) - y_max = min(m1_y[-1][0], shp[0]) - m1_y = None + m1_y: Array2D[np.integer] = np.argwhere(np.sum(m1, axis=1) != 0) + y_min: int = max(m1_y[0][0], 0) + y_max: int = min(m1_y[-1][0], lats.shape[0]) + del m1_y # big # X extent - m1_x = np.argwhere(np.sum(m1, axis=0) != 0) - x_min = max(m1_x[0][0], 0) - x_max = min(m1_x[-1][0], shp[1]) - m1_x = None - m1 = None + m1_x: Array2D[np.integer] = np.argwhere(np.sum(m1, axis=0) != 0) + x_min: int = max(m1_x[0][0], 0) + x_max: int = min(m1_x[-1][0], lats.shape[1]) return x_min, x_max, y_min, y_max diff --git a/tools/RAiDER/models/merra2.py b/tools/RAiDER/models/merra2.py index c35f0b4e9..c3c56990c 100755 --- a/tools/RAiDER/models/merra2.py +++ b/tools/RAiDER/models/merra2.py @@ -3,7 +3,6 @@ from pathlib import Path import numpy as np -import pydap.client import xarray as xr from pyproj import CRS @@ -12,7 +11,7 @@ LEVELS_137_HEIGHTS, ) from RAiDER.models.weatherModel import WeatherModel -from RAiDER.utilFcns import writeWeatherVarsXarray +from RAiDER.utilFcns import write_weather_vars_to_ds # Path to Netrc file, can be controlled by env var @@ -69,7 +68,11 @@ def __init__(self) -> None: self._proj = CRS.from_epsg(4326) 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.""" + """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 # check whether the file already exists @@ -100,28 +103,45 @@ def _fetch(self, out: Path) -> None: # open the dataset and pull the data url = ( - 'dap4://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' + f'dap4://goldsmr5.gesdisc.eosdis.nasa.gov/opendap/MERRA2/M2T3NVASM.5.12.4/' + f'{time.strftime('%Y/%m')}/' + f'MERRA2_{url_sub}.tavg3_3d_asm_Nv.{time.strftime('%Y%m%d')}.nc4' ) - stream = pydap.client.open_url(url) - - 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() + # pydap engine required for password-protected datasets + ds = xr.open_dataset(url, decode_times=False, engine='pydap') + + q = ds['QV'][ + 0, + :, + lat_min_ind : lat_max_ind + 1, + lon_min_ind : lon_max_ind + 1, + ].data.squeeze() + p = ds['PL'][ + 0, + :, + lat_min_ind : lat_max_ind + 1, + lon_min_ind : lon_max_ind + 1, + ].data.squeeze() + t = ds['T'][ + 0, + :, + lat_min_ind : lat_max_ind + 1, + lon_min_ind : lon_max_ind + 1, + ].data.squeeze() + h = ds['H'][ + 0, + :, + lat_min_ind : lat_max_ind + 1, + lon_min_ind : lon_max_ind + 1, + ].data.squeeze() try: - writeWeatherVarsXarray(lat, lon, h, q, p, t, time, self._proj, out_path=out) + write_weather_vars_to_ds(lat, lon, h, q, p, t, time, self._proj, out_path=out) except Exception as e: logger.debug(e) - logger.exception('MERRA-2: Unable to save weathermodel to file') - raise RuntimeError(f'MERRA-2 failed with the following error: {e}') + logger.exception('MERRA-2: Unable to save weather model query to file') + raise def load_weather(self, f=None, *args, **kwargs) -> None: """ @@ -134,7 +154,7 @@ def load_weather(self, f=None, *args, **kwargs) -> None: self._load_model_level(f) def _load_model_level(self, filename) -> None: - """Get the variables from the GMAO link using OpenDAP.""" + """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 diff --git a/tools/RAiDER/models/ncmr.py b/tools/RAiDER/models/ncmr.py index 58338f398..3918017da 100755 --- a/tools/RAiDER/models/ncmr.py +++ b/tools/RAiDER/models/ncmr.py @@ -18,7 +18,7 @@ from RAiDER.utilFcns import ( read_NCMR_loginInfo, show_progress, - writeWeatherVarsXarray, + write_weather_vars_to_ds, ) @@ -193,7 +193,7 @@ def _download_ncmr_file(self, out: Path, date_time, bounding_box) -> None: ######################################################################################################################## try: - writeWeatherVarsXarray(lats, lons, hgt, q, p, t, self._time, self._proj, out_path=out) + write_weather_vars_to_ds(lats, lons, hgt, q, p, t, self._time, self._proj, out_path=out) except: logger.exception('Unable to save weathermodel to file') diff --git a/tools/RAiDER/models/template.py b/tools/RAiDER/models/template.py index e7b2cf7b0..8f0df1c7f 100644 --- a/tools/RAiDER/models/template.py +++ b/tools/RAiDER/models/template.py @@ -77,7 +77,7 @@ def _fetch(self, out) -> None: # download dataset of the custom weather model "ABCD" from a server and then save it to a file named out. # This function needs to be writen by the users. For download from the weather model server, the weather model # name, time and bounding box may be needed to retrieve the dataset; for cases where no file is actually - # downloaded, e.g. the GMAO and MERRA-2 models using OpenDAP, this function can be omitted leaving the data + # downloaded, e.g. the GMAO and MERRA-2 models using OPeNDAP, this function can be omitted leaving the data # retrieval to the following "load_weather" function. self._files = self._download_abcd_file(out, 'abcd', self._time, self._ll_bounds) @@ -91,7 +91,7 @@ def load_weather(self, filename) -> None: # read individual variables (in 3-D cube format with exactly the same dimension) from downloaded file # This function needs to be writen by the users. For downloaded file from the weather model server, # this function extracts the individual variables from the saved file named filename; - # for cases where no file is actually downloaded, e.g. the GMAO and MERRA-2 models using OpenDAP, + # for cases where no file is actually downloaded, e.g. the GMAO and MERRA-2 models using OPeNDAP, # this function retrieves the individual variables directly from the weblink of the weather model. lats, lons, xs, ys, t, q, p, hgt = self._makeDataCubes(filename) diff --git a/tools/RAiDER/models/weatherModel.py b/tools/RAiDER/models/weatherModel.py index 240e184dd..b8d7c7103 100755 --- a/tools/RAiDER/models/weatherModel.py +++ b/tools/RAiDER/models/weatherModel.py @@ -472,22 +472,20 @@ def checkValidBounds( def checkContainment(self, ll_bounds: Union[List, Tuple,np.ndarray], buffer_deg: float = 1e-5) -> bool: """ - Checks containment of weather model bbox. + Checks containment of weather model bbox. Args: - ---------- - weather_model : WeatherModel - ll_bounds: an array of floats (SNWE) demarcating bbox of targets - buffer_deg : float - For x-translates for extents that lie outside of world bounding box, - this ensures that translates have some overlap. The default is 1e-5 - or ~11.1 meters. + weather_model: WeatherModel + ll_bounds: an array of floats (SNWE) demarcating bbox of targets + buffer_deg: + For x-translates for extents that lie outside of world bounding box, + this ensures that translates have some overlap. The default is 1e-5 + or ~11.1 meters. Returns: - ------- - bool - True if weather model contains bounding box of OutLats and outLons - and False otherwise. + bool: + True if weather model contains bounding box of OutLats and + outLons and False otherwise. """ # Parse the input ymin_input, ymax_input, xmin_input, xmax_input = ll_bounds @@ -495,7 +493,7 @@ def checkContainment(self, ll_bounds: Union[List, Tuple,np.ndarray], buffer_deg: # Parse the weather model bounding box input_box = box(xmin_input, ymin_input, xmax_input, ymax_input) - xmin, ymin, xmax, ymax = self.bbox + xmin, ymin, xmax, ymax = self._ll_bounds weather_model_box = box(xmin, ymin, xmax, ymax) # Logger @@ -522,7 +520,7 @@ def checkContainment(self, ll_bounds: Union[List, Tuple,np.ndarray], buffer_deg: if weather_model_box.contains(world_box): # Handle the case where the whole world is requested - self.bbox = (-180, -90, 180, 90) + self._ll_bounds = (-180, -90, 180, 90) return True else: if weather_model_box.contains(input_box): diff --git a/tools/RAiDER/processWM.py b/tools/RAiDER/processWM.py index d8ab4b80d..0debdbd70 100755 --- a/tools/RAiDER/processWM.py +++ b/tools/RAiDER/processWM.py @@ -18,11 +18,16 @@ ExistingWeatherModelTooSmall, TryToKeepGoingError, ) -from RAiDER.models.weatherModel import checkContainment_raw, make_raw_weather_data_filename, make_weather_model_filename +from RAiDER.models.weatherModel import ( + WeatherModel, + checkContainment_raw, + make_raw_weather_data_filename, + make_weather_model_filename, +) def prepareWeatherModel( - weather_model, + weather_model: WeatherModel, time, ll_bounds, download_only: bool = False, diff --git a/tools/RAiDER/types/__init__.py b/tools/RAiDER/types/__init__.py index 58d00736b..036cbdd98 100644 --- a/tools/RAiDER/types/__init__.py +++ b/tools/RAiDER/types/__init__.py @@ -1,10 +1,32 @@ -"""Types specific to RAiDER.""" +"""Helper types used throughout RAiDER.""" -from typing import Literal, Union +from typing import Literal, TypeVar, Union +import numpy as np from pyproj import CRS LookDir = Literal['right', 'left'] TimeInterpolationMethod = Literal['none', 'center_time', 'azimuth_time_grid'] CRSLike = Union[CRS, str, int] + + +# From numpy +_ScalarT_co = TypeVar("_ScalarT_co", bound=np.generic, covariant=True) + +# Helpers for type annotating the dimensions of a numpy array. +# When you use these, your type checker will alert you when an annotated array +# is used in a way inconsistent with its dimensions. +Array1D = np.ndarray[tuple[int], np.dtype[_ScalarT_co]] +Array2D = np.ndarray[tuple[int, int], np.dtype[_ScalarT_co]] +Array3D = np.ndarray[tuple[int, int, int], np.dtype[_ScalarT_co]] +# ... (repeat the pattern as needed for higher dimensions) + +# Any number of dimensions -- when ndim is not able to be known ahead of time +ArrayND = np.ndarray[tuple[int, ...], np.dtype[_ScalarT_co]] + + +FloatArray1D = Array1D[np.floating] +FloatArray2D = Array2D[np.floating] +FloatArray3D = Array3D[np.floating] +FloatArrayND = ArrayND[np.floating] diff --git a/tools/RAiDER/utilFcns.py b/tools/RAiDER/utilFcns.py index 0dee94634..8a538f04d 100644 --- a/tools/RAiDER/utilFcns.py +++ b/tools/RAiDER/utilFcns.py @@ -22,7 +22,7 @@ from RAiDER.constants import _g1 as G1 from RAiDER.llreader import AOI from RAiDER.logger import logger -from RAiDER.types import BB, RIO, CRSLike +from RAiDER.types import BB, RIO, CRSLike, FloatArray2D, FloatArray3D # Optional imports @@ -638,68 +638,76 @@ def requests_retry_session(retries: int=10, session=None): # noqa: ANN001, ANN2 return session -def writeWeatherVarsXarray( - lat: npt.NDArray[np.floating], - lon: npt.NDArray[np.floating], - h: npt.NDArray[np.floating], - q: npt.NDArray[np.floating], - p: npt.NDArray[np.floating], - t: npt.NDArray[np.floating], +def write_weather_vars_to_ds( + lat_range: FloatArray2D, + lon_range: FloatArray2D, + h: FloatArray3D, + q: FloatArray3D, + p: FloatArray3D, + t: FloatArray3D, datetime: dt.datetime, crs: CRS, out_path: Path, NoDataValue: int = -9999, - chunksize: Tuple[int, ...] = (1, 91, 144), + chunksize: Tuple[int, int, int] = (1, 91, 144), ) -> None: - assert len(h.shape) == 3, "Invalid h array dimensions" - assert len(q.shape) == 3, "Invalid q array dimensions" - assert len(p.shape) == 3, "Invalid p array dimensions" - assert len(t.shape) == 3, "Invalid t array dimensions" - dataset_dict = { - 'h': (('z', 'y', 'x'), h), - 'q': (('z', 'y', 'x'), q), - 'p': (('z', 'y', 'x'), p), - 't': (('z', 'y', 'x'), t), - } - assert len(lat.shape) == 2, "Invalid lats array dimensions" - assert len(lon.shape) == 2, "Invalid lons array dimensions" - dimension_dict = { - 'latitude': (('y', 'x'), lat), - 'longitude': (('y', 'x'), lon), - } - # I added datetime as an input to the function and just copied these two lines from merra2 for the attrs_dict - attrs_dict = { - 'datetime': datetime.strftime('%Y_%m_%dT%H_%M_%S'), - 'date_created': datetime.now().strftime('%Y_%m_%dT%H_%M_%S'), - 'NoDataValue': NoDataValue, - 'chunksize': chunksize, - # 'mapping_name': mapping_name, - } - - ds = xr.Dataset( - data_vars=dataset_dict, - coords=dimension_dict, - attrs=attrs_dict, - ) - - ds['h'].attrs['standard_name'] = 'mid_layer_heights' - ds['p'].attrs['standard_name'] = 'mid_level_pressure' - ds['q'].attrs['standard_name'] = 'specific_humidity' - ds['t'].attrs['standard_name'] = 'air_temperature' - - ds['h'].attrs['units'] = 'm' - ds['p'].attrs['units'] = 'Pa' - ds['q'].attrs['units'] = 'kg kg-1' - ds['t'].attrs['units'] = 'K' - - ds['proj'] = 0 - for k, v in crs.to_cf().items(): - ds.proj.attrs[k] = v - for var in ds.data_vars: - ds[var].attrs['grid_mapping'] = 'proj' - - ds.to_netcdf(out_path) - del ds + DIMS_3D = ('z', 'y', 'x') + DIMS_2D = ('y', 'x') + with xr.Dataset( + data_vars={ + 'h': xr.Variable( + dims=DIMS_3D, + data=h, + attrs={ + 'standard_name': 'mid_layer_heights', + 'units': 'm', + }, + ), + 'q': xr.Variable( + dims=DIMS_3D, + data=q, + attrs={ + 'standard_name': 'specific_humidity', + 'units': 'kg kg-1', + }, + ), + 'p': xr.Variable( + dims=DIMS_3D, + data=p, + attrs={ + 'standard_name': 'mid_level_pressure', + 'units': 'Pa', + }, + ), + 't': xr.Variable( + dims=DIMS_3D, + data=t, + attrs={ + 'standard_name': 'air_temperature', + 'units': 'K', + }, + ), + 'proj': 0, + }, + coords={ + 'latitude': xr.Variable(dims=DIMS_2D, data=lat_range), + 'longitude': xr.Variable(dims=DIMS_2D, data=lon_range), + }, + attrs={ + # I added datetime as an input to the function and just copied these + # two lines from merra2 for the attrs_dict + 'datetime': datetime.strftime('%Y_%m_%dT%H_%M_%S'), + 'date_created': dt.datetime.now().strftime('%Y_%m_%dT%H_%M_%S'), + 'NoDataValue': NoDataValue, + 'chunksize': chunksize, + # 'mapping_name': mapping_name, + }, + ) as ds: + for k, v in crs.to_cf().items(): + ds['proj'].attrs[k] = v + for var in ds.data_vars: + ds[var].attrs['grid_mapping'] = 'proj' + ds.to_netcdf(out_path) def convertLons(inLons: np.ndarray) -> np.ndarray: