Release 0.5.0

[mcgarth]

This PR completes the refactor of the PyRate workflow to enable easier analysis of correction data and quicker re-computation with different parameter settings

Added

• New functionality "linear_rate" to calculate linear regression of cumulative displacement time series for every pixel as part of the timeseries step.
• Script for plotting timeseries and linear_rate output geotiff products using Matplotlib. To use, additional dependencies listed in requirements-plot.txt are required.
• Correction data (except maxvar and vcmt) applied to the ifg data is saved to disk and re-used on subsequent repeat runs. Corrections are only re-calculated if config parameters change between runs.
• MPI parallelisation of APS spatio-temporal filter correction.
• Unit test coverage for refpixel lat/lon to x/y conversion and aps module.

Fixed

• Re-enable ifglksx and ifglksy to be different values, resulting in different resolutions in x and y dimensions in multi-looked interferograms.
• Re-enable orbfitlksx and orbfitlksy to be different values, resulting in different resolutions in x and y dimensions during network orbit correction.
• Screen messages from main process only during MPI runs.

Changed

• process step has been renamed correct. Stacking and timeseries have been removed from this step and are now invoked by separate timeseries and stack command line options.
• Processing of coherence files by conv2tif and prepifg is now triggered by the presence of cohfilelist in the config file. If the list is present, multilooked/cropped coherence files are saved to disk, regardless of whether cohmask is 0 or 1.
• Parallelisation capability is refactored - MPI and multiprocessing both now use a common tiling framework for stack, timeseries and mst algorithms.
• Introduced a simplified and standardised file naming format for files produced by the prepifg step. Info previously saved in the filename (e.g. ifglksx, ifglksy, ifgcropopt) is now saved to the geotiff header instead.

Removed

• Redundant tscal config parameter was deprecated - not needed now there is a timeseries step invokable on the command line.
• Pillow, cython and glob2 dependencies - not used in main workflow code
• Deprecated function pyrate.prepifg_helper.prepare_ifgs, which is no longer needed.

[codecov-commenter]

Codecov Report

Merging #295 into master will decrease coverage by 7.48%. The diff coverage is 70.21%.

@@            Coverage Diff             @@
##           master     #295      +/-   ##
==========================================
- Coverage   90.71%   83.22%   -7.49%     
==========================================
  Files          26       26              
  Lines        3305     3470     +165     
  Branches      514      544      +30     
==========================================
- Hits         2998     2888     -110     
- Misses        215      486     +271     
- Partials       92       96       +4     

Impacted Files	Coverage Δ
pyrate/default_parameters.py	100.00% <ø> (ø)
pyrate/merge.py	16.16% <5.45%> (-79.04%)	:arrow_down:
pyrate/main.py	23.25% <19.44%> (-75.14%)	:arrow_down:
pyrate/core/logger.py	36.36% <33.33%> (-16.58%)	:arrow_down:
pyrate/core/stack.py	75.00% <40.00%> (-22.81%)	:arrow_down:
pyrate/core/prepifg_helper.py	93.52% <50.00%> (-1.65%)	:arrow_down:
pyrate/prepifg.py	50.74% <59.09%> (-2.43%)	:arrow_down:
pyrate/core/timeseries.py	77.24% <59.67%> (-12.33%)	:arrow_down:
pyrate/core/shared.py	93.14% <64.28%> (-2.12%)	:arrow_down:
pyrate/conv2tif.py	94.44% <66.66%> (-5.56%)	:arrow_down:
... and 22 more

---

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 f77ad6e...df7cb1e. Read the comment docs.

[tfuhrmann]

I've done some testing of the orbit correction today. First of all: it's great that the results are now saved into .npy files and we can actually have a look at the correction values!

Three issues to be solved or discussed before we can release:

1. The setting of the new parameter 'orbfitoffset' is not used in the file naming in the orb_error directory. That means, when running two successive runs one with 'orbfitoffset: 0' and one with 'orbfitoffset: 1', the orbit correction is not re-done and the following warning appears on the screen: 15:02:07 orbital 269 31029 WARNING 0/7 Reusing orbfit errors from previous run!!!
2. The outputs of the independent correction (method 1) are fine and what I would expect for the different polynomial degrees (see images below). I would like to understand though what the 'orbfitoffset' does and why it changes the direction of the planar phase ramp and the orientation of the quadratic "ramp".
3. The outputs of the network method are not what I would expect. All npy-files in the orb_error directory contain zero elements only (in case 'orbfitoffset: 0') or equal values for all elements (in case 'orbfitoffset: 1'). FYI: I've used the following python code snippet to verify this observation:

   
   import numpy as np
   import glob

for file in glob.glob('./out/orb_error/*.npy'): print(file) orb = np.load(file) if np.all(orb == orb[0]): if np.all(orb == 0): print('all elements in numpy array are zero') else: print('all elements in numpy array are equal')

Plots of the orbit correction values using different polynomial degrees a) planar, 'orbfitoffset: 0', b) planar, 'orbfitoffset: 1', c) quadratic, 'orbfitoffset: 0', d) quadratic, 'orbfitoffset: 1', e) part-cubic, 'orbfitoffset: 0', f) part-cubic, 'orbfitoffset: 1':
![orb_corr_1planar](https://user-images.githubusercontent.com/20388928/91400218-a16eb100-e882-11ea-8825-d04a0fa5cc82.png)
![orb_corr_1planar_offset0](https://user-images.githubusercontent.com/20388928/91400228-a59ace80-e882-11ea-9627-f937e29b0345.png)
![orb_corr_2quadratic](https://user-images.githubusercontent.com/20388928/91400232-a6cbfb80-e882-11ea-91bc-dd2f2a668fba.png)
![orb_corr_2quadratic_offset0](https://user-images.githubusercontent.com/20388928/91400233-a7fd2880-e882-11ea-8e2e-39fff76bdde9.png)
![orb_corr_3part-cubic](https://user-images.githubusercontent.com/20388928/91400236-a895bf00-e882-11ea-801b-7dfec2caab81.png)
![orb_corr_3part-cubic_offset0](https://user-images.githubusercontent.com/20388928/91400240-a9c6ec00-e882-11ea-9a53-b94f95e4af72.png)

**EDIT** 

Min/max values of the different orbit corrections:

method 1, offset 1 (IFG1)
planar:      -0.7990   -0.2056 
quadratic:    3.7605   19.9497 
part-cubic: -12.3915    1.6210 

method 1, offset 0 (IFG1)
planar:      -0.3140    0.0990
quadratic:   -9.4310    1.4995
part-cubic: -12.3913    1.6211

method 2, offset 1
IFG1: 0.0672    0.0672 
IFG2: 27.6266   27.6266
-> identical values for all three polynomial methods

method 2, offset 0
-> zero for all methods and IFGs

[tfuhrmann]

I've tested the new plot_time_series.py tool and can confirm that it works. As I already said in our meeting: great work, @chandra2ga ! That's how I ran it (inside a PyRateVenv):

pip install -r requirements-plot.txt
python3 utils/plot_time_series.py /g/data/dg9/INSAR_ANALYSIS/CAMDEN/RSAT2/GAMMA/T081D/pyrate/out

Has a description of this tool been added to the PyRate documentation?

I'm happy to have the tool released as is. Some items/ideas for implementation in the next release:

• change coordinates shown on the spatial plot from pixel row/col to lat/lon.
• change date format in the time bar below the spatial plot and in the time series plot (e.g. use a MM/YY format such as Jan16 Jul16 Jan17 Jul17 ...).
• always use a colour scale centered around zero (i.e. pixels with no/small movement are displayed in white)
• I've tried the right-drag to change the reference area. While the Ref area X/Y values left of the spatial image change, the velocities in the plot remain the same. So, I guess this feature is currently not fully implemented. Can that be easily done or should we remove it altogether?

[tfuhrmann]

Last comment from me. I've now also tested the APS correction and can confirm that it works as expected. Below is a visualisation of the original interferometric phase and the APS correction values for the same interferogram. As you can see, the APS correction does filter out some real deformation in this case (i.e. the red spots which are located at underground coal mines). The good thing is: we can now visualise the APS correction and play around with the parameters used for spatio-temporal filtering to help avoid filtering out signals which can be clearly attributed to deformation. Note that all images shown in my comments were generated with Matlab (mainly because I'm so much quicker in doing such things in Matlab :-). For the next release, it would be great if we had a Python tool to generate such plots for particular/all interferograms from the .npy files.

[mcgarth]

Thanks for the review @adeane-ga . I have just pushed two commits that address:

• enabling different orbitfitlksx/y values (error messages and process fail), see below ERROR MESSAGES 1 and 2.
• perhaps there could be more explanation in the terminal -h display regarding the timeseries and stack options. At the moment, a user might assume they are both required and in the stated order.

I believe the test errors you report (as does @tfuhrmann ) are a problem with running the full test suite on gadi because the full test suite is passing for four builds in Travis. Is this your thinking @basaks ?

[mcgarth]

Last comment from me. I've now also tested the APS correction and can confirm that it works as expected. Below is a visualisation of the original interferometric phase and the APS correction values for the same interferogram. As you can see, the APS correction does filter out some real deformation in this case (i.e. the red spots which are located at underground coal mines). The good thing is: we can now visualise the APS correction and play around with the parameters used for spatio-temporal filtering to help avoid filtering out signals which can be clearly attributed to deformation. Note that all images shown in my comments were generated with Matlab (mainly because I'm so much quicker in doing such things in Matlab :-). For the next release, it would be great if we had a Python tool to generate such plots for particular/all interferograms from the .npy files.

Good to hear the APS correction is doing something sensible. An action for the next release is to further test this functionality to identify reasonable parameter values and also improve unit test coverage.

[basaks]

• enabling different orbitfitlksx/y values (error messages and process fail), see below ERROR MESSAGES 1 and 2.
• perhaps there could be more explanation in the terminal -h display regarding the timeseries and stack options. At the moment, a user might assume they are both required and in the stated order.

I will address this.

I believe the test errors you report (as does @tfuhrmann ) are a problem with running the full test suite on gadi because the full test suite is passing for four builds in Travis. Is this your thinking @basaks ?

Yes, that's right. We have certain restriction in gadi, which throws those errors in gadi becasue of the way the tests are written.

[mcgarth]

I will address this.

No need @basaks! I did it in be39624 and 55aee31 ;-)

[mcgarth]

Thanks for this thorough testing @tfuhrmann! It has uncovered an unknown issue with the network orbital method. I am documenting here what we discussed verbally as a team

1. The setting of the new parameter 'orbfitoffset' is not used in the file naming in the orb_error directory. That means, when running two successive runs one with 'orbfitoffset: 0' and one with 'orbfitoffset: 1', the orbit correction is not re-done and the following warning appears on the screen: 15:02:07 orbital 269 31029 WARNING 0/7 Reusing orbfit errors from previous run!!!

We agreed that orbfitoffset parameter should remain hidden from the user (i.e. removed from config files). Furthermore, orbfitoffset should be hard-coded off for the independent orbit method and on for the network orbit method. Because it will not be a variable, the parameter orbfitoffset can stay out of the *orbfit.npy filenames.

2. The outputs of the independent correction (method 1) are fine and what I would expect for the different polynomial degrees (see images below). I would like to understand though what the 'orbfitoffset' does and why it changes the direction of the planar phase ramp and the orientation of the quadratic "ramp".

This observation is a bit worrying. I think to resolve this we need a new synthetic interferogram dataset that contains only simulated orbital errors and include that as a new unit test dataset that is used in the test suite. I'll add a ticket...

3. The outputs of the network method are not what I would expect. All npy-files in the orb_error directory contain zero elements only (in case 'orbfitoffset: 0') or equal values for all elements (in case 'orbfitoffset: 1').

I think there are some deep-seated issues here that need to be investigated. In the interim, we agreed that network orbital method will be disabled for release-0.5.0. Users selecting network method will be informed (via log.warning) of the reason for disablement and PyRate will automatically use independent method instead. The above-suggested synthetic interferograms could help us fix the issues with network method.

[adeane-ga]

@mcgarth : I have just tested your commits relating to the orbfitlksx/y parameter and the changes to the pyrate -h display - both of these changes are working and I am not getting the process fail or error messages that I was getting before. Thanks.

[tfuhrmann]

I haven't had a look at the orbfit code in PyRate, but I guess it would be similar to the below Matlab code snippet from StaMPS. It's bascially a 2D regression, e.g. like this for the case of a 2D linear phase ramp with the m_i parameters describing the spatial trend:

The code might be good for a cross-check with what's done in PyRate:

if scla_deramp>0
    fprintf('   deramping ifgs...\n')
    ph_ramp=zeros(ps.n_ps,n_ifg,'single');
    x=ps.xy(:,2);
    y=ps.xy(:,3);
    if scla_deramp==1
     % linear ramp
     fprintf('   linear ramp is estimated...\n')
     G=double([ones(ps.n_ps,1),x,y]);
    end
    if scla_deramp==2
     % quadratic "ramp"
     fprintf('   quadratic function is estimated...\n')
     G=double([ones(ps.n_ps,1),x,y,x.*y,x.^2,y.^2]);
    end
    if scla_deramp==3
     % cubic "ramp"
     fprintf('   cubic function is estimated...\n')
     G=double([ones(ps.n_ps,1),x,y,x.*y,x.^2,y.^2,x.^2.*y,x.*y.^2,x.^3,y.^3]);
    end
    for i=1:length(unwrap_ifg_index)
        d=uw.ph_uw(:,unwrap_ifg_index(i));
        % don't use NaN values for ramp estimation
        nanix = find(isnan(d));
        if ~isempty(nanix)
            d(nanix) = 0; % use zero for matrix multiplication
        end        
        m=G\double(d(:));
        ph_this_ramp=G*m;
        uw.ph_uw(:,unwrap_ifg_index(i))=uw.ph_uw(:,unwrap_ifg_index(i))-ph_this_ramp; % subtract ramp
        ph_ramp(:,unwrap_ifg_index(i))=ph_this_ramp;
        uw.ph_uw(nanix,unwrap_ifg_index(i)) = NaN;
        ph_ramp(nanix,unwrap_ifg_index(i)) = NaN;
    end
else
    ph_ramp=[];
end

[mcgarth]

The name offset is then a bit misleading. Rather it should be intercept. I think it is clear that we want this parameter to be always on (True). Actions:

• remove orbfitoffset from user control
• set orbfitoffset to True for both network and independent methods.
• consider deprecating orbfitoffset as a (albeit hidden) parameter in the code altogether in the future.

First two points are actioned in the merged PR commit e6cee7a

[mcgarth]

@tfuhrmann @chandra2ga - commit ee39524 adds some information to the documentation on how to use the new plotting script

[mcgarth]

Thanks for your contributions to this release @basaks @chandra2ga @adeane-ga @tfuhrmann !