Azimuth time interpolation for Weather Models in GUNW Workflow #571

cmarshak · 2023-07-25T00:08:17Z

Ports the azimuth time grid interpolation demonstrated here to RAiDER's GUNW workflow. Formerly the center time was used to interpolate across the entire weather cube for two weather model times. Now, it is based on the time in which the sensor passes the pixel at zero-doppler (using ISCE3 geo2rdr).

The main test will be to show that the center_time interpolation and the azimuth_time_grid interpolation doe not deviate more than 1 mm.

There is a change in the API for yaml (see here) and CLI (see here). Specifically, updates interpolate_time options to:

'none' (formerly False),
'center_time' (formerly True and no longer default), and
azimuth_time_grid (not implemented previously and now default only working for HRRR model)

This will make the default behavior for the GUNW workflow azimuth_time_grid, which only works for HRRR.

cmarshak · 2023-07-25T21:22:27Z

@dbekaert

@sssangha and I worked on this repo together: https://github.com/ACCESS-Cloud-Based-InSAR/s1_azimuth_time_grid
I have integrated into RAiDER for HRRR in the PR. It performs an inverse weighting based on a per pixel difference between the azimuth acquisition time (zero-doppler) and the model time of the weather model. The difference between center-time interpolation and this per pixel interpolation for the troposphere phase are on the order of 1 mm (not sure if this is too large given the model times are much greater than the azimuth timing I am putting in). The test to verify this for all relevant layers for a single GUNW is illustrated here:

RAiDER/test/test_GUNW.py

Lines 121 to 150 in 4e00730

    
           @pytest.mark.parametrize('model', ['HRRR']) 
        
           def test_azimuth_timing_against_interpolation(model, tmp_path, gunw_azimuth_test): 
        
               """This test shows that the azimuth timing interpolation does not deviate from 
        
               the center time by more than 1 mm for the HRRR model. This is expected since the model times are 
        
               6 hours apart and a the azimuth time is changing the interpolation weights for a given pixel at the order 
        
               of seconds and thus these two approaches are quite similar.""" 
        
               out_0 = gunw_azimuth_test.name.replace('.nc', '__ct_interp.nc') 
        
               out_1 = gunw_azimuth_test.name.replace('.nc', '__az_interp.nc') 
        
               out_path_0 = shutil.copy(gunw_azimuth_test, tmp_path / out_0) 
        
               out_path_1 = shutil.copy(gunw_azimuth_test, tmp_path / out_1) 
        
               cmd = f'raider.py ++process calcDelaysGUNW -f {out_path_0} -m HRRR -interp center_time' 
        
               subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True) 
        
               cmd = f'raider.py ++process calcDelaysGUNW -f {out_path_1} -m HRRR -interp azimuth_time_grid' 
        
               subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True) 
        
               for ifg_type in ['reference', 'secondary']: 
        
                   for var in ['troposphereHydrostatic', 'troposphereWet']: 
        
                       group = f'science/grids/corrections/external/troposphere/{model}/{ifg_type}' 
        
                       with xr.open_dataset(out_0, group=group) as ds: 
        
                           da_0 = ds[var] 
        
                       with xr.open_dataset(out_1, group=group) as ds: 
        
                           da_1 = ds[var] 
        
                       # diff * wavelength / (4 pi) transforms to meters; then x 1000 to mm 
        
                       diff_mm = (da_1 - da_0).data * 0.055465761572122574 / (4 * np.pi) * 1_000 
        
                       # Differences in mm are bounded by 1 
        
                       assert np.all(diff_mm < 1)

Let us know your thoughts.

Note: I updated the test above so all the data is mocked. The HRRR API (AWS s3) is not required.

cmarshak · 2023-07-28T00:32:14Z

@sssangha - the easiest way to use this PR is:

raider.py ++process calcDelaysGUNW -f <gunw_path> -m HRRR -interp azimuth_time_grid

cmarshak · 2023-08-14T21:14:48Z

Here is an example GUNW product that is also used in the test.

The command to generate is:

raider.py ++process calcDelaysGUNW -f S1-GUNW-A-R-064-tops-20210723_20210711-015000-00119W_00033N-PP-6267-v2_0_6.nc --weather-model HRRR -interp azimuth_time_grid

Below are some screenshots of the reference delay in radians for the above product.

bbuzz31 · 2023-08-15T16:08:44Z

Here is an example GUNW product that is also used in the test.

The command to generate is:
raider.py ++process calcDelaysGUNW -f S1-GUNW-A-R-064-tops-20210723_20210711-015000-00119W_00033N-PP-6267-v2_0_6.nc --weather-model HRRR -interp azimuth_time_grid
Below are some screenshots of the reference delay in radians for the above product.

Could you please show the residual between this and the original (uncorrected for azimuth time)?

cmarshak · 2023-08-15T16:22:09Z

I can do this in the afternoon. Please note the test suite for GUNWs include a test that shows the residual of said product above is below 1 mm. Do you want to see spatial pattern of said residual?

cmarshak · 2023-08-15T20:53:39Z

When I run the command line argument for GUNWs, note I get the following output in the log:

Processing slice 1 / 20: -500.0
Processing slice 2 / 20: 0.0
Processing slice 3 / 20: 500.0
Processing slice 4 / 20: 1000.0
Processing slice 5 / 20: 1500.0
Processing slice 6 / 20: 2000.0
Processing slice 7 / 20: 2500.0
Processing slice 8 / 20: 3000.0
Processing slice 9 / 20: 3500.0
Processing slice 10 / 20: 4000.0
Processing slice 11 / 20: 4500.0
Processing slice 12 / 20: 5000.0
Processing slice 13 / 20: 5500.0
Processing slice 14 / 20: 6000.0
Processing slice 15 / 20: 6500.0
Processing slice 16 / 20: 7000.0
Processing slice 17 / 20: 7500.0
Processing slice 18 / 20: 8000.0
Processing slice 19 / 20: 8500.0
Processing slice 20 / 20: 9000.0
CRITICAL: There are missing delay values. Check your inputs.

I also have weird nodata values in the HRRR.

Will open an issue ticket to illustrate.

cmarshak · 2023-08-16T16:58:25Z

From @sssangha:

Here is a series of quick figures. I’ll also paste the input GUNWs + sequence of few commands I used to extract these layers. I will go ahead and pass the latter into a formal notebook later.

Troposphere hydrostatic layer at 500 m elevation directly from the overlapping GUNWs, for the reference date:

Troposphere wet layer at 500 m elevation directly from the overlapping GUNWs, for the reference date:

Reference date extracted and stitched tropo total via ARIA-tools:

Secondary date extracted and stitched tropo total via ARIA-tools:

Everything overall looks seamless, but upon closer inspection in the GUNWs, I found surprisingly large mean/min values.

E.g. for the overlap region of the troposphere hydro layer:
STATISTICS_MAXIMUM=2.3955078125
STATISTICS_MEAN=-0.000690513815912
STATISTICS_MINIMUM=-3.829345703125
STATISTICS_STDDEV=0.067453400601181

For the overlap region of the troposphere wet layer:
STATISTICS_MAXIMUM=2.4646682739258
STATISTICS_MEAN=5.98180550712e-06
STATISTICS_MINIMUM=-1.6872520446777
STATISTICS_STDDEV=0.03005129298574

When I looked at the actual residual fields for the hydro and wet layers, respectively, I found a few strange, very narrow rows with some non-zero residual:

Here are the commands I used:

GUNWs used:

S1-GUNW-A-R-064-tops-20210723_20210711-015001-35393N_33512N-PP-6267-v2_0_4.nc
S1-GUNW-A-R-064-tops-20210723_20210711-014936-33900N_32016N-PP-53c5-v2_0_4.nc

# using raider to populate the GUNWs
raider.py ++process calcDelaysGUNW -f products/S1-GUNW-A-R-064-tops-20210723_20210711-015001-35393N_33512N-PP-6267-v2_0_4.nc -m HRRR -interp azimuth_time_grid
 
raider.py ++process calcDelaysGUNW -f products/S1-GUNW-A-R-064-tops-20210723_20210711-014936-33900N_32016N-PP-53c5-v2_0_4.nc -m HRRR -interp azimuth_time_grid  
 
# using ARIA-tools to extract and stitch the tropo layers
ariaExtract.py -f 'products/S1*' -l 'unwrappedPhase,troposphereTotal' -tm HRRR -d Download -w output_dir -of ISCE

bbuzz31 · 2023-08-17T22:18:32Z

@sssangha do you have a function handy for cropping to the overlap?

bbuzz31 · 2023-08-18T23:29:16Z

I am not seeing this residual in the delays within the GUNW from RAiDER (prior to running ARIA-tools and the intersection). I'll see if I can recreate after running AT and debug, but as the problem does not appear to be in RAiDER, I think this PR can be merged. Note, I also tested with GMAO (no azimuth interpolation) and there is no weird striping difference.
What do you think @dbekaert?

Before #575

After #575

bbuzz31

LGTM nice work @cmarshak @sssangha

bbuzz31

on closer inspection, there needs to be a fail safe to set a default value for interpolation grid in the case that one does not exist in the users .yaml file (when running raider not in the GUNW workfloW)

cmarshak · 2023-08-21T16:39:23Z

on closer inspection, there needs to be a fail safe to set a default value for interpolation grid in the case that one does not exist in the users .yaml file (when running raider not in the GUNW workfloW)

Hey @bbuzz31 - I feel like this ask would be a separate PR. The above PR is really focused on the GUNW. Did this PR change an existing functionality for the non-GUNW workflow? I don't think interpolation is an option (yet) in the non-GUNW workflow. Regardless, this is a pretty big PR largely because of testing.

If you wanted specific changes, please give some reference in the code of what you wanted.

bbuzz31 · 2023-08-21T16:59:12Z

on closer inspection, there needs to be a fail safe to set a default value for interpolation grid in the case that one does not exist in the users .yaml file (when running raider not in the GUNW workfloW)

Hey @bbuzz31 - I feel like this ask would be a separate PR. The above PR is really focused on the GUNW. Did this PR change an existing functionality for the non-GUNW workflow? I don't think interpolation is an option (yet) in the non-GUNW workflow. Regardless, this is a pretty big PR largely because of testing.

If you wanted specific changes, please give some reference in the code of what you wanted.

No it should be in here. There is no default behavior if there isn't an interp_method specified by the user, and no clear warning given. So prior config files will fail for nonGUNW raider workflow. Sections starting on Line 353 in raider.py should set a default method (I forget what we decided on) with a warning message perhaps.

cmarshak · 2023-08-21T17:47:51Z

tools/RAiDER/cli/raider.py

@@ -245,12 +246,27 @@ def calcDelays(iargs=None):
        # 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]


@bbuzz31 - this is the default behavior currently i.e. when the interpolate_time is True or False. If we change this to None by not specifying in yaml, the same behavior will occur. What more are you asking?

cmarshak · 2023-08-21T17:48:30Z

tools/RAiDER/cli/raider.py

@@ -245,12 +246,27 @@ def calcDelays(iargs=None):
        # 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]
+        interp_method = params['interpolate_time']


@bbuzz31 - when you read a yaml file with a given field and it is not specified - then the dictionary returns None for the said blank-valued key so the default behavior is not changed. I will change this to get for further clarity.

cmarshak · 2023-08-21T17:59:03Z

on closer inspection, there needs to be a fail safe to set a default value for interpolation grid in the case that one does not exist in the users .yaml file (when running raider not in the GUNW workfloW)

Hey @bbuzz31 - I feel like this ask would be a separate PR. The above PR is really focused on the GUNW. Did this PR change an existing functionality for the non-GUNW workflow? I don't think interpolation is an option (yet) in the non-GUNW workflow. Regardless, this is a pretty big PR largely because of testing.
If you wanted specific changes, please give some reference in the code of what you wanted.

No it should be in here. There is no default behavior if there isn't an interp_method specified by the user, and no clear warning given. So prior config files will fail for nonGUNW raider workflow. Sections starting on Line 353 in raider.py should set a default method (I forget what we decided on) with a warning message perhaps.

All the tests are passing. I am still unclear what you mean by default for the non-GUNW workflow or how that is called. I am hoping that the non-GUNW workflow just reads or modifies the yaml file here? Is there an example?

My opinion is there should be a test for the control flow that you are asking and therefore a separate PR. I assume you mean the NotImplementedError in Raider here for 353?

This seems problematic to me: https://github.com/dbekaert/RAiDER/blob/dev/tools/RAiDER/cli/raider.py#L286-L289 - unexpected behavior when you specify a non-nearest interpolation and only one file is returned.

I will add some additional control flow, but again - please be as specific as possible with what you want here. It's very intricate control flow and I tried to modify it as little as possible while including these changes.

bbuzz31 · 2023-08-21T18:00:48Z

Lines 251 to 264 in raider.py are as follows

interp_method = params['interpolate_time']
if interp_method in ['none', 'center_time']:
    times = get_nearest_wmtimes(t, [model.dtime() if \
                                model.dtime() is not None else 6][0]) if interp_method == 'center_time' else [t]
elif interp_method == 'azimuth_time_grid':
    n_target_dates = 3
    step = model.dtime()
    time_step_hours = step if step is not None else 6
    # Will yield 3 dates
    times = get_n_closest_datetimes(t,
                                    n_target_dates,
                                    time_step_hours)
else:
    raise NotImplementedError

interp_method just needs to be a string, currently its a boolean when it doesn't exist in the user specified yaml file.
.get('interp_method', 'none') as you suggested should be fine

bbuzz31 · 2023-08-22T17:55:28Z

@cmarshak thanks for taking care of these last bits. feel free to merge if you are finished

cmarshak added 4 commits July 24, 2023 15:42

include azimuth-timing

ab36c3f

relative import

aed06bb

update interp options

1a01b5b

changelog

af0440b

cmarshak changed the title ~~Azimuth time interpolation for Weather ModelS~~ Azimuth time interpolation for Weather Models in GUNW Workflow Jul 25, 2023

cmarshak marked this pull request as draft July 25, 2023 00:08

cmarshak added 2 commits July 25, 2023 13:53

end2end test

afd8495

fix azi test

4e00730

update tests

260b132

cmarshak added 3 commits July 27, 2023 17:35

resolve conflicts

a08c445

weather files

9f34b38

mock weather model files

7ac486d

error handling

ccdaaf6

cmarshak requested review from bbuzz31 and jlmaurer August 14, 2023 21:26

cmarshak marked this pull request as ready for review August 14, 2023 21:27

cmarshak added 2 commits August 14, 2023 14:30

lint

e74c638

mock orbit download from prepGUNW

89a2bbb

cmarshak requested a review from dbekaert August 14, 2023 22:17

cmarshak added 2 commits August 14, 2023 17:30

default azimuth and include hrrrak

a487c59

specify interpolation for GUNW basic test

57e8775

bbuzz31 approved these changes Aug 18, 2023

View reviewed changes

bbuzz31 requested changes Aug 19, 2023

View reviewed changes

cmarshak commented Aug 21, 2023

View reviewed changes

more explicit control flow for interp_method

f20d81b

bbuzz31 approved these changes Aug 21, 2023

View reviewed changes

cmarshak added 2 commits August 21, 2023 12:15

update test for explicit lack of interpolation

1a4476e

more updates for default none time

dd89667

add time to log message with date

e98ca9c

cmarshak merged commit a4679e7 into dev Aug 22, 2023
4 checks passed

cmarshak deleted the azimuth-time-interpolation branch August 22, 2023 19:51

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Azimuth time interpolation for Weather Models in GUNW Workflow #571

Azimuth time interpolation for Weather Models in GUNW Workflow #571

cmarshak commented Jul 25, 2023 •

edited

Loading

cmarshak commented Jul 25, 2023 •

edited

Loading

cmarshak commented Jul 28, 2023

cmarshak commented Aug 14, 2023

bbuzz31 commented Aug 15, 2023

cmarshak commented Aug 15, 2023

cmarshak commented Aug 15, 2023

cmarshak commented Aug 16, 2023

bbuzz31 commented Aug 17, 2023

bbuzz31 commented Aug 18, 2023 •

edited

Loading

bbuzz31 left a comment

bbuzz31 left a comment

cmarshak commented Aug 21, 2023 •

edited

Loading

bbuzz31 commented Aug 21, 2023 •

edited

Loading

cmarshak Aug 21, 2023 •

edited

Loading

cmarshak Aug 21, 2023 •

edited

Loading

cmarshak commented Aug 21, 2023

bbuzz31 commented Aug 21, 2023 •

edited

Loading

bbuzz31 commented Aug 22, 2023

Azimuth time interpolation for Weather Models in GUNW Workflow #571

Azimuth time interpolation for Weather Models in GUNW Workflow #571

Conversation

cmarshak commented Jul 25, 2023 • edited Loading

cmarshak commented Jul 25, 2023 • edited Loading

cmarshak commented Jul 28, 2023

cmarshak commented Aug 14, 2023

bbuzz31 commented Aug 15, 2023

cmarshak commented Aug 15, 2023

cmarshak commented Aug 15, 2023

cmarshak commented Aug 16, 2023

GUNWs used:

bbuzz31 commented Aug 17, 2023

bbuzz31 commented Aug 18, 2023 • edited Loading

bbuzz31 left a comment

Choose a reason for hiding this comment

bbuzz31 left a comment

Choose a reason for hiding this comment

cmarshak commented Aug 21, 2023 • edited Loading

bbuzz31 commented Aug 21, 2023 • edited Loading

cmarshak Aug 21, 2023 • edited Loading

Choose a reason for hiding this comment

cmarshak Aug 21, 2023 • edited Loading

Choose a reason for hiding this comment

cmarshak commented Aug 21, 2023

bbuzz31 commented Aug 21, 2023 • edited Loading

bbuzz31 commented Aug 22, 2023

cmarshak commented Jul 25, 2023 •

edited

Loading

cmarshak commented Jul 25, 2023 •

edited

Loading

bbuzz31 commented Aug 18, 2023 •

edited

Loading

cmarshak commented Aug 21, 2023 •

edited

Loading

bbuzz31 commented Aug 21, 2023 •

edited

Loading

cmarshak Aug 21, 2023 •

edited

Loading

cmarshak Aug 21, 2023 •

edited

Loading

bbuzz31 commented Aug 21, 2023 •

edited

Loading