From 02c720a911f01879669d9c0404568169965e5039 Mon Sep 17 00:00:00 2001 From: bbuzz31 Date: Mon, 23 Jan 2023 18:32:57 -0800 Subject: [PATCH 01/19] Setup structure for interpolating wm times --- pyproject.toml | 1 + tools/RAiDER/checkArgs.py | 1 + tools/RAiDER/cli/__main__.py | 4 +- tools/RAiDER/cli/raider.py | 173 +++++++++++++++++++++++++++++++++ tools/RAiDER/cli/validators.py | 2 + tools/RAiDER/constants.py | 1 + tools/RAiDER/delay.py | 17 +++- tools/RAiDER/utilFcns.py | 40 +++++++- 8 files changed, 233 insertions(+), 6 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 3a69e4d44..57a61db1b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -43,6 +43,7 @@ repository = "https://github.com/dbekaert/RAiDER" [project.scripts] "raider.py" = "RAiDER.cli.__main__:main" "calcDelays.py" = "RAiDER.cli.__main__:calcDelays" +"calcDelaysInterp.py" = "RAiDER.cli.__main__:calcDelaysInterp" "calcDelaysGUNW.py" = "RAiDER.cli.__main__:calcDelaysGUNW" "raiderDownloadGNSS.py" = "RAiDER.cli.__main__:downloadGNSS" "downloadGNSS.py" = "RAiDER.cli.__main__:downloadGNSS" diff --git a/tools/RAiDER/checkArgs.py b/tools/RAiDER/checkArgs.py index 661ca4151..46464164c 100644 --- a/tools/RAiDER/checkArgs.py +++ b/tools/RAiDER/checkArgs.py @@ -18,6 +18,7 @@ from RAiDER.logger import logger +## put it in here def checkArgs(args): ''' Helper fcn for checking argument compatibility and returns the diff --git a/tools/RAiDER/cli/__main__.py b/tools/RAiDER/cli/__main__.py index 3b195d957..fd863cd3b 100644 --- a/tools/RAiDER/cli/__main__.py +++ b/tools/RAiDER/cli/__main__.py @@ -2,7 +2,7 @@ import sys from importlib.metadata import entry_points -from RAiDER.cli.raider import calcDelays, downloadGNSS, calcDelaysGUNW +from RAiDER.cli.raider import calcDelays, calcDelaysInterp, downloadGNSS, calcDelaysGUNW def main(): @@ -11,7 +11,7 @@ def main(): formatter_class=argparse.ArgumentDefaultsHelpFormatter ) parser.add_argument( - '++process', choices=['calcDelays', 'downloadGNSS', 'calcDelaysGUNW'], + '++process', choices=['calcDelays', 'calcDelaysInterp' 'downloadGNSS', 'calcDelaysGUNW'], default='calcDelays', help='Select the entrypoint to use' ) diff --git a/tools/RAiDER/cli/raider.py b/tools/RAiDER/cli/raider.py index 3f471a867..f33027ef4 100644 --- a/tools/RAiDER/cli/raider.py +++ b/tools/RAiDER/cli/raider.py @@ -218,6 +218,8 @@ def calcDelays(iargs=None): logger.debug('Starting to run the weather model calculation') logger.debug(f'Date: {t.strftime("%Y%m%d")}') logger.debug('Beginning weather model pre-processing') + + try: weather_model_file = prepareWeatherModel( model, t, @@ -287,6 +289,177 @@ def calcDelays(iargs=None): return wet_filenames +def calcDelaysInterp(iargs=None): + """ + Same as calcDelays, except will calculate at nearest weather model hour then interpolate + """ + import RAiDER + from RAiDER.delay import tropo_delay, avg_delays + from RAiDER.checkArgs import checkArgs + from RAiDER.processWM import prepareWeatherModel + from RAiDER.utilFcns import writeDelays, get_nearest_wmtimes + examples = 'Examples of use:' \ + '\n\t raider.py ++process calcDelaysInterp customTemplatefile.cfg' \ + '\n\t calcDelaysInterp.py customTemplateFile' + + p = argparse.ArgumentParser( + description = + 'Command line interface for RAiDER processing with a configure file.' + 'Default options can be found by running: raider.py --generate_config', + epilog=examples, formatter_class=argparse.RawDescriptionHelpFormatter) + + p.add_argument( + 'customTemplateFile', nargs='?', + help='custom template with option settings.\n' + + "ignored if the default smallbaselineApp.cfg is input." + ) + + p.add_argument( + '--download_only', action='store_true', + help='only download a weather model.' + ) + + ## if not None, will replace first argument (customTemplateFile) + args = p.parse_args(args=iargs) + + # default input file + template_file = os.path.join(os.path.dirname(RAiDER.__file__), + 'cli', 'raider.yaml') + + # check: existence of input template files + if (not args.customTemplateFile + and not os.path.isfile(os.path.basename(template_file)) + and not args.generate_template): + msg = "No template file found! It requires that either:" + msg += "\n a custom template file, OR the default template " + msg += "\n file 'raider.yaml' exists in current directory." + + p.print_usage() + print(examples) + raise SystemExit(f'ERROR: {msg}') + + if args.customTemplateFile: + # check the existence + if not os.path.isfile(args.customTemplateFile): + raise FileNotFoundError(args.customTemplateFile) + + args.customTemplateFile = os.path.abspath(args.customTemplateFile) + else: + args.customTemplateFile = template_file + + # Read the template file + params = read_template_file(args.customTemplateFile) + + # Argument checking + params = checkArgs(params) + dl_only = True if params['download_only'] or args.download_only else False + + if not params.verbose: + logger.setLevel(logging.INFO) + + wet_filenames = [] + for t, w, f in zip( + params['date_list'], + params['wetFilenames'], + params['hydroFilenames'] + ): + + los = params['los'] + aoi = params['aoi'] + model = params['weather_model'] + + # add a buffer for raytracing + if los.ray_trace(): + ll_bounds = aoi.add_buffer(buffer=1) + else: + ll_bounds = aoi.add_buffer(buffer=1) + + ########################################################### + # weather model calculation + logger.debug('Starting to run the weather model calculation') + logger.debug(f'Date: {t.strftime("%Y%m%d")}') + logger.debug('Beginning weather model pre-processing') + + + + time1, time2 = get_nearest_wmtimes(t, model._Name) + wetDelays, hydroDelays = [], [] + for tt in [time1, time2]: + try: + weather_model_file = prepareWeatherModel( + model, tt, + ll_bounds=ll_bounds, # SNWE + wmLoc=params['weather_model_directory'], + makePlots=params['verbose'], + ) + except RuntimeError: + logger.exception("Date %s failed", tt) + continue + + + # dont process the delays for download only + if dl_only: + continue + + # Now process the delays + try: + wet_delay, hydro_delay = tropo_delay( + t, weather_model_file, aoi, los, + height_levels = params['height_levels'], + out_proj = params['output_projection'], + look_dir = params['look_dir'], + cube_spacing_m = params['cube_spacing_in_m'], + ) + except RuntimeError: + logger.exception("Date %s failed", t) + continue + + ########################################################### + # store in memory + wetDelays.append(wet_delay) + hydroDelays.append(hydro_delay) + + wet_delay, hydro_delay = avg_delays(wetDelays, hydroDelays, t, time1, time2) + + + # Write the delays to file + # Different options depending on the inputs + if los.is_Projected(): + out_filename = w.replace("_ztd", "_std") + f = f.replace("_ztd", "_std") + elif los.ray_trace(): + out_filename = w.replace("_std", "_ray") + f = f.replace("_std", "_ray") + else: + out_filename = w + + if hydro_delay is None: + # means that a dataset was returned + ds = wet_delay + ext = os.path.splitext(out_filename)[1] + out_filename = out_filename.replace('wet', 'tropo') + + if ext not in ['.nc', '.h5']: + out_filename = f'{os.path.splitext(out_filename)[0]}.nc' + + + if out_filename.endswith(".nc"): + ds.to_netcdf(out_filename, mode="w") + elif out_filename.endswith(".h5"): + ds.to_netcdf(out_filename, engine="h5netcdf", invalid_netcdf=True) + logger.info('Wrote delays to: %s', out_filename) + + else: + if aoi.type() == 'station_file': + out_filename = f'{os.path.splitext(out_filename)[0]}.csv' + + if aoi.type() in ['station_file', 'radar_rasters', 'geocoded_file']: + writeDelays(aoi, wet_delay, hydro_delay, out_filename, f, outformat=params['raster_format']) + + wet_filenames.append(out_filename) + + return wet_filenames + ## ------------------------------------------------------ downloadGNSSDelays.py def downloadGNSS(): """Parse command line arguments using argparse.""" diff --git a/tools/RAiDER/cli/validators.py b/tools/RAiDER/cli/validators.py index e07a778b1..6d6dc0b93 100755 --- a/tools/RAiDER/cli/validators.py +++ b/tools/RAiDER/cli/validators.py @@ -202,6 +202,8 @@ def parse_dates(arg_dict): l = arg_dict['date_list'] if isinstance(l, str): l = re.findall('[0-9]+', l) + elif isinstance(l, int): + l = [l] L = [enforce_valid_dates(d) for d in l] else: diff --git a/tools/RAiDER/constants.py b/tools/RAiDER/constants.py index 7d7c5cbe7..8b9c84913 100644 --- a/tools/RAiDER/constants.py +++ b/tools/RAiDER/constants.py @@ -18,3 +18,4 @@ R_EARTH_MAX = 6378137 R_EARTH_MIN = 6356752 + diff --git a/tools/RAiDER/delay.py b/tools/RAiDER/delay.py index da986ce92..d4f5837c4 100755 --- a/tools/RAiDER/delay.py +++ b/tools/RAiDER/delay.py @@ -43,7 +43,7 @@ def tropo_delay( look_dir: str='right', ): """ - Calculate integrated delays on query points. Options are: + Calculate integrated delays on query points. Options are: 1. Zenith delays (ZTD) 2. Zenith delays projected to the line-of-sight (STD-projected) 3. Slant delays integrated along the raypath (STD-raytracing) @@ -70,7 +70,7 @@ def tropo_delay( wm_proj = CRS.from_epsg(4326) else: wm_proj = CRS.from_wkt(wm_proj.to_wkt()) - + # get heights if height_levels is None: if aoi.type() == 'Geocube': @@ -128,7 +128,7 @@ def _get_delays_on_cube(dt, weather_model_file, wm_proj, ll_bounds, heights, los snwe = transform_bbox(ll_bounds, src_crs=4326, dest_crs=crs) out_spacing = get_output_spacing(cube_spacing_m, weather_model_file, wm_proj, crs) out_snwe = clip_bbox(snwe, out_spacing) - + logger.debug(f"Output SNWE: {out_snwe}") logger.debug(f"Output cube spacing: {out_spacing}") @@ -412,3 +412,14 @@ def _build_cube_ray( if output_created_here: return outputArrs + + +def avg_delays(wetDelays, hydroDelays, t0, t1, t2): + if hydroDelays[0] is None: + # do something with wet + pass + else: + # do something with both + pass + + return wetDelays[0], hydroDelays[0] \ No newline at end of file diff --git a/tools/RAiDER/utilFcns.py b/tools/RAiDER/utilFcns.py index 987a9cbfc..533d848e8 100755 --- a/tools/RAiDER/utilFcns.py +++ b/tools/RAiDER/utilFcns.py @@ -291,7 +291,7 @@ def writeResultsToXarray(dt, xpts, ypts, zpts, crs, wetDelay, hydroDelay, weathe ds.x.attrs["standard_name"] = "projection_x_coordinate" ds.x.attrs["long_name"] = "x-coordinate in projected coordinate system" ds.x.attrs["units"] = "m" - + return ds @@ -1135,3 +1135,41 @@ def transform_coords(proj1, proj2, x, y): """ transformer = Transformer.from_crs(proj1, proj2, always_xy=True) return transformer.transform(x, y) + + +def get_nearest_wmtimes(dt0, wm): + """" Get the nearest model times""" + import pandas as pd + from dateutil.relativedelta import relativedelta + from datetime import time as dtime + ## hourly time step availability - not sure about ECMWF / HRRR + wm_hour = {'GMAO': 3, 'ERA5':1, 'ERA5T': 1, 'HRES':1, 'HRRR': 1, 'NCMR': '?'} + wm_hours = np.arange(0, 24, wm_hour[wm]) + + dates = [] + times = [] + + ref_time = dt0 + ## sort the model and GUNW times + # date = dt0.date() + mod_times = [datetime.strptime(f'{str(dt0.date())}-{dt}', '%Y-%m-%d-%H') for dt in wm_hours] + [str(dt0)] + ser = pd.Series(mod_times, name='datetimes', dtype='datetime64[ns]') + ix_mod = ser.index[-1] # to get the nearest times + df = ser.sort_values().reset_index() + ix_mod_s = df[df['index']==ix_mod].index.item() + + ## case 1: index is not at the end + if ix_mod_s < df.index[-1]: + dt1, dt2 = df.loc[ix_mod_s-1, 'datetimes'], df.loc[ix_mod_s+1, 'datetimes'] + dt2 = dt2.to_pydatetime() + + ## case 2: close to midnight next day; add one to the date and start at midnight + else: + dt1 = df.loc[ix_mod_s-1, 'datetimes'] + + # handle leap years + days = 2 if (dt1.year % 4 == 0 and dt1.month == 2 and dt1.day == 28) else 1 + dt2 = datetime.combine(dt1 + relativedelta(days=days), dtime()) + + return dt1.to_pydatetime(), dt2 + From 89d5476202f94651d8da6011702166cf0d1b53a4 Mon Sep 17 00:00:00 2001 From: bbuzz31 Date: Mon, 23 Jan 2023 18:52:10 -0800 Subject: [PATCH 02/19] Implement weighting for wm cube --- tools/RAiDER/checkArgs.py | 1 - tools/RAiDER/cli/raider.py | 1 - tools/RAiDER/delay.py | 29 +++++++++++++++++++++++------ 3 files changed, 23 insertions(+), 8 deletions(-) diff --git a/tools/RAiDER/checkArgs.py b/tools/RAiDER/checkArgs.py index 46464164c..661ca4151 100644 --- a/tools/RAiDER/checkArgs.py +++ b/tools/RAiDER/checkArgs.py @@ -18,7 +18,6 @@ from RAiDER.logger import logger -## put it in here def checkArgs(args): ''' Helper fcn for checking argument compatibility and returns the diff --git a/tools/RAiDER/cli/raider.py b/tools/RAiDER/cli/raider.py index f33027ef4..f514a8df8 100644 --- a/tools/RAiDER/cli/raider.py +++ b/tools/RAiDER/cli/raider.py @@ -421,7 +421,6 @@ def calcDelaysInterp(iargs=None): wet_delay, hydro_delay = avg_delays(wetDelays, hydroDelays, t, time1, time2) - # Write the delays to file # Different options depending on the inputs if los.is_Projected(): diff --git a/tools/RAiDER/delay.py b/tools/RAiDER/delay.py index d4f5837c4..0ef6ad470 100755 --- a/tools/RAiDER/delay.py +++ b/tools/RAiDER/delay.py @@ -414,12 +414,29 @@ def _build_cube_ray( return outputArrs -def avg_delays(wetDelays, hydroDelays, t0, t1, t2): +def avg_delays(wetDelays, hydroDelays, dtref, dt1, dt2): + import xarray as xr if hydroDelays[0] is None: - # do something with wet - pass + ## difference in seconds between the model dates and the reference date + wgts = xr.DataArray([(dtref-dt1).seconds, (dt2-dtref).seconds], dims=['time'], coords={'time': [1,2]}) + ds_out_wet = wetDelays[0].copy() # store the updated data in here + for kind in 'hydro wet'.split(): + da1 = wetDelays[0][kind].assign_coords(time=1) + da2 = wetDelays[1][kind].assign_coords(time=2) + da = xr.concat([da1, da2], 'time').weighted(wgts) + da_mu = da.mean('time') + ds_out_wet[kind] = da_mu + + ds_out_wet.attrs['source'] = [wetDelays[i].attrs['source'] for i in range(2)] + ds_out_hydro = None + ## may be useful for testing + # da_sec1 = da_sec1.isel(z=0, y=slice(0,1), x=slice(0, 2)) + # da_sec2 = da_sec2.isel(z=0, y=slice(0,1), x=slice(0, 2)) + # np.average(da_sec.data, weights=wgts_sec.data, axis=0) ## equals weighted xarray + # dat = np.stack([da_sec1.data, da_sec2.data]).mean(0) ## equivalent to xr concat + else: - # do something with both - pass + for kind in 'hydro wet'.split(): + raise Exception('Not supported yet.') - return wetDelays[0], hydroDelays[0] \ No newline at end of file + return ds_out_wet, ds_out_hydro \ No newline at end of file From 10e7bfd69736e3ddf9ada56be7c9974569e956a2 Mon Sep 17 00:00:00 2001 From: bbuzz31 Date: Tue, 24 Jan 2023 09:37:54 -0800 Subject: [PATCH 03/19] Implement a test for time interp --- tools/RAiDER/cli/__main__.py | 2 +- tools/RAiDER/delay.py | 9 +++++---- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/tools/RAiDER/cli/__main__.py b/tools/RAiDER/cli/__main__.py index fd863cd3b..1bf62e118 100644 --- a/tools/RAiDER/cli/__main__.py +++ b/tools/RAiDER/cli/__main__.py @@ -11,7 +11,7 @@ def main(): formatter_class=argparse.ArgumentDefaultsHelpFormatter ) parser.add_argument( - '++process', choices=['calcDelays', 'calcDelaysInterp' 'downloadGNSS', 'calcDelaysGUNW'], + '++process', choices=['calcDelays', 'calcDelaysInterp', 'downloadGNSS', 'calcDelaysGUNW'], default='calcDelays', help='Select the entrypoint to use' ) diff --git a/tools/RAiDER/delay.py b/tools/RAiDER/delay.py index 0ef6ad470..7ef86438e 100755 --- a/tools/RAiDER/delay.py +++ b/tools/RAiDER/delay.py @@ -429,10 +429,11 @@ def avg_delays(wetDelays, hydroDelays, dtref, dt1, dt2): ds_out_wet.attrs['source'] = [wetDelays[i].attrs['source'] for i in range(2)] ds_out_hydro = None - ## may be useful for testing - # da_sec1 = da_sec1.isel(z=0, y=slice(0,1), x=slice(0, 2)) - # da_sec2 = da_sec2.isel(z=0, y=slice(0,1), x=slice(0, 2)) - # np.average(da_sec.data, weights=wgts_sec.data, axis=0) ## equals weighted xarray + ## useful for testing + # da1_crop = da1.isel(z=0, y=slice(0,1), x=slice(0, 2)) + # da2_crop = da2.isel(z=0, y=slice(0,1), x=slice(0, 2)) + + # np.average(da_1_cr.data, weights=wgts_sec.data, axis=0) ## equals weighted xarray # dat = np.stack([da_sec1.data, da_sec2.data]).mean(0) ## equivalent to xr concat else: From bba597313c58d77d0f2da4004a6d20f181bfbcb3 Mon Sep 17 00:00:00 2001 From: bbuzz31 Date: Tue, 24 Jan 2023 10:04:30 -0800 Subject: [PATCH 04/19] add another test for the weighting calculation --- test/test_temporal_interpolate.py | 186 ++++++++++++++++++++++++++++++ tools/RAiDER/delay.py | 6 - 2 files changed, 186 insertions(+), 6 deletions(-) create mode 100644 test/test_temporal_interpolate.py diff --git a/test/test_temporal_interpolate.py b/test/test_temporal_interpolate.py new file mode 100644 index 000000000..3f774c749 --- /dev/null +++ b/test/test_temporal_interpolate.py @@ -0,0 +1,186 @@ +import glob +import os +import pytest +import shutil +import subprocess +import yaml + +import numpy as np +import pandas as pd +import xarray as xr + +from test import TEST_DIR, WM + +import RAiDER +from RAiDER.logger import logger + +wm = 'ERA-5' if WM == 'ERA5' else WM + + +def makeLatLonGrid(bbox, reg, out_dir, spacing=0.1): + """ Make lat lons at a specified spacing """ + S, N, W, E = bbox + lat_st, lat_en = S, N + lon_st, lon_en = W, E + + lats = np.arange(lat_st, lat_en, spacing) + lons = np.arange(lon_st, lon_en, spacing) + Lat, Lon = np.meshgrid(lats, lons) + da_lat = xr.DataArray(Lat.T, name='data', coords={'lon': lons, 'lat': lats}, dims='lat lon'.split()) + da_lon = xr.DataArray(Lon.T, name='data', coords={'lon': lons, 'lat': lats}, dims='lat lon'.split()) + + dst_lat = os.path.join(out_dir, f'lat_{reg}.nc') + dst_lon = os.path.join(out_dir, f'lon_{reg}.nc') + da_lat.to_netcdf(dst_lat) + da_lon.to_netcdf(dst_lon) + + return dst_lat, dst_lon + + +def update_yaml(dct_cfg:dict, dst:str='temp.yaml'): + """ Write a new yaml file from a dictionary. + + Updates parameters in the default 'raider.yaml' file. + Each key:value pair will in 'dct_cfg' will overwrite that in the default + """ + + template_file = os.path.join( + os.path.dirname(RAiDER.__file__), 'cli', 'raider.yaml') + + with open(template_file, 'r') as f: + try: + params = yaml.safe_load(f) + except yaml.YAMLError as exc: + print(exc) + raise ValueError(f'Something is wrong with the yaml file {template_file}') + + params = {**params, **dct_cfg} + + with open(dst, 'w') as fh: + yaml.safe_dump(params, fh, default_flow_style=False) + + return dst + + +def test_cube_timemean(): + """ Test the mean interpolation by computing cube delays at 12 noon (exact) and 11 AM / 1PM for GMAO """ + SCENARIO_DIR = os.path.join(TEST_DIR, "INTERP_TIME") + os.makedirs(SCENARIO_DIR, exist_ok=True) + ## make the lat lon grid + S, N, W, E = 34, 35, -117, -116 + date = 20200130 + hr1, hr2 = 12, 15 + + grp = { + 'date_group': {'date_start': date}, + 'weather_model': 'GMAO', + 'aoi_group': {'bounding_box': [S, N, W, E]}, + 'runtime_group': {'output_directory': SCENARIO_DIR}, + } + + + ## run raider original for two exact weather model times + for hr in [hr1, hr2]: + grp['time_group'] = {'time': f'{hr}:00:00'} + ## generate the default template file and overwrite it with new parms + cfg = update_yaml(grp) + + ## run raider for the default date + cmd = f'raider.py {cfg}' + proc = subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True) + assert np.isclose(proc.returncode, 0) + + ## run interpolation in the middle of the two + grp['time_group'] = {'time': f'13:30:00'} + cfg = update_yaml(grp) + + cmd = f'raider.py ++process calcDelaysInterp {cfg}' + proc = subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True) + assert np.isclose(proc.returncode, 0) + + + wm = 'GMAO' + with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{wm}_tropo_{date}T{hr1}0000_ztd.nc')) as ds: + da1_tot = ds['wet'] + ds['hydro'] + + with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{wm}_tropo_{date}T{hr2}0000_ztd.nc')) as ds: + da2_tot = ds['wet'] + ds['hydro'] + + with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{wm}_tropo_{date}T133000_ztd.nc')) as ds: + da_interp_tot = ds['wet'] + ds['hydro'] + + da_mu = (da1_tot + da2_tot) / 2 + assert np.allclose(da_mu, da_interp_tot) + + + # # Clean up files + shutil.rmtree(SCENARIO_DIR) + [os.remove(f) for f in glob.glob(f'{wm}*')] + os.remove('temp.yaml') + + return + +def test_cube_weighting(): + from datetime import datetime + SCENARIO_DIR = os.path.join(TEST_DIR, "INTERP_TIME") + os.makedirs(SCENARIO_DIR, exist_ok=True) + ## make the lat lon grid + S, N, W, E = 34, 35, -117, -116 + date = 20200130 + hr1, hr2 = 12, 15 + t_ref = '12:05:00' + + grp = { + 'date_group': {'date_start': date}, + 'weather_model': 'GMAO', + 'aoi_group': {'bounding_box': [S, N, W, E]}, + 'runtime_group': {'output_directory': SCENARIO_DIR}, + } + + + ## run raider original for two exact weather model times + for hr in [hr1, hr2]: + grp['time_group'] = {'time': f'{hr}:00:00'} + ## generate the default template file and overwrite it with new parms + cfg = update_yaml(grp) + + ## run raider for the default date + cmd = f'raider.py {cfg}' + proc = subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True) + assert np.isclose(proc.returncode, 0) + + ## run interpolation very near the first + grp['time_group'] = {'time': t_ref} + cfg = update_yaml(grp) + + cmd = f'raider.py ++process calcDelaysInterp {cfg}' + proc = subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True) + + ## double check on weighting + + wm = 'GMAO' + with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{wm}_tropo_{date}T{hr1}0000_ztd.nc')) as ds: + da1_tot = ds['wet'] + ds['hydro'] + + with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{wm}_tropo_{date}T{hr2}0000_ztd.nc')) as ds: + da2_tot = ds['wet'] + ds['hydro'] + + with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{wm}_tropo_{date}T{t_ref.replace(":", "")}_ztd.nc')) as ds: + da_interp_tot = ds['wet'] + ds['hydro'] + + assert None + dt1 = datetime.strptime(f'{date}{hr1}', '%Y%m%d%H') + dt2 = datetime.strptime(f'{date}{hr2}', '%Y%m%d%H') + dt_ref = datetime.strptime(f'{date}{t_ref}', '%Y%m%d%H:%M:%S') + + wgts = (dt_ref-dt1).seconds, (dt2-dt_ref).seconds + logger.info ('Weights: %s', wgts) + logger.info ('Tstart: %s, Tend: %s, Tref: %s', dt1, dt1, dt_ref) + da1_crop = da1_tot.isel(z=0, y=slice(0,1), x=slice(0, 2)) + da2_crop = da2_tot.isel(z=0, y=slice(0,1), x=slice(0, 2)) + da_out_crop = da_interp_tot.isel(z=0, y=slice(0,1), x=slice(0,2)) + + dat = np.vstack([da1_crop.data, da2_crop.data]) + np.allclose(da_out_crop, np.average(dat, weights=wgts, axis=0)) + return + diff --git a/tools/RAiDER/delay.py b/tools/RAiDER/delay.py index 7ef86438e..0c9a90f45 100755 --- a/tools/RAiDER/delay.py +++ b/tools/RAiDER/delay.py @@ -429,12 +429,6 @@ def avg_delays(wetDelays, hydroDelays, dtref, dt1, dt2): ds_out_wet.attrs['source'] = [wetDelays[i].attrs['source'] for i in range(2)] ds_out_hydro = None - ## useful for testing - # da1_crop = da1.isel(z=0, y=slice(0,1), x=slice(0, 2)) - # da2_crop = da2.isel(z=0, y=slice(0,1), x=slice(0, 2)) - - # np.average(da_1_cr.data, weights=wgts_sec.data, axis=0) ## equals weighted xarray - # dat = np.stack([da_sec1.data, da_sec2.data]).mean(0) ## equivalent to xr concat else: for kind in 'hydro wet'.split(): From 46b2cf79af80dbe5e91f2d8875259a0e5d74a317 Mon Sep 17 00:00:00 2001 From: bbuzz31 Date: Tue, 24 Jan 2023 10:18:41 -0800 Subject: [PATCH 05/19] correctly assign weights --- test/test_temporal_interpolate.py | 30 ++++++++++++++---------------- tools/RAiDER/delay.py | 3 ++- 2 files changed, 16 insertions(+), 17 deletions(-) diff --git a/test/test_temporal_interpolate.py b/test/test_temporal_interpolate.py index 3f774c749..e2329ed8c 100644 --- a/test/test_temporal_interpolate.py +++ b/test/test_temporal_interpolate.py @@ -14,7 +14,7 @@ import RAiDER from RAiDER.logger import logger -wm = 'ERA-5' if WM == 'ERA5' else WM +WM = 'GMAO' def makeLatLonGrid(bbox, reg, out_dir, spacing=0.1): @@ -63,7 +63,7 @@ def update_yaml(dct_cfg:dict, dst:str='temp.yaml'): def test_cube_timemean(): - """ Test the mean interpolation by computing cube delays at 12 noon (exact) and 11 AM / 1PM for GMAO """ + """ Test the mean interpolation by computing cube delays at 1:30PM vs mean of 12 PM / 3PM for GMAO """ SCENARIO_DIR = os.path.join(TEST_DIR, "INTERP_TIME") os.makedirs(SCENARIO_DIR, exist_ok=True) ## make the lat lon grid @@ -73,7 +73,7 @@ def test_cube_timemean(): grp = { 'date_group': {'date_start': date}, - 'weather_model': 'GMAO', + 'weather_model': WM, 'aoi_group': {'bounding_box': [S, N, W, E]}, 'runtime_group': {'output_directory': SCENARIO_DIR}, } @@ -99,14 +99,13 @@ def test_cube_timemean(): assert np.isclose(proc.returncode, 0) - wm = 'GMAO' - with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{wm}_tropo_{date}T{hr1}0000_ztd.nc')) as ds: + with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{hr1}0000_ztd.nc')) as ds: da1_tot = ds['wet'] + ds['hydro'] - with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{wm}_tropo_{date}T{hr2}0000_ztd.nc')) as ds: + with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{hr2}0000_ztd.nc')) as ds: da2_tot = ds['wet'] + ds['hydro'] - with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{wm}_tropo_{date}T133000_ztd.nc')) as ds: + with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T133000_ztd.nc')) as ds: da_interp_tot = ds['wet'] + ds['hydro'] da_mu = (da1_tot + da2_tot) / 2 @@ -115,12 +114,13 @@ def test_cube_timemean(): # # Clean up files shutil.rmtree(SCENARIO_DIR) - [os.remove(f) for f in glob.glob(f'{wm}*')] + [os.remove(f) for f in glob.glob(f'{WM}*')] os.remove('temp.yaml') return def test_cube_weighting(): + """ Test the weighting by comparing a small crop with numpy directly """ from datetime import datetime SCENARIO_DIR = os.path.join(TEST_DIR, "INTERP_TIME") os.makedirs(SCENARIO_DIR, exist_ok=True) @@ -132,7 +132,7 @@ def test_cube_weighting(): grp = { 'date_group': {'date_start': date}, - 'weather_model': 'GMAO', + 'weather_model': WM, 'aoi_group': {'bounding_box': [S, N, W, E]}, 'runtime_group': {'output_directory': SCENARIO_DIR}, } @@ -158,22 +158,20 @@ def test_cube_weighting(): ## double check on weighting - wm = 'GMAO' - with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{wm}_tropo_{date}T{hr1}0000_ztd.nc')) as ds: + with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{hr1}0000_ztd.nc')) as ds: da1_tot = ds['wet'] + ds['hydro'] - with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{wm}_tropo_{date}T{hr2}0000_ztd.nc')) as ds: + with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{hr2}0000_ztd.nc')) as ds: da2_tot = ds['wet'] + ds['hydro'] - with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{wm}_tropo_{date}T{t_ref.replace(":", "")}_ztd.nc')) as ds: + with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{t_ref.replace(":", "")}_ztd.nc')) as ds: da_interp_tot = ds['wet'] + ds['hydro'] - assert None dt1 = datetime.strptime(f'{date}{hr1}', '%Y%m%d%H') dt2 = datetime.strptime(f'{date}{hr2}', '%Y%m%d%H') dt_ref = datetime.strptime(f'{date}{t_ref}', '%Y%m%d%H:%M:%S') - wgts = (dt_ref-dt1).seconds, (dt2-dt_ref).seconds + wgts = np.array([(dt_ref-dt1).seconds, (dt2-dt_ref).seconds]) logger.info ('Weights: %s', wgts) logger.info ('Tstart: %s, Tend: %s, Tref: %s', dt1, dt1, dt_ref) da1_crop = da1_tot.isel(z=0, y=slice(0,1), x=slice(0, 2)) @@ -181,6 +179,6 @@ def test_cube_weighting(): da_out_crop = da_interp_tot.isel(z=0, y=slice(0,1), x=slice(0,2)) dat = np.vstack([da1_crop.data, da2_crop.data]) - np.allclose(da_out_crop, np.average(dat, weights=wgts, axis=0)) + np.allclose(da_out_crop, np.average(dat, weights=1/wgts, axis=0)) return diff --git a/tools/RAiDER/delay.py b/tools/RAiDER/delay.py index 0c9a90f45..79df68b4e 100755 --- a/tools/RAiDER/delay.py +++ b/tools/RAiDER/delay.py @@ -418,7 +418,8 @@ def avg_delays(wetDelays, hydroDelays, dtref, dt1, dt2): import xarray as xr if hydroDelays[0] is None: ## difference in seconds between the model dates and the reference date - wgts = xr.DataArray([(dtref-dt1).seconds, (dt2-dtref).seconds], dims=['time'], coords={'time': [1,2]}) + dsecs = np.array([(dtref-dt1).seconds, (dt2-dtref).seconds]) + wgts = xr.DataArray(1/dsecs, dims=['time'], coords={'time': [1,2]}) ds_out_wet = wetDelays[0].copy() # store the updated data in here for kind in 'hydro wet'.split(): da1 = wetDelays[0][kind].assign_coords(time=1) From b730d5ab1cf6c0871340c9c9ec840a1ab1bbbd62 Mon Sep 17 00:00:00 2001 From: bbuzz31 Date: Tue, 24 Jan 2023 10:58:33 -0800 Subject: [PATCH 06/19] support ERA5 tinterp test --- test/test_temporal_interpolate.py | 16 +++++++--------- tools/RAiDER/cli/validators.py | 2 +- tools/RAiDER/models/ecmwf.py | 7 ++++++- tools/RAiDER/models/gmao.py | 2 +- tools/RAiDER/utilFcns.py | 2 +- 5 files changed, 16 insertions(+), 13 deletions(-) diff --git a/test/test_temporal_interpolate.py b/test/test_temporal_interpolate.py index e2329ed8c..52d0531c7 100644 --- a/test/test_temporal_interpolate.py +++ b/test/test_temporal_interpolate.py @@ -14,9 +14,6 @@ import RAiDER from RAiDER.logger import logger -WM = 'GMAO' - - def makeLatLonGrid(bbox, reg, out_dir, spacing=0.1): """ Make lat lons at a specified spacing """ S, N, W, E = bbox @@ -69,7 +66,8 @@ def test_cube_timemean(): ## make the lat lon grid S, N, W, E = 34, 35, -117, -116 date = 20200130 - hr1, hr2 = 12, 15 + dct_hrs = {'GMAO': [12, 15, '13:30:00'], 'ERA5': [13, 14, '13:30:00']} + hr1, hr2, ti = dct_hrs[WM] grp = { 'date_group': {'date_start': date}, @@ -91,7 +89,7 @@ def test_cube_timemean(): assert np.isclose(proc.returncode, 0) ## run interpolation in the middle of the two - grp['time_group'] = {'time': f'13:30:00'} + grp['time_group'] = {'time': ti} cfg = update_yaml(grp) cmd = f'raider.py ++process calcDelaysInterp {cfg}' @@ -105,7 +103,7 @@ def test_cube_timemean(): with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{hr2}0000_ztd.nc')) as ds: da2_tot = ds['wet'] + ds['hydro'] - with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T133000_ztd.nc')) as ds: + with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{ti.replace(":", "")}_ztd.nc')) as ds: da_interp_tot = ds['wet'] + ds['hydro'] da_mu = (da1_tot + da2_tot) / 2 @@ -127,8 +125,8 @@ def test_cube_weighting(): ## make the lat lon grid S, N, W, E = 34, 35, -117, -116 date = 20200130 - hr1, hr2 = 12, 15 - t_ref = '12:05:00' + dct_hrs = {'GMAO': [12, 15, '12:05:00'], 'ERA5': [13, 14, '13:05:00']} + hr1, hr2, t_ref = dct_hrs[WM] grp = { 'date_group': {'date_start': date}, @@ -173,7 +171,7 @@ def test_cube_weighting(): wgts = np.array([(dt_ref-dt1).seconds, (dt2-dt_ref).seconds]) logger.info ('Weights: %s', wgts) - logger.info ('Tstart: %s, Tend: %s, Tref: %s', dt1, dt1, dt_ref) + logger.info ('Tstart: %s, Tend: %s, Tref: %s', dt1, dt2, dt_ref) da1_crop = da1_tot.isel(z=0, y=slice(0,1), x=slice(0, 2)) da2_crop = da2_tot.isel(z=0, y=slice(0,1), x=slice(0, 2)) da_out_crop = da_interp_tot.isel(z=0, y=slice(0,1), x=slice(0,2)) diff --git a/tools/RAiDER/cli/validators.py b/tools/RAiDER/cli/validators.py index 6d6dc0b93..b5d999148 100755 --- a/tools/RAiDER/cli/validators.py +++ b/tools/RAiDER/cli/validators.py @@ -114,7 +114,7 @@ def get_heights(args, out, station_file, bounding_box=None): out['height_levels'] = np.array([float(ll) for ll in l]) if np.any(out['height_levels'] < 0): logger.warning('Weather model only extends to the surface topography; ' - 'height levels below the topography will be interpolated from the surface' + 'height levels below the topography will be interpolated from the surface ' 'and may be inaccurate.') return out diff --git a/tools/RAiDER/models/ecmwf.py b/tools/RAiDER/models/ecmwf.py index 532abdd93..193d2179b 100755 --- a/tools/RAiDER/models/ecmwf.py +++ b/tools/RAiDER/models/ecmwf.py @@ -164,6 +164,8 @@ def _get_from_ecmwf(self, lat_min, lat_max, lat_step, lon_min, lon_max, server = ecmwfapi.ECMWFDataServer() corrected_date = util.round_date(time, datetime.timedelta(hours=6)) + # if not time1 == time: + # logger.warning('Rounded given hour from %d to %d', time1.hour, time.hour) server.retrieve({ "class": self._classname, # ERA-Interim @@ -204,6 +206,7 @@ def _get_from_cds( acqTime, outname ): + """ Used for ERA5 """ import cdsapi c = cdsapi.Client(verify=0) @@ -219,7 +222,7 @@ def _get_from_cds( # round to the closest legal time corrected_date = util.round_date(acqTime, datetime.timedelta(hours=self._time_res)) if not corrected_date == acqTime: - logger.warning('Rounded given datetime from %s to %s', acqTime, corrected_date) + logger.warning('Rounded given hour from %d to %d', acqTime.hour, corrected_date.hour) # I referenced https://confluence.ecmwf.int/display/CKB/How+to+download+ERA5 dataDict = { @@ -247,7 +250,9 @@ def _get_from_cds( logger.exception(e) raise Exception + def _download_ecmwf(self, lat_min, lat_max, lat_step, lon_min, lon_max, lon_step, time, out): + """ Used for HRES """ from ecmwfapi import ECMWFService server = ECMWFService("mars") diff --git a/tools/RAiDER/models/gmao.py b/tools/RAiDER/models/gmao.py index ad24655ed..ec7e21245 100755 --- a/tools/RAiDER/models/gmao.py +++ b/tools/RAiDER/models/gmao.py @@ -37,7 +37,7 @@ def __init__(self): self._k2 = 0.233 # [K/Pa] self._k3 = 3.75e3 # [K^2/Pa] - self._time_res = 1 + self._time_res = 3 # horizontal grid spacing self._lat_res = 0.25 diff --git a/tools/RAiDER/utilFcns.py b/tools/RAiDER/utilFcns.py index 533d848e8..ff16600e7 100755 --- a/tools/RAiDER/utilFcns.py +++ b/tools/RAiDER/utilFcns.py @@ -1143,7 +1143,7 @@ def get_nearest_wmtimes(dt0, wm): from dateutil.relativedelta import relativedelta from datetime import time as dtime ## hourly time step availability - not sure about ECMWF / HRRR - wm_hour = {'GMAO': 3, 'ERA5':1, 'ERA5T': 1, 'HRES':1, 'HRRR': 1, 'NCMR': '?'} + wm_hour = {'GMAO': 3, 'ERA5':1, 'ERA-5': 1, 'ERA5T': 1, 'HRES':6, 'HRRR': 1, 'NCMR': '?'} wm_hours = np.arange(0, 24, wm_hour[wm]) dates = [] From 441288489cf2bee890cbb517d558d4fcef163324 Mon Sep 17 00:00:00 2001 From: bbuzz31 Date: Tue, 24 Jan 2023 11:01:04 -0800 Subject: [PATCH 07/19] few more warnings --- tools/RAiDER/models/ecmwf.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tools/RAiDER/models/ecmwf.py b/tools/RAiDER/models/ecmwf.py index 193d2179b..fde5fef3b 100755 --- a/tools/RAiDER/models/ecmwf.py +++ b/tools/RAiDER/models/ecmwf.py @@ -164,8 +164,8 @@ def _get_from_ecmwf(self, lat_min, lat_max, lat_step, lon_min, lon_max, server = ecmwfapi.ECMWFDataServer() corrected_date = util.round_date(time, datetime.timedelta(hours=6)) - # if not time1 == time: - # logger.warning('Rounded given hour from %d to %d', time1.hour, time.hour) + if not corrected_date == time: + logger.warning('Rounded given hour from %d to %d', time.hour, corrected_date.hour) server.retrieve({ "class": self._classname, # ERA-Interim @@ -260,7 +260,7 @@ def _download_ecmwf(self, lat_min, lat_max, lat_step, lon_min, lon_max, lon_step # round to the closest legal time corrected_date = util.round_date(time, datetime.timedelta(hours=self._time_res)) if not corrected_date == time: - logger.warning('Rounded given datetime from %s to %s', time, corrected_date) + logger.warning('Rounded given hour from %d to %d', time.hour, corrected_date.hour) if self._model_level_type == 'ml': param = "129/130/133/152" From d756bab60a8412e503eca6f3042b7a458f66a697 Mon Sep 17 00:00:00 2001 From: bbuzz31 Date: Tue, 24 Jan 2023 14:34:40 -0800 Subject: [PATCH 08/19] Embed temporal interpolation into delay.py --- pyproject.toml | 1 - test/test_temporal_interpolate.py | 33 ++++-- tools/RAiDER/aria/prepFromGUNW.py | 2 +- tools/RAiDER/cli/__init__.py | 1 + tools/RAiDER/cli/__main__.py | 4 +- tools/RAiDER/cli/raider.py | 181 ++---------------------------- tools/RAiDER/cli/raider.yaml | 1 + 7 files changed, 37 insertions(+), 186 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 57a61db1b..3a69e4d44 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -43,7 +43,6 @@ repository = "https://github.com/dbekaert/RAiDER" [project.scripts] "raider.py" = "RAiDER.cli.__main__:main" "calcDelays.py" = "RAiDER.cli.__main__:calcDelays" -"calcDelaysInterp.py" = "RAiDER.cli.__main__:calcDelaysInterp" "calcDelaysGUNW.py" = "RAiDER.cli.__main__:calcDelaysGUNW" "raiderDownloadGNSS.py" = "RAiDER.cli.__main__:downloadGNSS" "downloadGNSS.py" = "RAiDER.cli.__main__:downloadGNSS" diff --git a/test/test_temporal_interpolate.py b/test/test_temporal_interpolate.py index 52d0531c7..5b132b7ea 100644 --- a/test/test_temporal_interpolate.py +++ b/test/test_temporal_interpolate.py @@ -14,6 +14,7 @@ import RAiDER from RAiDER.logger import logger + def makeLatLonGrid(bbox, reg, out_dir, spacing=0.1): """ Make lat lons at a specified spacing """ S, N, W, E = bbox @@ -89,10 +90,10 @@ def test_cube_timemean(): assert np.isclose(proc.returncode, 0) ## run interpolation in the middle of the two - grp['time_group'] = {'time': ti} + grp['time_group'] = {'time': ti, 'interpolate_wmtimes': True} cfg = update_yaml(grp) - cmd = f'raider.py ++process calcDelaysInterp {cfg}' + cmd = f'raider.py {cfg}' proc = subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True) assert np.isclose(proc.returncode, 0) @@ -110,13 +111,14 @@ def test_cube_timemean(): assert np.allclose(da_mu, da_interp_tot) - # # Clean up files + # Clean up files shutil.rmtree(SCENARIO_DIR) [os.remove(f) for f in glob.glob(f'{WM}*')] os.remove('temp.yaml') return + def test_cube_weighting(): """ Test the weighting by comparing a small crop with numpy directly """ from datetime import datetime @@ -126,7 +128,7 @@ def test_cube_weighting(): S, N, W, E = 34, 35, -117, -116 date = 20200130 dct_hrs = {'GMAO': [12, 15, '12:05:00'], 'ERA5': [13, 14, '13:05:00']} - hr1, hr2, t_ref = dct_hrs[WM] + hr1, hr2, ti = dct_hrs[WM] grp = { 'date_group': {'date_start': date}, @@ -148,10 +150,10 @@ def test_cube_weighting(): assert np.isclose(proc.returncode, 0) ## run interpolation very near the first - grp['time_group'] = {'time': t_ref} + grp['time_group'] = {'time': ti, 'interpolate_wmtimes': True} cfg = update_yaml(grp) - cmd = f'raider.py ++process calcDelaysInterp {cfg}' + cmd = f'raider.py {cfg}' proc = subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True) ## double check on weighting @@ -162,21 +164,30 @@ def test_cube_weighting(): with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{hr2}0000_ztd.nc')) as ds: da2_tot = ds['wet'] + ds['hydro'] - with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{t_ref.replace(":", "")}_ztd.nc')) as ds: + with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{ti.replace(":", "")}_ztd.nc')) as ds: da_interp_tot = ds['wet'] + ds['hydro'] dt1 = datetime.strptime(f'{date}{hr1}', '%Y%m%d%H') dt2 = datetime.strptime(f'{date}{hr2}', '%Y%m%d%H') - dt_ref = datetime.strptime(f'{date}{t_ref}', '%Y%m%d%H:%M:%S') + dt_ref = datetime.strptime(f'{date}{ti}', '%Y%m%d%H:%M:%S') wgts = np.array([(dt_ref-dt1).seconds, (dt2-dt_ref).seconds]) - logger.info ('Weights: %s', wgts) - logger.info ('Tstart: %s, Tend: %s, Tref: %s', dt1, dt2, dt_ref) da1_crop = da1_tot.isel(z=0, y=slice(0,1), x=slice(0, 2)) da2_crop = da2_tot.isel(z=0, y=slice(0,1), x=slice(0, 2)) da_out_crop = da_interp_tot.isel(z=0, y=slice(0,1), x=slice(0,2)) dat = np.vstack([da1_crop.data, da2_crop.data]) - np.allclose(da_out_crop, np.average(dat, weights=1/wgts, axis=0)) + + logger.info ('Tstart: %s, Tend: %s, Tref: %s', dt1, dt2, dt_ref) + logger.info ('Weights: %s', wgts) + logger.info ('Data from two dates: %s', dat) + logger.info ('Weighted mean: %s', da_out_crop.data) + assert np.allclose(da_out_crop, np.average(dat, weights=1/wgts, axis=0)) + + # Clean up files + shutil.rmtree(SCENARIO_DIR) + [os.remove(f) for f in glob.glob(f'{WM}*')] + os.remove('temp.yaml') + return diff --git a/tools/RAiDER/aria/prepFromGUNW.py b/tools/RAiDER/aria/prepFromGUNW.py index e50161fe7..b0a465dd9 100644 --- a/tools/RAiDER/aria/prepFromGUNW.py +++ b/tools/RAiDER/aria/prepFromGUNW.py @@ -245,7 +245,7 @@ def main(args): 'aoi_group' : {'bounding_box': GUNWObj.SNWE}, 'height_group' : {'height_levels': GUNWObj.heights}, 'date_group': {'date_list': str(GUNWObj.dates)}, - 'time_group': {'time': GUNWObj.ref_time}, + 'time_group': {'time': GUNWObj.ref_time, 'interpolate_wmtimes': True}, 'los_group' : {'ray_trace': True, 'orbit_file': GUNWObj.OrbitFiles, 'wavelength': GUNWObj.wavelength, diff --git a/tools/RAiDER/cli/__init__.py b/tools/RAiDER/cli/__init__.py index 8a0715ae1..bbc0ebbc6 100644 --- a/tools/RAiDER/cli/__init__.py +++ b/tools/RAiDER/cli/__init__.py @@ -39,5 +39,6 @@ class AttributeDict(dict): output_directory='.', weather_model_directory=None, output_projection='EPSG:4326', + interpolate_wmtimes=False ) ) diff --git a/tools/RAiDER/cli/__main__.py b/tools/RAiDER/cli/__main__.py index 1bf62e118..3b195d957 100644 --- a/tools/RAiDER/cli/__main__.py +++ b/tools/RAiDER/cli/__main__.py @@ -2,7 +2,7 @@ import sys from importlib.metadata import entry_points -from RAiDER.cli.raider import calcDelays, calcDelaysInterp, downloadGNSS, calcDelaysGUNW +from RAiDER.cli.raider import calcDelays, downloadGNSS, calcDelaysGUNW def main(): @@ -11,7 +11,7 @@ def main(): formatter_class=argparse.ArgumentDefaultsHelpFormatter ) parser.add_argument( - '++process', choices=['calcDelays', 'calcDelaysInterp', 'downloadGNSS', 'calcDelaysGUNW'], + '++process', choices=['calcDelays', 'downloadGNSS', 'calcDelaysGUNW'], default='calcDelays', help='Select the entrypoint to use' ) diff --git a/tools/RAiDER/cli/raider.py b/tools/RAiDER/cli/raider.py index f514a8df8..287b2cd2a 100644 --- a/tools/RAiDER/cli/raider.py +++ b/tools/RAiDER/cli/raider.py @@ -121,10 +121,10 @@ def drop_nans(d): def calcDelays(iargs=None): """ Parse command line arguments using argparse. """ import RAiDER - from RAiDER.delay import tropo_delay + from RAiDER.delay import tropo_delay, avg_delays from RAiDER.checkArgs import checkArgs from RAiDER.processWM import prepareWeatherModel - from RAiDER.utilFcns import writeDelays + from RAiDER.utilFcns import writeDelays, get_nearest_wmtimes examples = 'Examples of use:' \ '\n\t raider.py customTemplatefile.cfg' \ '\n\t raider.py -g' @@ -220,171 +220,9 @@ def calcDelays(iargs=None): logger.debug('Beginning weather model pre-processing') - try: - weather_model_file = prepareWeatherModel( - model, t, - ll_bounds=ll_bounds, # SNWE - wmLoc=params['weather_model_directory'], - makePlots=params['verbose'], - ) - except RuntimeError: - logger.exception("Date %s failed", t) - continue - - # dont process the delays for download only - if dl_only: - continue - - # Now process the delays - try: - wet_delay, hydro_delay = tropo_delay( - t, weather_model_file, aoi, los, - height_levels = params['height_levels'], - out_proj = params['output_projection'], - look_dir = params['look_dir'], - cube_spacing_m = params['cube_spacing_in_m'], - ) - except RuntimeError: - logger.exception("Date %s failed", t) - continue - - ########################################################### - # Write the delays to file - # Different options depending on the inputs - - if los.is_Projected(): - out_filename = w.replace("_ztd", "_std") - f = f.replace("_ztd", "_std") - elif los.ray_trace(): - out_filename = w.replace("_std", "_ray") - f = f.replace("_std", "_ray") - else: - out_filename = w - - if hydro_delay is None: - # means that a dataset was returned - ds = wet_delay - ext = os.path.splitext(out_filename)[1] - out_filename = out_filename.replace('wet', 'tropo') - - if ext not in ['.nc', '.h5']: - out_filename = f'{os.path.splitext(out_filename)[0]}.nc' - - - if out_filename.endswith(".nc"): - ds.to_netcdf(out_filename, mode="w") - elif out_filename.endswith(".h5"): - ds.to_netcdf(out_filename, engine="h5netcdf", invalid_netcdf=True) - logger.info('Wrote delays to: %s', out_filename) - - else: - if aoi.type() == 'station_file': - out_filename = f'{os.path.splitext(out_filename)[0]}.csv' - - if aoi.type() in ['station_file', 'radar_rasters', 'geocoded_file']: - writeDelays(aoi, wet_delay, hydro_delay, out_filename, f, outformat=params['raster_format']) - - wet_filenames.append(out_filename) - - return wet_filenames - - -def calcDelaysInterp(iargs=None): - """ - Same as calcDelays, except will calculate at nearest weather model hour then interpolate - """ - import RAiDER - from RAiDER.delay import tropo_delay, avg_delays - from RAiDER.checkArgs import checkArgs - from RAiDER.processWM import prepareWeatherModel - from RAiDER.utilFcns import writeDelays, get_nearest_wmtimes - examples = 'Examples of use:' \ - '\n\t raider.py ++process calcDelaysInterp customTemplatefile.cfg' \ - '\n\t calcDelaysInterp.py customTemplateFile' - - p = argparse.ArgumentParser( - description = - 'Command line interface for RAiDER processing with a configure file.' - 'Default options can be found by running: raider.py --generate_config', - epilog=examples, formatter_class=argparse.RawDescriptionHelpFormatter) - - p.add_argument( - 'customTemplateFile', nargs='?', - help='custom template with option settings.\n' + - "ignored if the default smallbaselineApp.cfg is input." - ) - - p.add_argument( - '--download_only', action='store_true', - help='only download a weather model.' - ) - - ## if not None, will replace first argument (customTemplateFile) - args = p.parse_args(args=iargs) - - # default input file - template_file = os.path.join(os.path.dirname(RAiDER.__file__), - 'cli', 'raider.yaml') - - # check: existence of input template files - if (not args.customTemplateFile - and not os.path.isfile(os.path.basename(template_file)) - and not args.generate_template): - msg = "No template file found! It requires that either:" - msg += "\n a custom template file, OR the default template " - msg += "\n file 'raider.yaml' exists in current directory." - - p.print_usage() - print(examples) - raise SystemExit(f'ERROR: {msg}') - - if args.customTemplateFile: - # check the existence - if not os.path.isfile(args.customTemplateFile): - raise FileNotFoundError(args.customTemplateFile) - - args.customTemplateFile = os.path.abspath(args.customTemplateFile) - else: - args.customTemplateFile = template_file - - # Read the template file - params = read_template_file(args.customTemplateFile) - - # Argument checking - params = checkArgs(params) - dl_only = True if params['download_only'] or args.download_only else False - - if not params.verbose: - logger.setLevel(logging.INFO) - - wet_filenames = [] - for t, w, f in zip( - params['date_list'], - params['wetFilenames'], - params['hydroFilenames'] - ): - - los = params['los'] - aoi = params['aoi'] - model = params['weather_model'] - - # add a buffer for raytracing - if los.ray_trace(): - ll_bounds = aoi.add_buffer(buffer=1) - else: - ll_bounds = aoi.add_buffer(buffer=1) - - ########################################################### - # weather model calculation - logger.debug('Starting to run the weather model calculation') - logger.debug(f'Date: {t.strftime("%Y%m%d")}') - logger.debug('Beginning weather model pre-processing') - - - - time1, time2 = get_nearest_wmtimes(t, model._Name) + times = get_nearest_wmtimes(t, model._Name) if params['interpolate_wmtimes'] else [t] wetDelays, hydroDelays = [], [] - for tt in [time1, time2]: + for tt in times: try: weather_model_file = prepareWeatherModel( model, tt, @@ -393,10 +231,9 @@ def calcDelaysInterp(iargs=None): makePlots=params['verbose'], ) except RuntimeError: - logger.exception("Date %s failed", tt) + logger.exception("Date %s failed", t) continue - # dont process the delays for download only if dl_only: continue @@ -414,14 +251,15 @@ def calcDelaysInterp(iargs=None): logger.exception("Date %s failed", t) continue - ########################################################### # store in memory wetDelays.append(wet_delay) hydroDelays.append(hydro_delay) - wet_delay, hydro_delay = avg_delays(wetDelays, hydroDelays, t, time1, time2) - + ########################################################### # Write the delays to file + if params['interpolate_wmtimes']: + wet_delay, hydro_delay = avg_delays(wetDelays, hydroDelays, t, times[0], times[1]) + # Different options depending on the inputs if los.is_Projected(): out_filename = w.replace("_ztd", "_std") @@ -459,6 +297,7 @@ def calcDelaysInterp(iargs=None): return wet_filenames + ## ------------------------------------------------------ downloadGNSSDelays.py def downloadGNSS(): """Parse command line arguments using argparse.""" diff --git a/tools/RAiDER/cli/raider.yaml b/tools/RAiDER/cli/raider.yaml index 4c84e8724..bd854adaf 100644 --- a/tools/RAiDER/cli/raider.yaml +++ b/tools/RAiDER/cli/raider.yaml @@ -58,6 +58,7 @@ date_group: time_group: time: end_time: + interpolate_wmtimes: false # Interpolate between weather model times if true, else use nearest time ########## 4. Area of Interest From e10ffff7f120e5166049abc0f8507ade5a531c19 Mon Sep 17 00:00:00 2001 From: Joseph H Kennedy Date: Tue, 31 Jan 2023 13:17:53 -0900 Subject: [PATCH 09/19] Update CHANGELOG.md --- CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index fcd99f784..141742227 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,6 +19,11 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). + Only orbit file is used (even if multiple specified) to minimize errors and ensure consistency over region + GUNW packaging is restructed to store SLC (ref and sec) wet and tropo delays rather than the differential +##[0.4.2] + +### New/Updated Features ++ `calcDelaysGUNW` allows processing with any supported weather model as listed in [`RAiDER.models.allowed.ALLOWED_MODELS`](https://github.com/dbekaert/RAiDER/blob/dev/tools/RAiDER/models/allowed.py). + ## [0.4.1] ### New/Updated Features From 089afc808eeaa4452b3c3714a1a6e8bf632ee192 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Tue, 14 Feb 2023 10:10:43 -0800 Subject: [PATCH 10/19] final fixes --- CHANGELOG.md | 1 - 1 file changed, 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 141742227..7bfe5cc2c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,7 +24,6 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ### New/Updated Features + `calcDelaysGUNW` allows processing with any supported weather model as listed in [`RAiDER.models.allowed.ALLOWED_MODELS`](https://github.com/dbekaert/RAiDER/blob/dev/tools/RAiDER/models/allowed.py). -## [0.4.1] ### New/Updated Features + Reorder target points for intersection From 64cc5448855f64e05ee3857690ff7d733dcd1b68 Mon Sep 17 00:00:00 2001 From: Jeremy Maurer Date: Tue, 28 Feb 2023 22:52:30 -0600 Subject: [PATCH 11/19] move delay calc out of time loop --- tools/RAiDER/cli/raider.py | 67 +++++++++++++++++++++----------------- 1 file changed, 37 insertions(+), 30 deletions(-) diff --git a/tools/RAiDER/cli/raider.py b/tools/RAiDER/cli/raider.py index 287b2cd2a..9c97af33b 100644 --- a/tools/RAiDER/cli/raider.py +++ b/tools/RAiDER/cli/raider.py @@ -221,44 +221,51 @@ def calcDelays(iargs=None): times = get_nearest_wmtimes(t, model._Name) if params['interpolate_wmtimes'] else [t] - wetDelays, hydroDelays = [], [] + wfiles = [] for tt in times: try: - weather_model_file = prepareWeatherModel( - model, tt, - ll_bounds=ll_bounds, # SNWE - wmLoc=params['weather_model_directory'], - makePlots=params['verbose'], + wfiles.append( + prepareWeatherModel( + model, tt, + ll_bounds=ll_bounds, # SNWE + wmLoc=params['weather_model_directory'], + makePlots=params['verbose'], + ) ) except RuntimeError: logger.exception("Date %s failed", t) continue - # dont process the delays for download only - if dl_only: - continue - - # Now process the delays - try: - wet_delay, hydro_delay = tropo_delay( - t, weather_model_file, aoi, los, - height_levels = params['height_levels'], - out_proj = params['output_projection'], - look_dir = params['look_dir'], - cube_spacing_m = params['cube_spacing_in_m'], - ) - except RuntimeError: - logger.exception("Date %s failed", t) - continue - - # store in memory - wetDelays.append(wet_delay) - hydroDelays.append(hydro_delay) + # dont process the delays for download only + if dl_only: + continue + + if len(wfiles)>1: + import xarray + ds1 = xarray.open_dataset(wfiles[0]) + ds2 = xarray.open_dataset(wfiles[1]) + ds = ds1 + for var in ['wet', 'hydro', 'wet_total', 'hydro_total']: + ds[var] = (ds1[var] + ds2[var]) * 0.5 + ds.attrs['Date1'] = 0 + ds.attrs['Date2'] = 0 + weather_model_file = os.path.splitext(wfiles[0])[0] + '_interp.nc' + ds.to_netcdf(weather_model_file) + else: + weather_model_file = wfiles[0] - ########################################################### - # Write the delays to file - if params['interpolate_wmtimes']: - wet_delay, hydro_delay = avg_delays(wetDelays, hydroDelays, t, times[0], times[1]) + # Now process the delays + try: + wet_delay, hydro_delay = tropo_delay( + t, weather_model_file, aoi, los, + height_levels = params['height_levels'], + out_proj = params['output_projection'], + look_dir = params['look_dir'], + cube_spacing_m = params['cube_spacing_in_m'], + ) + except RuntimeError: + logger.exception("Date %s failed", t) + continue # Different options depending on the inputs if los.is_Projected(): From 4b69a2a0ad4f1e38b83284f420d6a57dfff33d32 Mon Sep 17 00:00:00 2001 From: Jeremy Maurer Date: Wed, 1 Mar 2023 08:44:34 -0600 Subject: [PATCH 12/19] update weights to be weighted by time nearness --- tools/RAiDER/cli/raider.py | 28 +++++++++++++++++++++++----- 1 file changed, 23 insertions(+), 5 deletions(-) diff --git a/tools/RAiDER/cli/raider.py b/tools/RAiDER/cli/raider.py index 9c97af33b..25184fa8f 100644 --- a/tools/RAiDER/cli/raider.py +++ b/tools/RAiDER/cli/raider.py @@ -1,9 +1,13 @@ import argparse +import datetime import os import shutil import sys import yaml +import numpy as np +import xarray as xr + from textwrap import dedent from RAiDER import aws @@ -241,15 +245,29 @@ def calcDelays(iargs=None): continue if len(wfiles)>1: - import xarray - ds1 = xarray.open_dataset(wfiles[0]) - ds2 = xarray.open_dataset(wfiles[1]) + ds1 = xr.open_dataset(wfiles[0]) + ds2 = xr.open_dataset(wfiles[1]) + + # calculate relative weights of each dataset + date1 = datetime.datetime.strptime(ds1.attrs['datetime'], '%Y_%m_%dT%H_%M_%S') + date2 = datetime.datetime.strptime(ds2.attrs['datetime'], '%Y_%m_%dT%H_%M_%S') + wgts = [np.abs((t - date1).seconds / (date2 - date1).seconds), np.abs((date2 - t).seconds / (date2 - date1).seconds)] + try: + assert np.sum(wgts)==1 + except AssertionError: + logging.error('Time interpolation weights do not sum to one, something is off with the dates') + raise ValueError('Time interpolation weights do not sum to one, something is off with the dates') + + # combine datasets ds = ds1 for var in ['wet', 'hydro', 'wet_total', 'hydro_total']: - ds[var] = (ds1[var] + ds2[var]) * 0.5 + ds[var] = (wgts[0] * ds1[var]) + (wgts[1] * ds2[var]) ds.attrs['Date1'] = 0 ds.attrs['Date2'] = 0 - weather_model_file = os.path.splitext(wfiles[0])[0] + '_interp.nc' + weather_model_file = os.path.join( + os.path.dirname(wfiles[0]), + os.path.basename(wfiles[0]).split('_')[0] + '_' + t.strftime('%Y_%m_%dT%H_%M_%S') + '_timeInterp_' + '_'.join(wfiles[0].split('_')[-4:]), + ) ds.to_netcdf(weather_model_file) else: weather_model_file = wfiles[0] From 49184ae0cfc2e97c3de4af96eb7b04d36317fd15 Mon Sep 17 00:00:00 2001 From: Jeremy Maurer Date: Tue, 7 Mar 2023 16:01:59 -0600 Subject: [PATCH 13/19] add utilfcn tests clean up and simplify time pulling --- test/test_temporal_interpolate.py | 4 +- test/test_util.py | 26 ++++++++- tools/RAiDER/aria/prepFromGUNW.py | 2 +- tools/RAiDER/cli/__init__.py | 2 +- tools/RAiDER/cli/raider.py | 11 ++-- tools/RAiDER/cli/raider.yaml | 2 +- tools/RAiDER/constants.py | 2 + tools/RAiDER/delay.py | 24 -------- tools/RAiDER/models/weatherModel.py | 3 + tools/RAiDER/utilFcns.py | 87 ++++++++++++++++++----------- 10 files changed, 95 insertions(+), 68 deletions(-) diff --git a/test/test_temporal_interpolate.py b/test/test_temporal_interpolate.py index 5b132b7ea..6bd93d4bf 100644 --- a/test/test_temporal_interpolate.py +++ b/test/test_temporal_interpolate.py @@ -90,7 +90,7 @@ def test_cube_timemean(): assert np.isclose(proc.returncode, 0) ## run interpolation in the middle of the two - grp['time_group'] = {'time': ti, 'interpolate_wmtimes': True} + grp['time_group'] = {'time': ti, 'interpolate_time': True} cfg = update_yaml(grp) cmd = f'raider.py {cfg}' @@ -150,7 +150,7 @@ def test_cube_weighting(): assert np.isclose(proc.returncode, 0) ## run interpolation very near the first - grp['time_group'] = {'time': ti, 'interpolate_wmtimes': True} + grp['time_group'] = {'time': ti, 'interpolate_time': True} cfg = update_yaml(grp) cmd = f'raider.py {cfg}' diff --git a/test/test_util.py b/test/test_util.py index 176c654eb..5539fbb70 100644 --- a/test/test_util.py +++ b/test/test_util.py @@ -13,7 +13,7 @@ _least_nonzero, cosd, rio_open, sind, writeArrayToRaster, rio_profile, rio_extents, getTimeFromFile, enu2ecef, ecef2enu, - transform_bbox, clip_bbox + transform_bbox, clip_bbox, get_nearest_wmtimes, ) @@ -462,3 +462,27 @@ def test_clip_bbox(): snwe_in = [34.005, 35.0006, -76.999, -76.0] assert clip_bbox(wesn, 0.01) == wesn assert clip_bbox(snwe_in, 0.01) == snwe + +def test_get_nearest_wmtimes(): + t0 = datetime.datetime(2020,1,1,11,35,0) + test_out = get_nearest_wmtimes(t0, 3) + true_out = [datetime.datetime(2020, 1, 1, 9, 0), datetime.datetime(2020, 1, 1, 12, 0)] + assert [t == t0 for t, t0 in zip(test_out, true_out)] + +def test_get_nearest_wmtimes_2(): + t0 = datetime.datetime(2020,1,1,11,3,0) + test_out = get_nearest_wmtimes(t0, 1) + true_out = [datetime.datetime(2020, 1, 1, 11, 0)] + assert [t == t0 for t, t0 in zip(test_out, true_out)] + +def test_get_nearest_wmtimes_3(): + t0 = datetime.datetime(2020,1,1,11,57,0) + test_out = get_nearest_wmtimes(t0, 3) + true_out = [datetime.datetime(2020, 1, 1, 12, 0)] + assert [t == t0 for t, t0 in zip(test_out, true_out)] + +def test_get_nearest_wmtimes_4(): + t0 = datetime.datetime(2020,1,1,11,25,0) + test_out = get_nearest_wmtimes(t0, 1) + true_out = [datetime.datetime(2020, 1, 1, 11, 0), datetime.datetime(2020, 1, 1, 12, 0)] + assert [t == t0 for t, t0 in zip(test_out, true_out)] diff --git a/tools/RAiDER/aria/prepFromGUNW.py b/tools/RAiDER/aria/prepFromGUNW.py index b0a465dd9..c701b3f6c 100644 --- a/tools/RAiDER/aria/prepFromGUNW.py +++ b/tools/RAiDER/aria/prepFromGUNW.py @@ -245,7 +245,7 @@ def main(args): 'aoi_group' : {'bounding_box': GUNWObj.SNWE}, 'height_group' : {'height_levels': GUNWObj.heights}, 'date_group': {'date_list': str(GUNWObj.dates)}, - 'time_group': {'time': GUNWObj.ref_time, 'interpolate_wmtimes': True}, + 'time_group': {'time': GUNWObj.ref_time, 'interpolate_time': True}, 'los_group' : {'ray_trace': True, 'orbit_file': GUNWObj.OrbitFiles, 'wavelength': GUNWObj.wavelength, diff --git a/tools/RAiDER/cli/__init__.py b/tools/RAiDER/cli/__init__.py index bbc0ebbc6..3b6156ef2 100644 --- a/tools/RAiDER/cli/__init__.py +++ b/tools/RAiDER/cli/__init__.py @@ -39,6 +39,6 @@ class AttributeDict(dict): output_directory='.', weather_model_directory=None, output_projection='EPSG:4326', - interpolate_wmtimes=False + interpolate_time=True, ) ) diff --git a/tools/RAiDER/cli/raider.py b/tools/RAiDER/cli/raider.py index 25184fa8f..cc4785a27 100644 --- a/tools/RAiDER/cli/raider.py +++ b/tools/RAiDER/cli/raider.py @@ -16,6 +16,7 @@ from RAiDER.cli.parser import add_out, add_cpus, add_verbose from RAiDER.cli.validators import DateListAction, date_type from RAiDER.models.allowed import ALLOWED_MODELS +from RAiDER.utilFcns import get_dt HELP_MESSAGE = """ @@ -125,7 +126,7 @@ def drop_nans(d): def calcDelays(iargs=None): """ Parse command line arguments using argparse. """ import RAiDER - from RAiDER.delay import tropo_delay, avg_delays + from RAiDER.delay import tropo_delay from RAiDER.checkArgs import checkArgs from RAiDER.processWM import prepareWeatherModel from RAiDER.utilFcns import writeDelays, get_nearest_wmtimes @@ -223,8 +224,10 @@ def calcDelays(iargs=None): logger.debug(f'Date: {t.strftime("%Y%m%d")}') logger.debug('Beginning weather model pre-processing') - - times = get_nearest_wmtimes(t, model._Name) if params['interpolate_wmtimes'] else [t] + # Grab the closest two times unless the user specifies 'nearest' + # If the model time_delta is not specified then use 6 + # The two datetimes will be combined to a single file and processed + times = get_nearest_wmtimes(t, [model.dtime() if model.dtime() is not None else 6][0]) if params['interpolate_time'] else [t] wfiles = [] for tt in times: try: @@ -251,7 +254,7 @@ def calcDelays(iargs=None): # calculate relative weights of each dataset date1 = datetime.datetime.strptime(ds1.attrs['datetime'], '%Y_%m_%dT%H_%M_%S') date2 = datetime.datetime.strptime(ds2.attrs['datetime'], '%Y_%m_%dT%H_%M_%S') - wgts = [np.abs((t - date1).seconds / (date2 - date1).seconds), np.abs((date2 - t).seconds / (date2 - date1).seconds)] + wgts = [ 1 - get_dt(t, date1) / get_dt(date2, date1), 1 - get_dt(date2, t) / get_dt(date2, date1)] try: assert np.sum(wgts)==1 except AssertionError: diff --git a/tools/RAiDER/cli/raider.yaml b/tools/RAiDER/cli/raider.yaml index bd854adaf..33fc5be05 100644 --- a/tools/RAiDER/cli/raider.yaml +++ b/tools/RAiDER/cli/raider.yaml @@ -58,7 +58,7 @@ date_group: time_group: time: end_time: - interpolate_wmtimes: false # Interpolate between weather model times if true, else use nearest time + interpolate_time: True # if True will interpolate two closest times, if False will take the nearest time ########## 4. Area of Interest diff --git a/tools/RAiDER/constants.py b/tools/RAiDER/constants.py index 8b9c84913..fe6fe1315 100644 --- a/tools/RAiDER/constants.py +++ b/tools/RAiDER/constants.py @@ -19,3 +19,5 @@ R_EARTH_MAX = 6378137 R_EARTH_MIN = 6356752 +_THRESHOLD_SECONDS = 5 * 60 # Threshold delta_time in seconds + diff --git a/tools/RAiDER/delay.py b/tools/RAiDER/delay.py index 79df68b4e..18d354ad5 100755 --- a/tools/RAiDER/delay.py +++ b/tools/RAiDER/delay.py @@ -412,27 +412,3 @@ def _build_cube_ray( if output_created_here: return outputArrs - - -def avg_delays(wetDelays, hydroDelays, dtref, dt1, dt2): - import xarray as xr - if hydroDelays[0] is None: - ## difference in seconds between the model dates and the reference date - dsecs = np.array([(dtref-dt1).seconds, (dt2-dtref).seconds]) - wgts = xr.DataArray(1/dsecs, dims=['time'], coords={'time': [1,2]}) - ds_out_wet = wetDelays[0].copy() # store the updated data in here - for kind in 'hydro wet'.split(): - da1 = wetDelays[0][kind].assign_coords(time=1) - da2 = wetDelays[1][kind].assign_coords(time=2) - da = xr.concat([da1, da2], 'time').weighted(wgts) - da_mu = da.mean('time') - ds_out_wet[kind] = da_mu - - ds_out_wet.attrs['source'] = [wetDelays[i].attrs['source'] for i in range(2)] - ds_out_hydro = None - - else: - for kind in 'hydro wet'.split(): - raise Exception('Not supported yet.') - - return ds_out_wet, ds_out_hydro \ No newline at end of file diff --git a/tools/RAiDER/models/weatherModel.py b/tools/RAiDER/models/weatherModel.py index 53abee761..4db2b473a 100755 --- a/tools/RAiDER/models/weatherModel.py +++ b/tools/RAiDER/models/weatherModel.py @@ -121,6 +121,9 @@ def __str__(self): def Model(self): return self._Name + + def dtime(self): + return self._time_res def fetch(self, out, ll_bounds, time): ''' diff --git a/tools/RAiDER/utilFcns.py b/tools/RAiDER/utilFcns.py index ff16600e7..bb46dbb7d 100755 --- a/tools/RAiDER/utilFcns.py +++ b/tools/RAiDER/utilFcns.py @@ -20,6 +20,7 @@ _g0 as g0, R_EARTH_MAX as Rmax, R_EARTH_MIN as Rmin, + _THRESHOLD_SECONDS, ) from RAiDER.logger import logger @@ -1137,39 +1138,57 @@ def transform_coords(proj1, proj2, x, y): return transformer.transform(x, y) -def get_nearest_wmtimes(dt0, wm): - """" Get the nearest model times""" - import pandas as pd - from dateutil.relativedelta import relativedelta - from datetime import time as dtime - ## hourly time step availability - not sure about ECMWF / HRRR - wm_hour = {'GMAO': 3, 'ERA5':1, 'ERA-5': 1, 'ERA5T': 1, 'HRES':6, 'HRRR': 1, 'NCMR': '?'} - wm_hours = np.arange(0, 24, wm_hour[wm]) - - dates = [] - times = [] - - ref_time = dt0 - ## sort the model and GUNW times - # date = dt0.date() - mod_times = [datetime.strptime(f'{str(dt0.date())}-{dt}', '%Y-%m-%d-%H') for dt in wm_hours] + [str(dt0)] - ser = pd.Series(mod_times, name='datetimes', dtype='datetime64[ns]') - ix_mod = ser.index[-1] # to get the nearest times - df = ser.sort_values().reset_index() - ix_mod_s = df[df['index']==ix_mod].index.item() - - ## case 1: index is not at the end - if ix_mod_s < df.index[-1]: - dt1, dt2 = df.loc[ix_mod_s-1, 'datetimes'], df.loc[ix_mod_s+1, 'datetimes'] - dt2 = dt2.to_pydatetime() - - ## case 2: close to midnight next day; add one to the date and start at midnight - else: - dt1 = df.loc[ix_mod_s-1, 'datetimes'] - - # handle leap years - days = 2 if (dt1.year % 4 == 0 and dt1.month == 2 and dt1.day == 28) else 1 - dt2 = datetime.combine(dt1 + relativedelta(days=days), dtime()) +def get_nearest_wmtimes(t0, time_delta): + """" + Get the nearest two available times to the requested time given a time step - return dt1.to_pydatetime(), dt2 + Args: + t0 - user-requested Python datetime + time_delta - time interval of weather model + + Returns: + tuple: list of datetimes representing the one or two closest + available times to the requested time + + Example: + >>> import datetime + >>> from RAiDER.utilFcns import get_nearest_wmtimes + >>> t0 = datetime.datetime(2020,1,1,11,35,0) + >>> get_nearest_wmtimes(t0, 3) + (datetime.datetime(2020, 1, 1, 9, 0), datetime.datetime(2020, 1, 1, 12, 0)) + """ + # get the closest time available + tclose = round_time(t0, roundTo = time_delta * 60 *60) + + # Just calculate both options and take the closest + t2_1 = tclose + timedelta(hours=time_delta) + t2_2 = tclose - timedelta(hours=time_delta) + t2 = [t2_1 if get_dt(t2_1, t0) < get_dt(t2_2, t0) else t2_2][0] + + # If you're within 5 minutes just take the closest time + if get_dt(tclose, t0) < _THRESHOLD_SECONDS: + return [tclose] + else: + if t2 > tclose: + return [tclose, t2] + else: + return [t2, tclose] +def get_dt(t1,t2): + ''' + Helper function for getting the absolute difference in seconds between + two python datetimes + + Args: + t1, t2 - Python datetimes + + Returns: + Absolute difference in seconds between the two inputs + + Examples: + >>> import datetime + >>> from RAiDER.utilFcns import get_dt + >>> get_dt(datetime.datetime(2020,1,1,5,0,0), datetime.datetime(2020,1,1,0,0,0)) + 18000.0 + ''' + return np.abs((t1 - t2).total_seconds()) From 2594e74094c1cd12483e952ee7b28612fef3fee6 Mon Sep 17 00:00:00 2001 From: Jeremy Maurer Date: Tue, 7 Mar 2023 20:58:53 -0600 Subject: [PATCH 14/19] update CHANGELOG --- CHANGELOG.md | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7bfe5cc2c..3b566c9f9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,8 @@ 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). +## [0.4.3] ++ Temporal interpolation of delays if the requested datetime is more than _THRESHOLD_SECONDS away from the closest weather model available time and `interpolate_time = True` (default behavior) + Add assert statement to raise error if the delay cube for each SAR date in a GUNW IFG is not written ## [0.4.2] @@ -19,11 +21,7 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). + Only orbit file is used (even if multiple specified) to minimize errors and ensure consistency over region + GUNW packaging is restructed to store SLC (ref and sec) wet and tropo delays rather than the differential -##[0.4.2] - -### New/Updated Features -+ `calcDelaysGUNW` allows processing with any supported weather model as listed in [`RAiDER.models.allowed.ALLOWED_MODELS`](https://github.com/dbekaert/RAiDER/blob/dev/tools/RAiDER/models/allowed.py). - +##[0.4.1] ### New/Updated Features + Reorder target points for intersection From 911be67773548f31e8a32c79e0242250e0712704 Mon Sep 17 00:00:00 2001 From: bbuzz31 Date: Wed, 8 Mar 2023 16:32:19 -0800 Subject: [PATCH 15/19] clean up util fns and wm time interval --- tools/RAiDER/constants.py | 9 +++--- tools/RAiDER/models/ecmwf.py | 4 +-- tools/RAiDER/models/gmao.py | 7 +++-- tools/RAiDER/models/hres.py | 4 +-- tools/RAiDER/models/hrrr.py | 6 ++-- tools/RAiDER/models/ncmr.py | 6 ++-- tools/RAiDER/models/weatherModel.py | 11 ++++++- tools/RAiDER/models/wrf.py | 4 ++- tools/RAiDER/utilFcns.py | 45 +++++++++++++++-------------- 9 files changed, 57 insertions(+), 39 deletions(-) diff --git a/tools/RAiDER/constants.py b/tools/RAiDER/constants.py index fe6fe1315..234fd9ef6 100644 --- a/tools/RAiDER/constants.py +++ b/tools/RAiDER/constants.py @@ -13,11 +13,12 @@ _CUBE_SPACING_IN_M = float(2000) # Horizontal spacing of cube -_g0 = np.float64(9.80665) +_g0 = np.float64(9.80665) # Standard gravitational constant +_g1 = np.float64(9.80616) # Gravitational constant @ 45° latitude used for corrections of earth's centrifugal force _RE = np.float64(6371008.7714) -R_EARTH_MAX = 6378137 -R_EARTH_MIN = 6356752 +R_EARTH_MAX = 6378137 # WGS84 +R_EARTH_MIN = 6356752 # WGS84 -_THRESHOLD_SECONDS = 5 * 60 # Threshold delta_time in seconds +_THRESHOLD_SECONDS = 1 * 60 # Threshold delta_time in seconds diff --git a/tools/RAiDER/models/ecmwf.py b/tools/RAiDER/models/ecmwf.py index fde5fef3b..f65af82ec 100755 --- a/tools/RAiDER/models/ecmwf.py +++ b/tools/RAiDER/models/ecmwf.py @@ -14,8 +14,8 @@ A_137_HRES, B_137_HRES, ) -from RAiDER.models.weatherModel import WeatherModel +from RAiDER.models.weatherModel import WeatherModel, TIME_RES class ECMWF(WeatherModel): ''' @@ -31,7 +31,7 @@ def __init__(self): self._k2 = 0.233 # [K/Pa] self._k3 = 3.75e3 # [K^2/Pa] - self._time_res = 1 + self._time_res = TIME_RES['ECMWF'] self._lon_res = 0.2 self._lat_res = 0.2 diff --git a/tools/RAiDER/models/gmao.py b/tools/RAiDER/models/gmao.py index ec7e21245..e0c9d3fd0 100755 --- a/tools/RAiDER/models/gmao.py +++ b/tools/RAiDER/models/gmao.py @@ -7,7 +7,7 @@ import pydap.client from pyproj import CRS -from RAiDER.models.weatherModel import WeatherModel +from RAiDER.models.weatherModel import WeatherModel, TIME_RES from RAiDER.logger import logger from RAiDER.utilFcns import writeWeatherVars2NETCDF4, round_time, requests_retry_session from RAiDER.models.model_levels import ( @@ -37,7 +37,8 @@ def __init__(self): self._k2 = 0.233 # [K/Pa] self._k3 = 3.75e3 # [K^2/Pa] - self._time_res = 3 + self._time_res = TIME_RES[self._dataset.upper()] + # horizontal grid spacing self._lat_res = 0.25 @@ -69,7 +70,7 @@ def _fetch(self, out): T0 = dt.datetime(2017, 12, 1, 0, 0, 0) # round time to nearest third hour time1 = time - time = round_time(time, 3 * 60 * 60) + time = round_time(time, self._time_res * 60 * 60) if not time1 == time: logger.warning('Rounded given hour from %d to %d. May be incorrect near beginning/end of day.', time1.hour, time.hour) diff --git a/tools/RAiDER/models/hres.py b/tools/RAiDER/models/hres.py index 4dfc45dd8..97be42c87 100755 --- a/tools/RAiDER/models/hres.py +++ b/tools/RAiDER/models/hres.py @@ -5,7 +5,7 @@ from pyproj import CRS from RAiDER.models.ecmwf import ECMWF -from RAiDER.models.weatherModel import WeatherModel +from RAiDER.models.weatherModel import WeatherModel, TIME_RES from RAiDER.models.model_levels import ( LEVELS_91_HEIGHTS, LEVELS_25_HEIGHTS, @@ -28,7 +28,6 @@ def __init__(self, level_type='ml'): self._k2 = 0.233 # [K/Pa] self._k3 = 3.75e3 # [K^2/Pa] - self._time_res = 6 # 9 km horizontal grid spacing. This is only used for extending the download-buffer, i.e. not in subsequent processing. self._lon_res = 9. / 111 # 0.08108115 @@ -44,6 +43,7 @@ def __init__(self, level_type='ml'): self._Name = 'HRES' self._proj = CRS.from_epsg(4326) + self._time_res = TIME_RES[self._dataset.upper()] # Tuple of min/max years where data is available. self._valid_range = (datetime.datetime(1983, 4, 20), "Present") # Availability lag time in days diff --git a/tools/RAiDER/models/hrrr.py b/tools/RAiDER/models/hrrr.py index 88599bf93..41d38f716 100644 --- a/tools/RAiDER/models/hrrr.py +++ b/tools/RAiDER/models/hrrr.py @@ -13,7 +13,9 @@ from RAiDER.logger import logger from RAiDER.utilFcns import rio_profile, rio_extents -from RAiDER.models.weatherModel import WeatherModel, transform_coords +from RAiDER.models.weatherModel import ( + WeatherModel, transform_coords, TIME_RES +) from RAiDER.models.model_levels import ( LEVELS_137_HEIGHTS, ) @@ -30,7 +32,7 @@ def __init__(self): self._classname = 'hrrr' self._dataset = 'hrrr' - self._time_res = 1 + self._time_res = TIME_RES[self._dataset.upper()] # Tuple of min/max years where data is available. self._valid_range = (datetime.datetime(2016, 7, 15), "Present") diff --git a/tools/RAiDER/models/ncmr.py b/tools/RAiDER/models/ncmr.py index 84fa62d35..86e3b418a 100755 --- a/tools/RAiDER/models/ncmr.py +++ b/tools/RAiDER/models/ncmr.py @@ -10,7 +10,7 @@ from pyproj import CRS -from RAiDER.models.weatherModel import WeatherModel +from RAiDER.models.weatherModel import WeatherModel, TIME_RES from RAiDER.logger import logger from RAiDER.utilFcns import ( writeWeatherVars2NETCDF4, @@ -36,12 +36,12 @@ def __init__(self): self._classname = 'ncmr' # name of the custom weather model self._dataset = 'ncmr' # same name as above self._Name = 'NCMR' # name of the new weather model (in Capital) + self._time_res = TIME_RES[self._dataset.upper()] # Tuple of min/max years where data is available. self._valid_range = (datetime.datetime(2015, 12, 1), "Present") # Availability lag time in days/hours self._lag_time = datetime.timedelta(hours=6) - self._time_res = 1 # model constants self._k1 = 0.776 # [K/Pa] @@ -67,7 +67,7 @@ def _fetch(self, out): Fetch weather model data from NCMR: note we only extract the lat/lon bounds for this weather model; fetching data is not needed here as we don't actually download data , data exist in same system ''' - time = self._time + time = self._time # Auxillary function: ''' diff --git a/tools/RAiDER/models/weatherModel.py b/tools/RAiDER/models/weatherModel.py index 4db2b473a..6238d9a28 100755 --- a/tools/RAiDER/models/weatherModel.py +++ b/tools/RAiDER/models/weatherModel.py @@ -22,6 +22,15 @@ robmax, robmin, write2NETCDF4core, calcgeoh, transform_coords ) +TIME_RES = {'GMAO': 3, + 'ECMWF': 1, + 'HRES': 6, + 'HRRR': 1, + 'WRF': 1, + 'NCMR': 1 + } + + class WeatherModel(ABC): ''' @@ -121,7 +130,7 @@ def __str__(self): def Model(self): return self._Name - + def dtime(self): return self._time_res diff --git a/tools/RAiDER/models/wrf.py b/tools/RAiDER/models/wrf.py index d35f795bd..827495a54 100644 --- a/tools/RAiDER/models/wrf.py +++ b/tools/RAiDER/models/wrf.py @@ -2,7 +2,7 @@ import scipy.io.netcdf as netcdf from pyproj import CRS, Transformer -from RAiDER.models.weatherModel import WeatherModel +from RAiDER.models.weatherModel import WeatherModel, TIME_RES # Need to incorporate this snippet into this part of the code. @@ -27,9 +27,11 @@ def __init__(self): self._k2 = 0.233 # K/Pa self._k3 = 3.75e3 # K^2/Pa + # Currently WRF is using RH instead of Q to get E self._humidityType = 'rh' self._Name = 'WRF' + self._time_res = TIME_RES[self._Name] def _fetch(self): pass diff --git a/tools/RAiDER/utilFcns.py b/tools/RAiDER/utilFcns.py index bb46dbb7d..ba62fe31f 100755 --- a/tools/RAiDER/utilFcns.py +++ b/tools/RAiDER/utilFcns.py @@ -18,6 +18,7 @@ from RAiDER.constants import ( _g0 as g0, + _g1 as G1, R_EARTH_MAX as Rmax, R_EARTH_MIN as Rmin, _THRESHOLD_SECONDS, @@ -397,30 +398,32 @@ def _get_g_ll(lats): ''' Compute the variation in gravity constant with latitude ''' - # TODO: verify these constants. In particular why is the reference g different from self._g0? - return 9.80616 * (1 - 0.002637 * cosd(2 * lats) + 0.0000059 * (cosd(2 * lats))**2) + return G1 * (1 - 0.002637 * cosd(2 * lats) + 0.0000059 * (cosd(2 * lats))**2) def _get_Re(lats): ''' Returns: the ellipsoid as a fcn of latitude + More precisely: Calculates the GEOID radius at a given latitude on WGS84 ellipsoid ''' - # TODO: verify constants, add to base class constants? return np.sqrt(1 / (((cosd(lats)**2) / Rmax**2) + ((sind(lats)**2) / Rmin**2))) def _geo_to_ht(lats, hts): - """Convert geopotential height to altitude.""" - # Convert geopotential to geometric height. This comes straight from - # TRAIN - # Map of g with latitude (I'm skeptical of this equation - Ray) + """ Returns geometric height above the reference ellipsoid + + Note that this formula technically computes height above geoid + but the geoid is actually a perfect sphere; + Thus returned heights are above a reference ellipsoid == sphere + + https://confluence.ecmwf.int/display/CKB/ERA5%3A+compute+pressure+and+geopotential+on+model+levels%2C+geopotential+height+and+geometric+height + + """ g_ll = _get_g_ll(lats) Re = _get_Re(lats) # Calculate Geometric Height, h h = (hts * Re) / (g_ll / g0 * Re - hts) - # from metpy - # return (geopotential * Re) / (g0 * Re - geopotential) return h @@ -1143,20 +1146,20 @@ def get_nearest_wmtimes(t0, time_delta): Get the nearest two available times to the requested time given a time step Args: - t0 - user-requested Python datetime + t0 - user-requested Python datetime time_delta - time interval of weather model - + Returns: - tuple: list of datetimes representing the one or two closest + tuple: list of datetimes representing the one or two closest available times to the requested time - - Example: + + Example: >>> import datetime >>> from RAiDER.utilFcns import get_nearest_wmtimes >>> t0 = datetime.datetime(2020,1,1,11,35,0) >>> get_nearest_wmtimes(t0, 3) (datetime.datetime(2020, 1, 1, 9, 0), datetime.datetime(2020, 1, 1, 12, 0)) - """ + """ # get the closest time available tclose = round_time(t0, roundTo = time_delta * 60 *60) @@ -1176,16 +1179,16 @@ def get_nearest_wmtimes(t0, time_delta): def get_dt(t1,t2): ''' - Helper function for getting the absolute difference in seconds between + Helper function for getting the absolute difference in seconds between two python datetimes - Args: + Args: t1, t2 - Python datetimes - - Returns: + + Returns: Absolute difference in seconds between the two inputs - - Examples: + + Examples: >>> import datetime >>> from RAiDER.utilFcns import get_dt >>> get_dt(datetime.datetime(2020,1,1,5,0,0), datetime.datetime(2020,1,1,0,0,0)) From d64d733c329451f8cb6ebce8548d5fb8270efcc7 Mon Sep 17 00:00:00 2001 From: bbuzz31 Date: Wed, 8 Mar 2023 16:36:50 -0800 Subject: [PATCH 16/19] update changelog --- CHANGELOG.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3b566c9f9..db82e7073 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,8 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [0.4.3] + Temporal interpolation of delays if the requested datetime is more than _THRESHOLD_SECONDS away from the closest weather model available time and `interpolate_time = True` (default behavior) + Add assert statement to raise error if the delay cube for each SAR date in a GUNW IFG is not written ++ Verify some constants / equations and remove the comments questioning them ++ Relocate the time resolution of wmodels to one spot ## [0.4.2] @@ -21,7 +23,7 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). + Only orbit file is used (even if multiple specified) to minimize errors and ensure consistency over region + GUNW packaging is restructed to store SLC (ref and sec) wet and tropo delays rather than the differential -##[0.4.1] +## [0.4.1] ### New/Updated Features + Reorder target points for intersection From 092d2941dcc0134ca61588629c0381354c534657 Mon Sep 17 00:00:00 2001 From: bbuzz31 Date: Wed, 8 Mar 2023 17:53:01 -0800 Subject: [PATCH 17/19] few bug fixes --- tools/RAiDER/aria/prepFromGUNW.py | 4 ++-- tools/RAiDER/models/gmao.py | 2 +- tools/RAiDER/models/hrrr.py | 10 ++++++++-- 3 files changed, 11 insertions(+), 5 deletions(-) diff --git a/tools/RAiDER/aria/prepFromGUNW.py b/tools/RAiDER/aria/prepFromGUNW.py index c701b3f6c..ad4d8b3d7 100644 --- a/tools/RAiDER/aria/prepFromGUNW.py +++ b/tools/RAiDER/aria/prepFromGUNW.py @@ -244,8 +244,8 @@ def main(args): 'cube_spacing_in_m': GUNWObj.spacing_m, 'aoi_group' : {'bounding_box': GUNWObj.SNWE}, 'height_group' : {'height_levels': GUNWObj.heights}, - 'date_group': {'date_list': str(GUNWObj.dates)}, - 'time_group': {'time': GUNWObj.ref_time, 'interpolate_time': True}, + 'date_group': {'date_list': GUNWObj.dates}, + 'time_group': {'time': GUNWObj.mid_time, 'interpolate_time': True}, 'los_group' : {'ray_trace': True, 'orbit_file': GUNWObj.OrbitFiles, 'wavelength': GUNWObj.wavelength, diff --git a/tools/RAiDER/models/gmao.py b/tools/RAiDER/models/gmao.py index e0c9d3fd0..f3aa5fdeb 100755 --- a/tools/RAiDER/models/gmao.py +++ b/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_time, requests_retry_session +from RAiDER.utilFcns import writeWeatherVars2NETCDF4, round_date, requests_retry_session from RAiDER.models.model_levels import ( LEVELS_137_HEIGHTS, ) diff --git a/tools/RAiDER/models/hrrr.py b/tools/RAiDER/models/hrrr.py index 41d38f716..5bcb7497e 100644 --- a/tools/RAiDER/models/hrrr.py +++ b/tools/RAiDER/models/hrrr.py @@ -12,7 +12,7 @@ from pyproj import CRS, Transformer from RAiDER.logger import logger -from RAiDER.utilFcns import rio_profile, rio_extents +from RAiDER.utilFcns import rio_profile, rio_extents, round_date from RAiDER.models.weatherModel import ( WeatherModel, transform_coords, TIME_RES ) @@ -264,8 +264,14 @@ def _download_hrrr_file(self, DATE, out, model='hrrr', product='prs', fxx=0, ver ''' Download a HRRR model ''' + ## TODO: Check how Herbie does temporal interpolation + # corrected_DT = round_date(DATE, datetime.timedelta(hours=self._time_res)) + # if not corrected_DT == DATE: + # logger.warning('Rounded given datetime from %s to %s', DATE, corrected_DT) + corrected_DT = DATE + H = Herbie( - DATE.strftime('%Y-%m-%d %H:%M'), + corrected_DT.strftime('%Y-%m-%d %H:%M'), model=model, product=product, overwrite=False, From 3d4d48be1af08cb78f61941688fd7ba0b78299d4 Mon Sep 17 00:00:00 2001 From: Jeremy Maurer Date: Thu, 9 Mar 2023 09:43:51 -0600 Subject: [PATCH 18/19] fix rebase bug and add a bunch of comments on the ellipsoidal height calculation --- tools/RAiDER/constants.py | 7 ++-- tools/RAiDER/models/ecmwf.py | 4 +- tools/RAiDER/models/weatherModel.py | 4 +- tools/RAiDER/utilFcns.py | 61 ++++++++++++++++++++++------- 4 files changed, 54 insertions(+), 22 deletions(-) diff --git a/tools/RAiDER/constants.py b/tools/RAiDER/constants.py index fe6fe1315..3b3343dc2 100644 --- a/tools/RAiDER/constants.py +++ b/tools/RAiDER/constants.py @@ -11,13 +11,12 @@ _ZREF = np.float64(15000) # maximum requierd height _STEP = np.float64(15.0) # integration step size in meters -_CUBE_SPACING_IN_M = float(2000) # Horizontal spacing of cube - _g0 = np.float64(9.80665) _RE = np.float64(6371008.7714) -R_EARTH_MAX = 6378137 -R_EARTH_MIN = 6356752 +R_EARTH_MAX_WGS84 = 6378137 +R_EARTH_MIN_WGS84 = 6356752 +_CUBE_SPACING_IN_M = float(2000) # Horizontal spacing of cube _THRESHOLD_SECONDS = 5 * 60 # Threshold delta_time in seconds diff --git a/tools/RAiDER/models/ecmwf.py b/tools/RAiDER/models/ecmwf.py index fde5fef3b..be71c827c 100755 --- a/tools/RAiDER/models/ecmwf.py +++ b/tools/RAiDER/models/ecmwf.py @@ -165,7 +165,7 @@ def _get_from_ecmwf(self, lat_min, lat_max, lat_step, lon_min, lon_max, corrected_date = util.round_date(time, datetime.timedelta(hours=6)) if not corrected_date == time: - logger.warning('Rounded given hour from %d to %d', time.hour, corrected_date.hour) + logger.warning('Rounded given datetime from %s to %s', time, corrected_date) server.retrieve({ "class": self._classname, # ERA-Interim @@ -222,7 +222,7 @@ def _get_from_cds( # round to the closest legal time corrected_date = util.round_date(acqTime, datetime.timedelta(hours=self._time_res)) if not corrected_date == acqTime: - logger.warning('Rounded given hour from %d to %d', acqTime.hour, corrected_date.hour) + logger.warning('Rounded given datetime from %s to %s', acqTime, corrected_date) # I referenced https://confluence.ecmwf.int/display/CKB/How+to+download+ERA5 dataDict = { diff --git a/tools/RAiDER/models/weatherModel.py b/tools/RAiDER/models/weatherModel.py index 4db2b473a..16b026b63 100755 --- a/tools/RAiDER/models/weatherModel.py +++ b/tools/RAiDER/models/weatherModel.py @@ -275,10 +275,10 @@ def _convertmb2Pa(self, pres): def _get_heights(self, lats, geo_hgt, geo_ht_fill=np.nan): ''' - Transform geo heights to actual heights + Transform geo heights to WGS84 ellipsoidal heights ''' geo_ht_fix = np.where(geo_hgt != geo_ht_fill, geo_hgt, np.nan) - self._zs = util._geo_to_ht(lats, geo_ht_fix) + self._zs = util.geo_to_ht(lats, geo_ht_fix) def _find_e(self): """Check the type of e-calculation needed""" diff --git a/tools/RAiDER/utilFcns.py b/tools/RAiDER/utilFcns.py index bb46dbb7d..c52f87b7b 100755 --- a/tools/RAiDER/utilFcns.py +++ b/tools/RAiDER/utilFcns.py @@ -18,8 +18,8 @@ from RAiDER.constants import ( _g0 as g0, - R_EARTH_MAX as Rmax, - R_EARTH_MIN as Rmin, + R_EARTH_MAX_WGS84 as Rmax, + R_EARTH_MIN_WGS84 as Rmin, _THRESHOLD_SECONDS, ) from RAiDER.logger import logger @@ -401,26 +401,59 @@ def _get_g_ll(lats): return 9.80616 * (1 - 0.002637 * cosd(2 * lats) + 0.0000059 * (cosd(2 * lats))**2) -def _get_Re(lats): +def get_Re(lats): ''' - Returns: the ellipsoid as a fcn of latitude + Returns: the ellipsoid as a fcn of latitude. + + The equation is a standard expression for g, see e.g., + https://arggit.usask.ca/ARGPackages/eratools/-/blob/b41e660c5392ede6b5ba93ffd1d61defe3ddcbc4/eratools/geopotential.py + + Args: + lats - ndarray of geodetic latitudes in degrees + + Returns: + ndarray of earth radius at each latitude + + Example: + >>> import numpy as np + >>> from RAiDER.utilFcns import get_Re + >>> output = get_Re(np.array([0, 30, 45, 60, 90])) + >>> output + array([6378137., 6372770.5219805, 6367417.56705189, 6362078.07851428, 6356752.]) + >>> assert output[0] == 6378137 # (Rmax) + >>> assert output[-1] == 6356752 # (Rmin) ''' - # TODO: verify constants, add to base class constants? return np.sqrt(1 / (((cosd(lats)**2) / Rmax**2) + ((sind(lats)**2) / Rmin**2))) -def _geo_to_ht(lats, hts): - """Convert geopotential height to altitude.""" - # Convert geopotential to geometric height. This comes straight from - # TRAIN - # Map of g with latitude (I'm skeptical of this equation - Ray) - g_ll = _get_g_ll(lats) - Re = _get_Re(lats) +def geo_to_ht(lats, hts): + """ + Convert geopotential height to ellipsoidal heights referenced to WGS84. + + Note that this formula technically computes height above geoid (geometric height) + but the geoid is actually a perfect sphere; + Thus returned heights are above a reference ellipsoid, which most assume to be + a sphere (e.g., ECMWF - see https://confluence.ecmwf.int/display/CKB/ERA5%3A+compute+pressure+and+geopotential+on+model+levels%2C+geopotential+height+and+geometric+height#ERA5:computepressureandgeopotentialonmodellevels,geopotentialheightandgeometricheight-Geopotentialheight + - "Geometric Height" and also https://confluence.ecmwf.int/display/CKB/ERA5%3A+data+documentation#ERA5:datadocumentation-Earthmodel). + However, by calculating the ellipsoid here we can directly reference to WGS84. + + Compare to MetPy: + (https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.geopotential_to_height.html) + # h = (geopotential * Re) / (g0 * Re - geopotential) + # Assumes a sphere instead of an ellipsoid + + Args: + lats - latitude of points of interest + hts - geopotential height at points of interest + + Returns: + ndarray: geometric heights. These are approximate ellipsoidal heights referenced to WGS84 + """ + g_ll = _get_g_ll(lats) # gravity function of latitude + Re = get_Re(lats) # Earth radius function of latitude # Calculate Geometric Height, h h = (hts * Re) / (g_ll / g0 * Re - hts) - # from metpy - # return (geopotential * Re) / (g0 * Re - geopotential) return h From 3df371eaee4c27470b517ecdcbb56d80b98ce749 Mon Sep 17 00:00:00 2001 From: Jeremy Maurer Date: Thu, 9 Mar 2023 15:33:59 -0600 Subject: [PATCH 19/19] fix import bug --- tools/RAiDER/models/gmao.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/RAiDER/models/gmao.py b/tools/RAiDER/models/gmao.py index f3aa5fdeb..e0c9d3fd0 100755 --- a/tools/RAiDER/models/gmao.py +++ b/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 writeWeatherVars2NETCDF4, round_time, requests_retry_session from RAiDER.models.model_levels import ( LEVELS_137_HEIGHTS, )