Fix stack processing in regression1d() and elevation_m() functions

pygmtsar/pygmtsar/Stack_detrend.py

@@ -1001,6 +1001,10 @@ def regression1d(self, data, dim='auto', degree=1):
         multi_index = None
         if 'stack' in data.dims and isinstance(data.coords['stack'].to_index(), pd.MultiIndex):
             multi_index = data['stack']
+            # detect unused coordinates
+            unused_coords = [d for d in multi_index.coords if not d in multi_index.dims and not d in multi_index.indexes]
+            # cleanup multiindex to merge it with the processed dataset later
+            multi_index = multi_index.drop_vars(unused_coords)
             data = data.reset_index('stack')
 
         stackdim = [_dim for _dim in ['date', 'pair'] if _dim in data.dims]
@@ -1018,10 +1022,14 @@ def regression1d(self, data, dim='auto', degree=1):
             else:
                 dim_da = xr.DataArray(dim, dims=[stackdim])
             data_dim = data.assign_coords(polyfit_coord=dim_da).swap_dims({'pair': 'polyfit_coord'})
+            # Polynomial coefficients, highest power first, see numpy.polyfit
             fit_coeff = data_dim.polyfit('polyfit_coord', degree).polyfit_coefficients.astype(np.float32)
             fit = xr.polyval(data_dim['polyfit_coord'], fit_coeff)\
                 .swap_dims({'polyfit_coord': stackdim}).drop_vars('polyfit_coord').astype(np.float32).rename('trend')
-            return xr.merge([fit, fit_coeff]).rename(polyfit_coefficients='coefficients')
+            out = xr.merge([fit, fit_coeff]).rename(polyfit_coefficients='coefficients')
+            if multi_index is not None:
+                return out.assign_coords(stack=multi_index)
+            return out
 
         # fit existing coordinate values
         fit_coeff = data.polyfit(dim, degree).polyfit_coefficients.astype(np.float32)

pygmtsar/pygmtsar/Stack_incidence.py

@@ -439,11 +439,20 @@ def elevation_m(self, data, baseline=1):
         # expected accuracy about 0.01%
         #wavelength, slant_range = self.PRM().get('radar_wavelength','SC_height')
         wavelength, slant_range_start,slant_range_end = self.PRM().get('radar_wavelength', 'SC_height_start', 'SC_height_end')
-        slant_range = xr.DataArray(np.linspace(slant_range_start,slant_range_end, data.shape[1]),
-                                   coords=[data.coords['x']])
-        incidence = self.incidence_angle().reindex_like(data, method='nearest')
+
+        incidence_angle = self.incidence_angle()
+        slant_range = xr.DataArray(np.linspace(slant_range_start,slant_range_end, incidence_angle.shape[1]),
+                                       coords=[incidence_angle.coords['x']])
+
+        if 'stack' in data.dims and 'y' in data.coords and 'x' in data.coords:
+            incidence = incidence_angle.interp(y=data.y, x=data.x, method='linear')
+            slant = slant_range.interp(x=data.x, method='linear')
+        else:
+            incidence = incidence_angle.reindex_like(data, method='nearest')
+            slant = slant_range.reindex_like(data.x, method='nearest')
+
         # sign corresponding to baseline and phase signs
-        return -(wavelength*data*slant_range*np.cos(incidence)/(4*np.pi*baseline)).rename('ele')
+        return -(wavelength*data*slant*np.cos(incidence)/(4*np.pi*baseline)).rename('ele')
 
     def compute_satellite_look_vector(self, interactive=False):
         #import dask