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 index b44fa488..a131c849 100644 Binary files a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/01_pca_variance_line.png and b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/01_pca_variance_line.png differ diff --git a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/02_PCA_sources_and_tcs.png b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/02_PCA_sources_and_tcs.png deleted file mode 100644 index 5c938c54..00000000 Binary files a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/02_PCA_sources_and_tcs.png and /dev/null differ diff --git a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/02_PCA_sources_correlations.png b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/02_PCA_sources_correlations.png new file mode 100644 index 00000000..b469963b Binary files /dev/null and b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/02_PCA_sources_correlations.png differ 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 new file mode 100644 index 00000000..d755c549 Binary files /dev/null and b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/02_PCA_sources_time.png differ diff --git a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/03_ICA_sources_correlations.png b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/03_ICA_sources_correlations.png new file mode 100644 index 00000000..b2fe4c98 Binary files /dev/null and b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/03_ICA_sources_correlations.png differ 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 new file mode 100644 index 00000000..44b4a1b8 Binary files /dev/null and b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/03_ICA_sources_time.png differ diff --git a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/03_PCA_source_correlations.png b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/03_PCA_source_correlations.png deleted file mode 100644 index 363b29ab..00000000 Binary files a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/03_PCA_source_correlations.png and /dev/null differ diff --git a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/04_clustering_and_manifold_results.png b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/04_clustering_and_manifold_results.png index 332d3bec..f6371092 100644 Binary files a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/04_clustering_and_manifold_results.png and b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/04_clustering_and_manifold_results.png differ diff --git a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/05_ICASAR_sourcs_and_tcs.png b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/05_ICASAR_sourcs_and_tcs.png deleted file mode 100644 index 644ee164..00000000 Binary files a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/05_ICASAR_sourcs_and_tcs.png and /dev/null differ diff --git a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/06_ICA_source_correlations.png b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/06_ICA_source_correlations.png deleted file mode 100644 index 367ed973..00000000 Binary files a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/06_ICA_source_correlations.png and /dev/null differ 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 index 488ed526..98a50fd6 100644 Binary files a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/ICASAR_results.pkl and b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/ICASAR_results.pkl differ diff --git a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/ICs.kmz b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/ICs.kmz index 971f16af..09ba901b 100644 Binary files a/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/ICs.kmz and b/LiCSAlert_01_Sierra_Negra_no_intermediate/ICASAR_outputs/ICs.kmz differ diff --git a/LiCSAlert_01_Sierra_Negra_no_intermediate/LiCSAlert_figure_with_063_monitoring_interferograms.png b/LiCSAlert_01_Sierra_Negra_no_intermediate/LiCSAlert_figure_with_063_monitoring_interferograms.png index 00c5ffc1..a0f18cce 100644 Binary files a/LiCSAlert_01_Sierra_Negra_no_intermediate/LiCSAlert_figure_with_063_monitoring_interferograms.png and b/LiCSAlert_01_Sierra_Negra_no_intermediate/LiCSAlert_figure_with_063_monitoring_interferograms.png differ 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 index 0e8bcdb4..d539d8ce 100644 Binary files a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/01_pca_variance_line.png and b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/01_pca_variance_line.png differ 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 deleted file mode 100644 index 107948a3..00000000 Binary files a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/02_PCA_sources_and_tcs.png and /dev/null differ 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 new file mode 100644 index 00000000..bebb2199 Binary files /dev/null and b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/02_PCA_sources_correlations.png differ 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 new file mode 100644 index 00000000..350874ae Binary files /dev/null and b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/02_PCA_sources_time.png differ 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 new file mode 100644 index 00000000..be6e1922 Binary files /dev/null and b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/03_ICA_sources_correlations.png differ 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 new file mode 100644 index 00000000..6e0b79d4 Binary files /dev/null and b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/03_ICA_sources_time.png differ 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 deleted file mode 100644 index 452958f6..00000000 Binary files a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/03_PCA_source_correlations.png and /dev/null differ 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 index 65828a1f..f104c524 100644 Binary files a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/04_clustering_and_manifold_results.png and b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/04_clustering_and_manifold_results.png differ 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 deleted file mode 100644 index b423a620..00000000 Binary files a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/05_ICASAR_sourcs_and_tcs.png and /dev/null differ 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 deleted file mode 100644 index 4d7b1432..00000000 Binary files a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/06_ICA_source_correlations.png and /dev/null differ diff --git a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/ICASAR_results.pkl b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/ICASAR_results.pkl index 380a09fb..94440312 100644 Binary files a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/ICASAR_results.pkl and b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/ICASAR_results.pkl differ diff --git a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/ICs.kmz b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/ICs.kmz index 72ca431b..a60cb05a 100644 Binary files a/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/ICs.kmz and b/LiCSAlert_03_Campi_Flegrei/ICASAR_outputs/ICs.kmz differ 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 index 92f97fa8..0f7fe25d 100644 Binary files a/LiCSAlert_03_Campi_Flegrei/LiCSAlert_figure_with_227_monitoring_interferograms.png and b/LiCSAlert_03_Campi_Flegrei/LiCSAlert_figure_with_227_monitoring_interferograms.png differ diff --git a/LiCSAlert_03_Campi_Flegrei/LiCSAlert_results.pkl b/LiCSAlert_03_Campi_Flegrei/LiCSAlert_results.pkl index 323b42e7..5d421060 100644 Binary files a/LiCSAlert_03_Campi_Flegrei/LiCSAlert_results.pkl and b/LiCSAlert_03_Campi_Flegrei/LiCSAlert_results.pkl differ diff --git a/LiCSAlert_batch_mode_examples.py b/LiCSAlert_batch_mode_examples.py index 9219bf8c..7a9b2c66 100755 --- 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 index 4394eb0a..5fcd2000 100644 Binary files a/licsalert/__pycache__/aux.cpython-37.pyc and b/licsalert/__pycache__/aux.cpython-37.pyc differ diff --git a/licsalert/__pycache__/licsalert.cpython-37.pyc b/licsalert/__pycache__/licsalert.cpython-37.pyc index 7f7cf44d..735e3ad9 100644 Binary files a/licsalert/__pycache__/licsalert.cpython-37.pyc and b/licsalert/__pycache__/licsalert.cpython-37.pyc differ diff --git a/licsalert/__pycache__/monitoring_functions.cpython-37.pyc b/licsalert/__pycache__/monitoring_functions.cpython-37.pyc index 4cb08d43..f19a24e9 100644 Binary files a/licsalert/__pycache__/monitoring_functions.cpython-37.pyc and b/licsalert/__pycache__/monitoring_functions.cpython-37.pyc differ diff --git a/licsalert/aux.py b/licsalert/aux.py index ca863569..8a709ccd 100755 --- 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) diff --git a/licsalert/licsalert.py b/licsalert/licsalert.py index 5eccd6d2..1fe22728 100644 --- a/licsalert/licsalert.py +++ b/licsalert/licsalert.py @@ -4,12 +4,13 @@ @author: Matthew Gaddes """ +import pdb #%% def LiCSAlert_batch_mode(displacement_r2, n_baseline_end, out_folder, ICASAR_settings, run_ICASAR = True, ICASAR_path = None, ic_classifying_model = None, - intermediate_figures = False, downsample_run = 1.0, downsample_plot = 0.5, t_recalculate = 10): + intermediate_figures = False, downsample_run = 1.0, downsample_plot = 0.5, t_recalculate = 10, residual_type = 'cumulative'): """ A function to run the LiCSAlert algorithm on a preprocssed time series. To run on a time series that is being updated, use LiCSAlert_monitoring_mode. @@ -33,6 +34,8 @@ def LiCSAlert_batch_mode(displacement_r2, n_baseline_end, out_folder, downsample_run | float | data can be downsampled to speed things up downsample_plot | float | and a 2nd time for fast plotting. Note this is applied to the restuls of the first downsampling, so is compound t_recalculate | int | rolling lines of best fit are recalcaluted every X times (nb done in number of data points, not time between them) + residual_type | str | 'cumulative' or 'window'. If cumulative, residual is the mean of the cumulative (ie summer in time for each pixel, then averaged in space across each ifg. ) + If window, the residual is the ratio of the max in a window (e.g. 20x20 pixels) over the mean for all windows for that time. Value at each point is the cumulative. Returns: out_folder with various items. History: @@ -113,9 +116,9 @@ def LiCSAlert_batch_mode(displacement_r2, n_baseline_end, out_folder, displacement_r2 = LiCSAlert_preprocessing(displacement_r2, downsample_run, downsample_plot) # mean centre and downsize the data if run_ICASAR: # set to True if we need to run ICASAR - baseline_data = {'mixtures_r2' : displacement_r2['incremental'][:n_baseline_end], # prepare a dictionary of data for ICASAR + baseline_data = {'ifgs_dc' : displacement_r2['incremental'][:n_baseline_end], # prepare a dictionary of data for ICASAR 'mask' : displacement_r2['mask'], - 'ifg_dates' : displacement_r2['ifg_dates'][:n_baseline_end]} # ifg dates are stored in displacement_r2, + 'ifg_dates_dc' : displacement_r2['ifg_dates'][:n_baseline_end]} # ifg dates are stored in displacement_r2, for optional_data in ['dem', 'lons', 'lats']: # these are not guaranteed to be supplied to LiCSAlert, to check if they have been, loop through each one. if optional_data in displacement_r2: # if it's in the main displacement_r2 dict @@ -168,8 +171,8 @@ def LiCSAlert_batch_mode(displacement_r2, n_baseline_end, out_folder, # 5b: Or just do LiCSAlet and the LiCSAlert figure for the final time step (much quicker, but only one LiCSAlert figure is created) else: sources_tcs_monitor, residual_monitor = LiCSAlert(sources, tbaseline_info['baselines_cumulative'], displacement_r2["incremental"][:n_baseline_end], # Run LiCSAlert once, on the whole time series. - displacement_r2["incremental"][n_baseline_end:], t_recalculate=t_recalculate, - out_file = out_folder / 'LiCSAlert_results.pkl') + displacement_r2['mask'], displacement_r2["incremental"][n_baseline_end:], t_recalculate=t_recalculate, + out_file = out_folder / 'LiCSAlert_results.pkl', residual_type = residual_type) LiCSAlert_figure(sources_tcs_monitor, residual_monitor, sources, displacement_r2, n_baseline_end, # and only make the plot once tbaseline_info['baselines_cumulative'], time_value_end=tbaseline_info['baselines_cumulative'][-1], day0_date = tbaseline_info['acq_dates'][0], @@ -178,17 +181,22 @@ def LiCSAlert_batch_mode(displacement_r2, n_baseline_end, out_folder, #%% -def LiCSAlert(sources, time_values, ifgs_baseline, ifgs_monitoring = None, t_recalculate = 10, verbose=False, out_file=None): +def LiCSAlert(sources, time_values, ifgs_baseline, mask, ifgs_monitoring = None, t_recalculate = 10, verbose=False, out_file=None, + n_pix_window = 20, residual_type = 'cumulative'): """ Main LiCSAlert algorithm for a daisy-chain timeseries of interferograms. Inputs: sources | r2 array | sources (from ICASAR) as row vectors, as per ICA, that can be turned back to interferograms with a rank 2 boolean mask of which pixels are masked and the col_to_ma function. time_values | r1 array | time values for each point in the time series, commonly (12,24,36) for Sentinel-1 data. Could also be described as the cumulative temporal baselines. ifgs_baseline | r2 array | ifgs used in baseline stage as row vectors + mask | r2 array boolean | mask to convert a row vector (e.g. a row in sources or ifgs) back to a rank2 array. ifgs_monitoring | r2 array or None | ifgs used in monitoring stage as row vectors t_recalculate | int | rolling lines of best fit are recalcaluted every X times (nb done in number of data points, not time between them) verbose | boolean | if True, various information is printed to screen. out_file | None or string | If not None, name of file that is used to save the LiCSALert outputs. + n_pix_window | int | side length of square windows used for the window residual calculation. + residual_type | str | 'cumulative' or 'window'. If cumulative, residual is the mean of the cumulative (ie summer in time for each pixel, then averaged in space across each ifg. ) + If window, the residual is the ratio of the max in a window (e.g. 20x20 pixels) over the mean for all windows for that time. Value at each point is the cumulative. Outputs sources_tcs_monitor | list of dicts | list, with item for each time course. Each dictionary contains: @@ -206,6 +214,7 @@ def LiCSAlert(sources, time_values, ifgs_baseline, ifgs_monitoring = None, t_rec 2019/12/XX | MEG | Written from existing script. 2020/02/16 | MEG | Update to work with no monitoring interferograms 2020/12/14 | MEG | Improve docs and add option to save outputs. + 2021_11_30 | MEG | Add new residual method (max residual of spatial window / mean of residual of all spatial windows for each ifg) """ import numpy as np import pickle @@ -247,10 +256,10 @@ def LiCSAlert(sources, time_values, ifgs_baseline, ifgs_monitoring = None, t_rec # 1: calculating time courses/distances etc for the baseline data tcs_c, _ = bss_components_inversion(sources, ifgs_baseline, cumulative=True) # compute cumulative time courses for baseline interferograms (ie simple inversion to fit each ifg in ifgs_baseline using sources. ) - sources_tcs = tcs_baseline(tcs_c, time_values[:n_times_baseline], t_recalculate) # lines, gradients, etc for time courses. Note done by tcs_baseline, which contrasts with tcs_monitoring (which is used lower down) - _, residual_cb = residual_for_pixels(sources, sources_tcs, ifgs_baseline) # get the cumulative residual for the baseline interferograms - residual_tcs = tcs_baseline(residual_cb, time_values[:n_times_baseline], t_recalculate) # lines, gradients. etc for residual - del tcs_c, residual_cb + sources_tcs = tcs_baseline(tcs_c, time_values[:n_times_baseline], t_recalculate) # lines, gradients, etc for time courses. Note done by tcs_baseline, which contrasts with tcs_monitoring (which is used lower down) + residual = residual_for_pixels(sources, sources_tcs, ifgs_baseline, mask, n_pix_window, residual_type) # get the cumulative residual for the baseline interferograms, sources are stored as row vectors, and each has an entry in sources_tcs which is a dict (ie a list of dicts). Ifgs_baseline are ifgs as row vectors, mean centered in space (ie mean of each row is 0) + residual_tcs = tcs_baseline(residual, time_values[:n_times_baseline], t_recalculate) # lines, gradients. etc for residual + del tcs_c, residual #2: Calculate time courses/distances etc for the monitoring data if ifgs_monitoring is not None: @@ -258,8 +267,8 @@ def LiCSAlert(sources, time_values, ifgs_baseline, ifgs_monitoring = None, t_rec sources_tcs_monitor = tcs_monitoring(tcs_c, sources_tcs, time_values) # update lines, gradients, etc for time courses Note done by tcs_monitoring, which contrasts with tcs_baseline (which is used before) #3: and update the residual stuff # which is handled slightly differently as must be recalcualted for baseline and monitoring data - _, residual_c_bm = residual_for_pixels(sources, sources_tcs_monitor, ifgs_all) # get the cumulative residual for baseline and monitoring (hence _cb) - residual_tcs_monitor = tcs_monitoring(residual_c_bm, residual_tcs, time_values, residual=True) # lines, gradients. etc for residual + residual_bm = residual_for_pixels(sources, sources_tcs_monitor, ifgs_all, mask, n_pix_window, residual_type) # get the cumulative residual for baseline and monitoring (hence _cb) + residual_tcs_monitor = tcs_monitoring(residual_bm, residual_tcs, time_values, residual=True) # lines, gradients. etc for residual if ifgs_monitoring is None: @@ -282,7 +291,7 @@ def LiCSAlert(sources, time_values, ifgs_baseline, ifgs_monitoring = None, t_rec #%% -def residual_for_pixels(sources, sources_tcs, ifgs, n_skip=None): +def residual_for_pixels(sources, sources_tcs, ifgs, mask, n_pix_window = 20, residual_type = 'cumulative', n_skip=None): """ Given spatial sources and their time courses, reconstruct the entire time series and calcualte: - RMS of the residual between each reconstructed and real ifg @@ -294,7 +303,11 @@ def residual_for_pixels(sources, sources_tcs, ifgs, n_skip=None): sources | r2 array | sources as row vectors sources_tcs | list of dicts | As per LiCSAlert, a list with an item for each sources, and each item is a dictionary of various itmes ifgs | r2 array | interferograms as row vectors + mask | r2 array of boolean | to convert a row vector back to a rank 2 masked array. n_skip | None or int | if an int, the first n_skip values of the timecourses will be skipped. + n_pix_window | int | side length of square windows used for the window residual calculation. + residual_type | str | 'cumulative' or 'window'. If cumulative, residual is the mean of the cumulative (ie summer in time for each pixel, then averaged in space across each ifg. ) + If window, the residual is the ratio of the max in a window (e.g. 20x20 pixels) over the mean for all windows for that time. Value at each point is the cumulative. Outputs: residual_ts | r2 array | Column vector of the RMS residual between that ifg, and its reconstruction @@ -306,9 +319,16 @@ def residual_for_pixels(sources, sources_tcs, ifgs, n_skip=None): 2019/12/06 | MEG | Comment and documentation 2020/01/02 | MEG | Update to use new LiCSAlert list of dictionaries 2020/02/06 | MEG | Fix bug as had forgotten to convert cumulative time courses to be incremental + 2021_11_30 | MEG | Add new residual method (ratio of maximum in window for mean of all windows) """ + # import warnings + # warnings.filterwarnings("error") + + import numpy as np + import numpy.ma as ma + from licsalert.aux import r2_to_r3 def list_dict_to_r2(sources_tcs): """ Extract the cumulative time courses from the sources_tcs list of dicts, and convert them to incremental @@ -326,19 +346,61 @@ def list_dict_to_r2(sources_tcs): (n_sources, n_pixs) = sources.shape # number of sources and number of pixels - tcs = list_dict_to_r2(sources_tcs) # get the incremental time courses as a rank 2 array + tcs = list_dict_to_r2(sources_tcs) # get the incremental time courses as a rank 2 array (ie. as column vectors, so something like 55x5, if there are 55 times and 5 sources) if n_skip is not None: # crop/remove the first ifgs tcs = tcs[n_skip:,] # usually the baseline ifgs when used with monitoring data - data_model_residual = ifgs - (tcs @ sources) # residual for each pixel at each time - data_model_residual_cs = np.cumsum(data_model_residual, axis = 0) # summing the residual for each pixel cumulatively through time + data_model_residual = ifgs - (tcs @ sources) # residual for each pixel at each time (ie the data - the reconstrutcion) + data_model_residual_cs = np.cumsum(data_model_residual, axis = 0) # summing the residual for each pixel cumulatively through time, same size as above. residual_ts = np.zeros((data_model_residual.shape[0], 1)) # initiate, n_ifgs x 1 array residual_cs = np.zeros((data_model_residual.shape[0], 1)) # initiate, n_ifgs x 1 array for row_n in range(data_model_residual.shape[0]): # loop through each ifg residual_ts[row_n, 0] = np.sqrt(np.sum(data_model_residual[row_n,:]**2)/n_pixs) # RMS of residual for each ifg residual_cs[row_n, 0] = np.sqrt(np.sum(data_model_residual_cs[row_n,:]**2)/n_pixs) # RMS of residual for cumulative ifgs. i.e. if atmosphere reverses, should cancel in 2nd cumulative ifg. + + # method 2: calculate the residual in spatial windows. + residuals_cs_r3 = r2_to_r3(data_model_residual_cs, mask) # convert from row vectors to rank 2 images. n_times x ny x nx + (n_residuals, ny, nx) = residuals_cs_r3.shape + n_windows_x = int(np.floor(nx / n_pix_window)) # the number of whole windows we can it in the x direction. + rem_x = nx % (n_pix_window * n_windows_x) # the number of pixels remaining if we do that. + border_x = int(rem_x / 2) # what half the border will be + n_windows_y = int(np.floor(ny / n_pix_window)) # as above but for y + rem_y = ny % (n_pix_window * n_windows_y) + border_y = int(rem_y / 2) - return residual_ts, residual_cs + residual_cs_window = np.zeros((data_model_residual.shape[0], 1)) # initiate, n_ifgs x 1 array, to store the ratio of the max residual to the mean residual for each ifg. + residual_cs_windows_r3 = np.zeros((n_residuals, n_windows_y, n_windows_x)) # initiate, n_ifgs x 1 array + for n_resid, residual_cs_r3 in enumerate(residuals_cs_r3): + for row_n in range(n_windows_y): + for col_n in range(n_windows_x): + window_region = residuals_cs_r3[n_resid, (border_y + row_n*n_pix_window):(border_y + (row_n+1)*n_pix_window), + (border_x + col_n*n_pix_window):(border_x + (col_n+1)*n_pix_window) ] # extract the window region + n_coherent_pixs = (n_pix_window**2) - np.sum(window_region.mask) # get the number of coherent pixels in the windo (maksed pixels are 1) + if False not in window_region.mask: # if there are no unmasked pixels (unmasked = False) + residual_cs_windows_r3[n_resid, row_n, col_n] = np.nan # then the residual for that window is just nan + else: + residual_cs_windows_r3[n_resid, row_n, col_n] = np.sqrt(ma.sum(window_region**2)/n_coherent_pixs) # ge the RMS error for the window, and assign to the array that stores them. + + residual_cs_window[n_resid] = np.nanmax(residual_cs_windows_r3[n_resid,]) / np.nanmean(residual_cs_windows_r3[n_resid,]) # calculate the ratio between the largest and mean residual + + # debugging/ examining plot + # import numpy.ma as ma + # import matplotlib.pyplot as plt + # n_plot = 15 + # f, axes = plt.subplots(2,n_plot, figsize = (30,7)) + # residuals = residuals_cs_r3[:n_plot,] + # window_residuals = residual_cs_windows_r3[:n_plot,] + # for i in range(n_plot): + # #axes[0,i].imshow(residuals[i,], vmin = ma.min(residuals), vmax = ma.max(residuals)) # share colorbar + # axes[0,i].imshow(residuals[i,]) # dont + # axes[1,i].imshow(window_residuals[i,], vmin = np.nanmin(window_residuals), vmax = np.nanmax(window_residuals)) # share colorbar + if residual_type == 'cumulative': + return residual_cs + elif residual_type == 'window': + return residual_cs_window + else: + raise Exception(f"'residual_type' must be either 'cumulative' or 'window'. Exiting as it's {residual_type}") + #%% @@ -693,7 +755,7 @@ def xticks_every_3months(ax_to_update, day0_date, time_values, include_tick_labe figtitle = f'LiCSAlert figure with {(n_ifgs-n_baseline_end):03d} monitoring interferograms' # 2 Initiate the figure - fig1 = plt.figure(figsize=(14,8)) + fig1 = plt.figure(figsize=(18,10)) fig1.canvas.manager.set_window_title(figtitle) grid = gridspec.GridSpec((n_ics + 3), 11, wspace=0.3, hspace=0.1) # divide into 2 sections, 1/5 for ifgs and 4/5 for components diff --git a/licsalert/monitoring_functions.py b/licsalert/monitoring_functions.py index b3fa8f2b..f2732458 100644 --- a/licsalert/monitoring_functions.py +++ b/licsalert/monitoring_functions.py @@ -5,7 +5,7 @@ @author: matthew """ - +import pdb #%% @@ -559,7 +559,7 @@ def LiCSAlert_dates_status(LiCSAlert_required_dates, LiCSAlert_dates, folder_LiC #%% def run_LiCSAlert_status(licsbas_dates, volcano_path, date_baseline_end): - """ When 'LiCSAlert_monitoring_mode' is run in a folder, this function determines which steps + """ When 'LiCSAlert_monitoring_mode' is run in a directory, this function determines which steps need to be done, ranging from nothing (everything is up to date), through to runnings LiCSBAS, ICASAR, and LiCSAlert. Inputs: