Sb/phase closure loops

[basaks]

This PR adds phase closure correction in PyRate .

This readme can also be found here: https://github.com/GeoscienceAustralia/PyRate/tree/6dcac13d26d059101ac67ae54f69b3a8b32b64e6/pyrate/core/phase_closure

Update the doc here to track changes in pyrate repo.

Phase Closure

To use phase closure correction, simply add the string _phaseclosure in the list of corrections as can be seen below:

[correct]
steps =
    orbfit
    refphase
    demerror
    phase_closure
    mst
    apscorrect
    maxvar

The following params are used from the pyrate config file:

# phase closure params
# large_deviation_threshold_for_pixel # pi/2, pi=3.1416
# threshold_to_remove_ifg: ifgs with more than this fraction of pixels with error will be dropped
# loop_count_for_threshold_to_remove_ifg: pixel with phase unwrap error in at least this many loops
# phase_unwrap_error_threshold: pixel with phase unwrap error in more than this many ifgs will be flagged
# max_loop_length: loops upto this many edges are considered for closure checks
# subtract_median_in_closure_check: whether to subtract median during closure checks
# max_loop_count_for_each_ifgs: loops with more than these many ifgs are discarded.
# max_loop_count_for_each_ifgs: Ifg count must be met for all ifgs in the loop for loop to be discarded
large_dev_thr: 1.5708
threshold_to_remove_ifg: 0.07
loop_count_for_threshold_to_remove_ifg: 2
phase_unwrap_error_threshold: 5
max_loop_length: 3
subtract_median_in_closure_check: 1
max_loop_count_for_each_ifgs: 2

Currently only max_loop_length: 3 produce repeatable results with Mexico CropA dataset.

Phase closure correction has the following main functionalities:

1. Compute the closed loops using networkx. Loops are assigned signs for each interferogram, and assigned a weight based on total weight of each loop, which is the sum of difference between the ifg second and first date. This is done in python file _mstclosure.py. We perform several steps in this stage:

   1. Discard loops that are more than _max_looplength.
   2. Sort each loop based on first date of each interferogram (lower weight first). In case of a tie, we sort based on the second date of the interferograms.
   3. Compute weight of each interferograms (=second date -first date).
   4. Then we sum the weights of interferograms in a loop to find the weight of each closed loop.
   5. Sort the loops based on weights. In case of ties, we further sort by primary dates, and then by secondary dates.
   6. Discard loops containing _max_loop_count_for_eachifgs. All ifgs in the loop must have contributed to at least _max_loop_count_for_eachifgs.
   7. Drop ifgs not part of any loop after this stage.

2. compute _sumclosure of each loop from stage 1 for each pixel. In addition, we perform the following steps in _sumclosure.py:

   1. Find pixels breaching _large_devthr. Note _large_devthr is specified in radians. Therefore at this stage we need to convert phase data (in millimeters) into radians (check functions _shared.convert_toradian and it's use in the Ifg class).
   2. In order to find Persistent Scatter, compute the _checkps for each pixel for each ifg.
   3. See use of _subtract_median_in_closurecheck in function _sumclosure.__compute_checkps.

3. _closurecheck.py is used for orchestration of the functionalities above. After stage 2, we drop ifgs exceeding _threshold_to_removeifg and _loop_count_for_threshold_to_removeifg. See docstring in function _closure_check.drop_ifgs_exceedingthreshold.

4. Steps 1-3 are repeated until a stable list of interferograms are returned (see
   _closure_check.filter_to_closure_checkedifgs).

5. Once a stable list of interferograms is found, in correct.py, write a new ifglist in the working directory, update params, and use the updated ifglist for further pyrate processing.

6. Also in correct.py, we detect persistent scatterers (pixels) breaching _phase_unwrap_errorthreshold.

7. Finally, we write the ifg phase data after assigning nulls to pixels breaching _phase_unwrap_errorthreshold. Phase closure correction is done at this stage.

[codecov-io]

Codecov Report

Merging #307 into develop will increase coverage by 0.23%. The diff coverage is 95.68%.

@@             Coverage Diff             @@
##           develop     #307      +/-   ##
===========================================
+ Coverage    83.22%   83.45%   +0.23%     
===========================================
  Files           27       28       +1     
  Lines         3845     3923      +78     
  Branches       588      599      +11     
===========================================
+ Hits          3200     3274      +74     
- Misses         522      526       +4     
  Partials       123      123              

Impacted Files	Coverage Δ
pyrate/core/phase_closure/mst_closure.py	91.52% <91.52%> (ø)
pyrate/core/mst.py	96.61% <100.00%> (-0.12%)	:arrow_down:
pyrate/core/orbital.py	93.65% <100.00%> (+0.20%)	:arrow_up:
pyrate/core/shared.py	92.42% <100.00%> (ø)
pyrate/core/prepifg_helper.py	94.96% <0.00%> (+1.43%)	:arrow_up:

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update b423291...92160cc. Read the comment docs.

[mcgarth]

I too have done some initial testing with the cropA Mexico data. Turning coherence masking on and off has a big impact on the result of closure check. With cohmask: 0 I get:

15:09:34 closure_check 68 904571 INFO 0/0 Performing closure check on original set of 30 ifgs
15:09:37 closure_check 78 904571 INFO 0/0 After closure check 27 ifgs are retained
15:09:38 correct 92 904571 INFO 0/0 Found 13 unique epochs in the 27 interferogram network

i.e. 3 ifgs are dropped.

But with cohmask: 1 and cohthresh: 0.3 I get:

14:57:36 closure_check 68 823640 INFO 0/0 Performing closure check on original set of 30 ifgs
14:57:39 closure_check 78 823640 INFO 0/0 After closure check 23 ifgs are retained
14:57:39 correct 92 823640 INFO 0/0 Found 11 unique epochs in the 23 interferogram network

i.e. 7 ifgs are dropped, and 2 epochs are lost.

I wonder if coherence masking should be happening after phase closure checking has been applied... thoughts?

[mcgarth]

Regarding a strategy for reducing the number of loops used, @tfuhrmann , @basaks and I agreed on implementing:

• a weight metric for each loop, suggested as the summation of temporal baselines for each interferogram forming that loop (e.g. 12 days + 12 days + 24 days = 48 days), being the minimum weight for a 12 day revisit sensor like Sentinel-1.
• then prioritise loops based on low weight. in this way, loops involving nearby epochs would be prioritised (relatively low weight), as would loops with the minimum number of interferograms (3).
• then use these loop weights in combination with user-definable parameters governing the min (e.g. 2) and max (e.g. 5) number of times an interferogram should feature in a loop, to select the subset of loops for the closure check.

[mcgarth]

We have agreed to close (merge) this PR once CI is passing in order to facilitate merging a bunch of backed-up features in to develop. Despite this, review of the phase closure functionality is not fully complete, and will be subject to a further PR and review prior to making it in to a PyRate release.