New residual method ("window" instead of "cumulative"). This calcault…

…es the residual as the ratio of the maximum RMS for all windows / mean RMS for all windows of the cumulative residual. Because of this, it doesn't grow in the same way the cumulative residual does. Note also that this could be changed to be as above but for the residual, rather than the cumulative residual.
matthew-gaddes · Dec 1, 2021 · 59142d3 · 59142d3
1 parent a083644
commit 59142d3
Show file tree

Hide file tree

Showing 34 changed files with 123 additions and 27 deletions.
diff --git a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/01_pca_variance_line.png b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/01_pca_variance_line.png
diff --git a/...Alert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/02_PCA_sources_and_tcs.png b/...Alert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/02_PCA_sources_and_tcs.png
diff --git a/..._01_Sierra_Negra_no_intermediate/ICASAR_outputs/02_PCA_sources_correlations.png b/..._01_Sierra_Negra_no_intermediate/ICASAR_outputs/02_PCA_sources_correlations.png
diff --git a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/02_PCA_sources_time.png b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/02_PCA_sources_time.png
diff --git a/..._01_Sierra_Negra_no_intermediate/ICASAR_outputs/03_ICA_sources_correlations.png b/..._01_Sierra_Negra_no_intermediate/ICASAR_outputs/03_ICA_sources_correlations.png
diff --git a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/03_ICA_sources_time.png b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/03_ICA_sources_time.png
diff --git a/...t_01_Sierra_Negra_no_intermediate/ICASAR_outputs/03_PCA_source_correlations.png b/...t_01_Sierra_Negra_no_intermediate/ICASAR_outputs/03_PCA_source_correlations.png
diff --git a/...rra_Negra_no_intermediate/ICASAR_outputs/04_clustering_and_manifold_results.png b/...rra_Negra_no_intermediate/ICASAR_outputs/04_clustering_and_manifold_results.png
diff --git a/...ert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/05_ICASAR_sourcs_and_tcs.png b/...ert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/05_ICASAR_sourcs_and_tcs.png
diff --git a/...t_01_Sierra_Negra_no_intermediate/ICASAR_outputs/06_ICA_source_correlations.png b/...t_01_Sierra_Negra_no_intermediate/ICASAR_outputs/06_ICA_source_correlations.png
diff --git a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/ICASAR_results.pkl b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/ICASAR_results.pkl
diff --git a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/ICs.kmz b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/ICs.kmz
diff --git a/...a_Negra_no_intermediate/LiCSAlert_figure_with_063_monitoring_interferograms.png b/...a_Negra_no_intermediate/LiCSAlert_figure_with_063_monitoring_interferograms.png
diff --git a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/01_pca_variance_line.png b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/01_pca_variance_line.png
diff --git a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/02_PCA_sources_and_tcs.png b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/02_PCA_sources_and_tcs.png
diff --git a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/02_PCA_sources_correlations.png b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/02_PCA_sources_correlations.png
diff --git a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/02_PCA_sources_time.png b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/02_PCA_sources_time.png
diff --git a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/03_ICA_sources_correlations.png b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/03_ICA_sources_correlations.png
diff --git a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/03_ICA_sources_time.png b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/03_ICA_sources_time.png
diff --git a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/03_PCA_source_correlations.png b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/03_PCA_source_correlations.png
diff --git a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/04_clustering_and_manifold_results.png b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/04_clustering_and_manifold_results.png
diff --git a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/05_ICASAR_sourcs_and_tcs.png b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/05_ICASAR_sourcs_and_tcs.png
diff --git a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/06_ICA_source_correlations.png b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/06_ICA_source_correlations.png
diff --git a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/ICASAR_results.pkl b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/ICASAR_results.pkl
diff --git a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/ICs.kmz b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/ICs.kmz
diff --git a/LiCSAlert_03_Campi_Flegrei/LiCSAlert_figure_with_227_monitoring_interferograms.png b/LiCSAlert_03_Campi_Flegrei/LiCSAlert_figure_with_227_monitoring_interferograms.png
diff --git a/LiCSAlert_03_Campi_Flegrei/LiCSAlert_results.pkl b/LiCSAlert_03_Campi_Flegrei/LiCSAlert_results.pkl
diff --git a/LiCSAlert_batch_mode_examples.py b/LiCSAlert_batch_mode_examples.py
@@ -15,8 +15,8 @@
 import licsalert
 from licsalert.licsalert import LiCSAlert_batch_mode
 
-ICASAR_path = Path("/home/matthew/university_work/15_my_software_releases/ICASAR-2.7.2/")                               # location of ICASAR functions
-#ICASAR_path = Path("/home/matthew/university_work/01_blind_signal_separation_python/13_ICASAR/ICASAR_GitHub")           # development version
+#ICASAR_path = Path("/home/matthew/university_work/15_my_software_releases/ICASAR-2.7.2/")                               # location of ICASAR functions
+ICASAR_path = Path("/home/matthew/university_work/01_blind_signal_separation_python/13_ICASAR/ICASAR_GitHub")           # development version
 
 
 #%% Load Sentinel-1 data for Sierra Negra, note that this wasn't processed with LiCSBAS, and doesn't include the DEM.  
@@ -49,7 +49,8 @@
                       "run_ICASAR" : True,                                           # If False, attempt to load results from previous run.  If True, run (which can be slow)
                       "intermediate_figures" : False,                                # if set to True, a figure is produced for all time steps in the monitoring data, which can be time consuming.  
                       "downsample_run" : 0.5,                                        # data can be downsampled to speed things up
-                      "downsample_plot" : 0.5}                                       # and a 2nd time for fast plotting.  Note this is applied to the restuls of the first downsampling, so is compound
+                      "downsample_plot" : 0.5,                                       # and a 2nd time for fast plotting.  Note this is applied to the restuls of the first downsampling, so is compound,
+                      "residual_type"        : 'cumulative'}                         # this is the default from the 2019 paper, but it could be switched to "window", but this hasn't been tested widely.  
 
 
 ICASAR_settings = {"n_comp" : 6,                                                     # number of components to recover with ICA (ie the number of PCA sources to keep)
@@ -59,7 +60,7 @@
                     "ica_param" : (1e-2, 150),                                       # (tolerance, max iterations)
                     "figures" : "png",                                               # if png, saved in a folder as .png.  If window, open as interactive matplotlib figures, if window+png then both.  
                     "create_all_ifgs_flag" : False,                                  # Creates all possible pairs of ifgs between all acquisitions.  Results can be more complex with this set to True, but also better at recovering small magnitude signals
-                    "load_fastICA_results" : False}                                   # If True, ICASAR will try to load the results of FastICA from previous runs.  This is useful if you wish to fine tune the hdbscan or tsne settings quickly.     
+                    "load_fastICA_results" : True}                                   # If True, ICASAR will try to load the results of FastICA from previous runs.  This is useful if you wish to fine tune the hdbscan or tsne settings quickly.     
 
 
 
@@ -74,6 +75,8 @@
 LiCSAlert_batch_mode(displacement_r2_copy, ICASAR_settings = ICASAR_settings, **LiCSAlert_settings,ICASAR_path = ICASAR_path)
 
 
+
+
 #%% Example 3, Running LiCSAlert with a smaller signal (Campi Flegrei, processed with LiCSBAS)
 
 LiCSBAS_out_folder_campi_flegrei = Path('./022D_04826_121209')
@@ -86,7 +89,8 @@
                       "run_ICASAR" : True,                                           # If False, attempt to load results from previous run.  If True, run (which can be slow)
                       "intermediate_figures" : False,                                # if set to True, a figure is produced for all time steps in the monitoring data, which can be time consuming.  
                       "downsample_run" : 0.5,                                        # data can be downsampled to speed things up
-                      "downsample_plot" : 0.5}                                       # and a 2nd time for fast plotting.  Note this is applied to the restuls of the first downsampling, so is compound
+                      "downsample_plot" : 0.5,                                       # and a 2nd time for fast plotting.  Note this is applied to the restuls of the first downsampling, so is compound
+                      "residual_type"        : 'cumulative'}                              # controls the type of residual used in the lower plot.  Either cumulative or window   
 
 
 ICASAR_settings = {"n_comp" : 5,                                         # number of components to recover with ICA (ie the number of PCA sources to keep)
@@ -95,11 +99,11 @@
                    "ica_param" : (1e-2, 150),                           # (tolerance, max iterations)
                    "hdbscan_param" : (100,10),                           # (min_cluster_size, min_samples) Discussed in more detail in Mcinnes et al. (2017). min_cluster_size sets the smallest collection of points that can be considered a cluster. min_samples sets how conservative the clustering is. With larger values, more points will be considered noise. 
                    "create_all_ifgs_flag" : True,                       # small signals are hard for ICA to extact from time series, so make it easier by creating all possible long temporal baseline ifgs from the incremental data.  
-                   "load_fastICA_results" : False,                      # If all the FastICA runs already exisit, setting this to True speeds up ICASAR as they don't need to be recomputed.  
+                   "load_fastICA_results" : True,                      # If all the FastICA runs already exisit, setting this to True speeds up ICASAR as they don't need to be recomputed.  
                    "figures" : "png+window"}                            # if png, saved in a folder as .png.  If window, open as interactive matplotlib figures,
 
 
-displacement_r2, tbaseline_info = LiCSBAS_to_ICASAR(LiCSBAS_out_folder_campi_flegrei, figures=True)        # open various LiCSBAS products, spatial ones in displacement_r2, temporal ones in tbaseline_info
+displacement_r2, tbaseline_info, _ = LiCSBAS_to_ICASAR(LiCSBAS_out_folder_campi_flegrei, figures=True)        # open various LiCSBAS products, spatial ones in displacement_r2, temporal ones in tbaseline_info
 displacement_r2['ifg_dates'] = tbaseline_info['ifg_dates']                                                  # Unlike ICASAR, LiCSAlert always needs the ifg_dates too.  
 
 LiCSAlert_batch_mode(displacement_r2, ICASAR_settings = ICASAR_settings, **LiCSAlert_settings, ICASAR_path = ICASAR_path)

diff --git a/licsalert/__pycache__/aux.cpython-37.pyc b/licsalert/__pycache__/aux.cpython-37.pyc
diff --git a/licsalert/__pycache__/licsalert.cpython-37.pyc b/licsalert/__pycache__/licsalert.cpython-37.pyc
diff --git a/licsalert/__pycache__/monitoring_functions.cpython-37.pyc b/licsalert/__pycache__/monitoring_functions.cpython-37.pyc
diff --git a/licsalert/aux.py b/licsalert/aux.py
@@ -8,6 +8,36 @@
 
 #%%
 
+
+
+def r2_to_r3(ifgs_r2, mask):
+    """ Given a rank2 of ifgs as row vectors, convert it to a rank3.   Copied from insar_tools.general to avoid making insar_tools a dependency.  
+    Inputs:
+        ifgs_r2 | rank 2 array | ifgs as row vectors 
+        mask | rank 2 array | to convert a row vector ifg into a rank 2 masked array        
+    returns:
+        phUnw | rank 3 array | n_ifgs x height x width
+    History:
+        2020/06/10 | MEG  | Written
+    """
+    import numpy as np
+    import numpy.ma as ma
+    from small_plot_functions import col_to_ma
+
+    n_ifgs = ifgs_r2.shape[0]
+    ny, nx = col_to_ma(ifgs_r2[0,], mask).shape                                   # determine the size of an ifg when it is converter from being a row vector
+
+    ifgs_r3 = np.zeros((n_ifgs, ny, nx))                                                # initate to store new ifgs
+    for ifg_n, ifg_row in enumerate(ifgs_r2):                                           # loop through all ifgs
+        ifgs_r3[ifg_n,] = col_to_ma(ifg_row, mask)                                  
+
+    mask_r3 = np.repeat(mask[np.newaxis,], n_ifgs, axis = 0)                            # expand the mask from r2 to r3
+    ifgs_r3_ma = ma.array(ifgs_r3, mask = mask_r3)                                      # and make a masked array    
+    return ifgs_r3_ma
+
+
+#%%
+
 def get_baseline_end_ifg_n(LiCSBAS_imdates, baseline_end):
     """ Given a list of the dates that there are steps in the LICSBAS time series for (i.e. when there was a Sentinel-1 acquisition),
     find which number is the last before the baseline stage ends.  Note the baseline stage can end on any date (i.e. not one when there's an acquisition)