Add get_pairs_matrix() function and use it to define least-squares an…

…d unwrapping matrices. Increase library version number.

pygmtsar/pygmtsar/Stack_base.py

@@ -251,3 +251,37 @@ def get_pairs(self, pairs, dates=False):
             dates = np.unique(pairs[['ref', 'rep']].astype(str).values.flatten())
             return (pairs, dates)
         return pairs
+
+    def get_pairs_matrix(self, pairs):
+        """
+        Create a matrix based on interferogram dates and pairs.
+
+        Parameters
+        ----------
+        pairs : pandas.DataFrame or xarray.DataArray or xarray.Dataset
+            DataFrame or DataArray containing interferogram date pairs.
+        
+        Returns
+        -------
+        numpy.ndarray
+            A matrix with one row for every interferogram and one column for every date.
+            Each element in the matrix is a float, with 1 indicating the start date,
+            -1 indicating the end date, 0 if the date is covered by the corresponding 
+            interferogram timeline, and NaN otherwise.
+
+        """
+        import numpy as np
+        import pandas as pd
+
+        # also define image capture dates from interferogram date pairs
+        pairs, dates = self.get_pairs(pairs, dates=True)
+        pairs = pairs[['ref', 'rep']].astype(str).values
+
+        # here are one row for every interferogram and one column for every date
+        matrix = []
+        for pair in pairs:
+            #mrow = [date>=pair[0] and date<=pair[1] for date in dates]
+            mrow = [(-1 if date==pair[0] else (1 if date==pair[1] else (0 if date>pair[0] and date<pair[1] else np.nan))) for date in dates]
+            matrix.append(mrow)
+        matrix = np.stack(matrix).astype(np.float32)
+        return matrix

pygmtsar/pygmtsar/Stack_lstsq.py

@@ -90,7 +90,7 @@ def lstsq_matrix(self, pairs):
 
         Parameters
         ----------
-        pairs : pandas.DataFrame
+        pairs : pandas.DataFrame or Xarray object
             DataFrame containing interferogram date pairs.
 
         Returns
@@ -100,34 +100,16 @@ def lstsq_matrix(self, pairs):
             Each element in the matrix is an integer, with 1 representing that the date
             is between the corresponding interferogram's reference and repeat dates, and
             0 otherwise.
-
-        Notes
-        -----
-        This function also calculates image capture dates from the interferogram date pairs.
-
         """
-        import numpy as np
-        import pandas as pd
-
-        # also define image capture dates from interferogram date pairs
-        pairs, dates = self.get_pairs(pairs, dates=True)
-        pairs = pairs[['ref', 'rep']].astype(str).values
-
-        # here are one row for every interferogram and one column for every date
-        matrix = []
-        for pair in pairs:
-            mrow = [date>pair[0] and date<=pair[1] for date in dates]
-            matrix.append(mrow)
-        matrix = np.stack(matrix).astype(int)
-        return matrix
+        return (self.get_pairs_matrix(pairs)>=0).astype(int)
 
     def lstsq_matrix_edge(self, pairs):
         """
         Create an edge matrix for use in the least squares computation based on interferogram date pairs.
 
         Parameters
         ----------
-        pairs : pandas.DataFrame
+        pairs : pandas.DataFrame or Xarray object
             DataFrame containing interferogram date pairs.
 
         Returns
@@ -136,33 +118,10 @@ def lstsq_matrix_edge(self, pairs):
             Matrix with one row for every interferogram and one column for every date.
             Each element in the matrix is an integer, with 1 representing the end date,
             -1 the start date and 0 otherwise.
-
-        Notes
-        -----
-        This function also calculates image capture dates from the interferogram date pairs.
-
         """
+        print ('NOTE: this function is not used in the code and created for test purposes only.')
         import numpy as np
-        import pandas as pd
-
-        # also define image capture dates from interferogram date pairs
-        pairs, dates = self.get_pairs(pairs, dates=True)
-        pairs = pairs[['ref', 'rep']].astype(str).values
-
-        def edge(date, date1, date2):
-            if date == date1:
-                return -1
-            if date == date2:
-                return 1
-            return 0
-
-        # here are one row for every interferogram and one column for every date
-        matrix = []
-        for pair in pairs:
-            mrow = [edge(date, pair[0], pair[1]) for date in dates]
-            matrix.append(mrow)
-        matrix = np.stack(matrix).astype(int)
-        return matrix
+        return np.nan_to_num(self.get_pairs_matrix(pairs)).astype(int)
 
     def lstsq(self, data, weight=None, matrix='auto', cumsum=True, debug=False):
         """

pygmtsar/pygmtsar/Stack_sbas.py

@@ -495,17 +495,20 @@ def plot_baseline_displacement(self, phase, corr=None, caption='LSQ and STL Disp
             print ("NOTE: Displacement is automatically set to 'True' because it is required for 'stl=True'.")
         assert isinstance(phase, xr.DataArray) and phase.dims == ('pair',), \
             'ERROR: Argument phase should be 1D Xarray with "pair" dimension'
+        assert phase.name is not None, \
+            'ERROR: Argument phase should be named Xarray array'
         plt.figure()
         colors = plt.colormaps[cmap]
 
         df = phase.to_dataframe()
         df['corr'] = corr.values if corr is not None else 1
         pairs, dates = self.get_pairs(phase, dates=True)
         dates = pd.DatetimeIndex(dates)
-        matrix = self.lstsq_matrix(pairs)
+        lstsq_matrix = self.lstsq_matrix(pairs)
+        unwrap_matrix = self.unwrap_matrix(pairs)
 
         if unwrap:
-            df['phase'] = self.unwrap_pairs(phase.values, df['corr'].values, matrix, tolerance)
+            df['phase'] = self.unwrap_pairs(phase.values, df['corr'].values, unwrap_matrix, tolerance)
 
         unit = 'rad'
         name = 'Phase'
@@ -515,7 +518,7 @@ def plot_baseline_displacement(self, phase, corr=None, caption='LSQ and STL Disp
             name = 'Displacement'
 
         if displacement or stl:
-            solution = self.lstsq1d(df['phase'].values, 0.999*df['corr'].values if corr is not None else None, matrix)
+            solution = self.lstsq1d(df['phase'].values, 0.999*df['corr'].values if corr is not None else None, lstsq_matrix)
             #print ('solution', solution)
             days = (dates - dates[0]).days
             #print ('days', days)

pygmtsar/pygmtsar/Stack_unwrap.py

@@ -382,12 +382,31 @@ def unwrap_pairs(data      : np.ndarray,
         #out[~nanmask] = np.where(mask[~nanmask], buffer[~nanmask], np.nan)
         return out
 
+    def unwrap_matrix(self, pairs):
+        """
+        Create a matrix for use in the least squares computation based on interferogram date pairs.
+
+        Parameters
+        ----------
+        pairs : pandas.DataFrame or xarray.DataArray or xarray.Dataset
+            DataFrame or DataArray containing interferogram date pairs.
+
+        Returns
+        -------
+        numpy.ndarray
+            A matrix with one row for every interferogram and one column for every date.
+            Each element in the matrix is a float, with 1 indicating the start date,
+            -1 indicating the end date, 0 if the date is covered by the corresponding
+            interferogram timeline, and NaN otherwise.
+        """
+        return self.get_pairs_matrix(pairs)
+
     def unwrap1d(self, data, weight=None, tolerance=np.pi/2):
         import xarray as xr
         import numpy as np
 
         pairs = self.get_pairs(data)
-        matrix = self.lstsq_matrix(pairs)
+        matrix = self.unwrap_matrix(pairs)
 
         if not 'stack' in data.dims:
             chunks_z, chunks_y, chunks_x = data.chunks if data.chunks is not None else np.inf, np.inf, np.inf

pygmtsar/pygmtsar/__init__.py

@@ -7,7 +7,7 @@
 # 
 # Licensed under the BSD 3-Clause License (see LICENSE for details)
 # ----------------------------------------------------------------------------
-__version__ = '2024.7.14'
+__version__ = '2024.8.2'
 
 # unified progress indicators
 from .tqdm_joblib import tqdm_joblib