Merge pull request #499 from jlmaurer/new_time_interp

Temporal interpolation of delays

CHANGELOG.md

@@ -6,7 +6,11 @@ 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 
++ Verify some constants / equations and remove the comments questioning them
++ Relocate the time resolution of wmodels to one spot
 
 ## [0.4.2]
 

test/test_temporal_interpolate.py

@@ -0,0 +1,193 @@
+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
+
+
+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 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
+    S, N, W, E = 34, 35, -117, -116
+    date       = 20200130
+    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},
+            'weather_model': WM,
+            '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': ti, 'interpolate_time': True}
+    cfg  = update_yaml(grp)
+
+    cmd  = f'raider.py {cfg}'
+    proc = subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True)
+    assert np.isclose(proc.returncode, 0)
+
+
+    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{ti.replace(":", "")}_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():
+    """ 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)
+    ## make the lat lon grid
+    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, ti = dct_hrs[WM]
+
+    grp = {
+            'date_group': {'date_start': date},
+            'weather_model': WM,
+            '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': ti, 'interpolate_time': True}
+    cfg  = update_yaml(grp)
+
+    cmd  = f'raider.py {cfg}'
+    proc = subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True)
+
+    ## double check on weighting
+
+    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{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}{ti}', '%Y%m%d%H:%M:%S')
+
+    wgts  = np.array([(dt_ref-dt1).seconds, (dt2-dt_ref).seconds])
+    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])
+
+    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
+

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)]

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},
+           '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,

tools/RAiDER/cli/__init__.py

@@ -39,5 +39,6 @@ class AttributeDict(dict):
             output_directory='.',
             weather_model_directory=None,
             output_projection='EPSG:4326',
+            interpolate_time=True,
         )
     )

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
@@ -12,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 = """
@@ -124,7 +129,7 @@ def calcDelays(iargs=None):
     from RAiDER.delay import tropo_delay
     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'
@@ -218,21 +223,58 @@ 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,
-                ll_bounds=ll_bounds, # SNWE
-                wmLoc=params['weather_model_directory'],
-                makePlots=params['verbose'],
-            )
-        except RuntimeError:
-            logger.exception("Date %s failed", t)
-            continue
+
+        # 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:
+                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
 
+        if len(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  = [ 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:
+                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] = (wgts[0] * ds1[var]) + (wgts[1] * ds2[var])
+            ds.attrs['Date1'] = 0
+            ds.attrs['Date2'] = 0
+            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]
+
         # Now process the delays
         try:
             wet_delay, hydro_delay = tropo_delay(
@@ -246,10 +288,7 @@ def calcDelays(iargs=None):
             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")

tools/RAiDER/cli/raider.yaml

@@ -58,6 +58,7 @@ date_group:
 time_group:
   time:
   end_time:
+  interpolate_time: True # if True will interpolate two closest times, if False will take the nearest time
 
 
 ########## 4. Area of Interest

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
@@ -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: