Unify data arguments name

pygmtsar/pygmtsar/Stack_detrend.py

@@ -323,6 +323,11 @@ def polyfit(self, data, weight=None, degree=0, days=None, count=None):
         import xarray as xr
         import pandas as pd
         import numpy as np
+        import warnings
+        # suppress Dask warning "RuntimeWarning: invalid value encountered in divide"
+        warnings.filterwarnings('ignore')
+        warnings.filterwarnings('ignore', module='dask')
+        warnings.filterwarnings('ignore', module='dask.core')
 
         pairs, dates = self.get_pairs(data, dates=True)
 
@@ -406,13 +411,13 @@ def polyfit(self, data, weight=None, degree=0, days=None, count=None):
 # 
 #         return out.rename(data.name)
 
-    def gaussian(self, grid, wavelength, truncate=3.0, resolution=60, debug=False):
+    def gaussian(self, data, wavelength, truncate=3.0, resolution=60, debug=False):
         """
         Apply a lazy Gaussian filter to an input 2D or 3D data array.
 
         Parameters
         ----------
-        dataarray : xarray.DataArray
+        data : xarray.DataArray
             The input data array with NaN values allowed.
         wavelength : float
             The cut-off wavelength for the Gaussian filter in meters.
@@ -454,12 +459,12 @@ def gaussian(self, grid, wavelength, truncate=3.0, resolution=60, debug=False):
 #         warnings.filterwarnings('ignore', module='dask')
 #         warnings.filterwarnings('ignore', module='dask.core')
 
-        assert self.is_ra(grid), 'ERROR: the processing requires grid in radar coordinates'
-        assert np.issubdtype(grid.dtype, np.floating), 'ERROR: expected float datatype input grid'
+        assert self.is_ra(data), 'ERROR: the processing requires grid in radar coordinates'
+        assert np.issubdtype(data.dtype, np.floating), 'ERROR: expected float datatype input data'
         assert wavelength is not None, 'ERROR: Gaussian filter cut-off wavelength is not defined'
 
         # ground pixel size
-        dy, dx = self.get_spacing(grid)
+        dy, dx = self.get_spacing(data)
         # downscaling
         yscale, xscale = int(np.round(resolution/dy)), int(np.round(resolution/dx))
         # gaussian kernel
@@ -469,7 +474,7 @@ def gaussian(self, grid, wavelength, truncate=3.0, resolution=60, debug=False):
             print (f'DEBUG: gaussian: ground pixel size in meters: y={dy:.1f}, x={dx:.1f}')
         if (xscale <=1 and yscale <=1) or (wavelength/resolution <= 3):
             # decimation is useless
-            return self.multilooking(grid, wavelength=wavelength, coarsen=None, debug=debug)
+            return self.multilooking(data, wavelength=wavelength, coarsen=None, debug=debug)
 
         # define filter on decimated grid, the correction value is typically small
         wavelength_dec = np.sqrt(wavelength**2 - resolution**2)
@@ -479,34 +484,34 @@ def gaussian(self, grid, wavelength, truncate=3.0, resolution=60, debug=False):
             print (f'DEBUG: gaussian: filtering on {resolution}m grid using wavelength {wavelength_dec:.1f}')
 
         # find stack dim
-        stackvar = grid.dims[0] if len(grid.dims) == 3 else 'stack'
+        stackvar = data.dims[0] if len(data.dims) == 3 else 'stack'
         #print ('stackvar', stackvar)
 
         # split coordinate grid to equal chunks and rest
-        ys_blocks = np.array_split(grid.y, np.arange(0, grid.y.size, self.chunksize)[1:])
-        xs_blocks = np.array_split(grid.x, np.arange(0, grid.x.size, self.chunksize)[1:])
+        ys_blocks = np.array_split(data.y, np.arange(0, data.y.size, self.chunksize)[1:])
+        xs_blocks = np.array_split(data.x, np.arange(0, data.x.size, self.chunksize)[1:])
 
-        grid_dec = self.multilooking(grid, wavelength=resolution, coarsen=(yscale,xscale), debug=debug)
-        grid_dec_gauss = self.multilooking(grid_dec, wavelength=wavelength_dec, debug=debug)
-        del grid_dec
+        data_dec = self.multilooking(data, wavelength=resolution, coarsen=(yscale,xscale), debug=debug)
+        data_dec_gauss = self.multilooking(data_dec, wavelength=wavelength_dec, debug=debug)
+        del data_dec
 
         stack = []
-        for stackval in grid[stackvar].values if len(grid.dims) == 3 else [None]:
-            grid_in = grid_dec_gauss.sel({stackvar: stackval}) if stackval is not None else grid_dec_gauss
-            grid_out = grid_in.reindex({'y': grid.y, 'x': grid.x}, method='nearest').chunk(self.chunksize)
-            del grid_in
-            stack.append(grid_out)
-            del grid_out
+        for stackval in data[stackvar].values if len(data.dims) == 3 else [None]:
+            data_in = data_dec_gauss.sel({stackvar: stackval}) if stackval is not None else data_dec_gauss
+            data_out = data_in.reindex({'y': data.y, 'x': data.x}, method='nearest').chunk(self.chunksize)
+            del data_in
+            stack.append(data_out)
+            del data_out
 
         # wrap lazy Dask array to Xarray dataarray
-        if len(grid.dims) == 2:
+        if len(data.dims) == 2:
             out = stack[0]
         else:
             out = xr.concat(stack, dim=stackvar)
         del stack
 
-        # append source grid coordinates excluding removed y, x ones
-        for (k,v) in grid.coords.items():
+        # append source data coordinates excluding removed y, x ones
+        for (k,v) in data.coords.items():
             if k not in ['y','x']:
                 out[k] = v
 

pygmtsar/pygmtsar/Stack_incidence.py

@@ -349,13 +349,13 @@ def plot_incidence_angle(self, incidence_angle='auto', caption='Incidence Angle
         plt.ylabel('Azimuth')
         plt.title(caption)
 
-    def vertical_displacement_mm(self, unwrap):
+    def vertical_displacement_mm(self, data):
         """
         Compute vertical displacement in millimeters in radar coordinates.
 
         Parameters
         ----------
-        unwrap : xarray.DataArray or xarray.Dataset
+        data : xarray.DataArray or xarray.Dataset
             Unwrapped phase grid(s) in radar coordinates.
 
         Returns
@@ -369,19 +369,19 @@ def vertical_displacement_mm(self, unwrap):
         """
         import numpy as np
 
-        assert self.is_ra(unwrap), 'ERROR: unwrapped phase needs to be defined in radar coordinates'
+        assert self.is_ra(data), 'ERROR: argument data needs to be defined in radar coordinates'
 
-        los_disp = self.los_displacement_mm(unwrap)
+        los_disp = self.los_displacement_mm(data)
         incidence = self.incidence_angle().reindex_like(los_disp, method='nearest')
         return los_disp/np.cos(incidence)
 
-    def eastwest_displacement_mm(self, unwrap):
+    def eastwest_displacement_mm(self, data):
         """
         Compute East-West displacement in millimeters.
 
         Parameters
         ----------
-        unwraps : xarray.DataArray or xarray.Dataset
+        data : xarray.DataArray or xarray.Dataset
             Unwrapped phase grid(s) in geographic coordinates.
 
         Returns
@@ -404,17 +404,17 @@ def eastwest_displacement_mm(self, unwrap):
         # this displacement is not symmetrical for the orbits due to scene geometries
         orbit = self.df.orbit.unique()[0]
         sign = 1 if orbit == 'D' else -1
-        los_disp = self.los_displacement_mm(unwrap)
+        los_disp = self.los_displacement_mm(data)
         incidence_ll = self.incidence_angle().reindex_like(los_disp, method='nearest')
         return sign * los_disp/np.sin(incidence_ll)
 
-    def elevation_m(self, unwrap, baseline):
+    def elevation_m(self, data, baseline):
         """
         Compute elevation in meters in radar coordinates.
 
         Parameters
         ----------
-        unwrap : xarray.DataArray or xarray.Dataset
+        data : xarray.DataArray or xarray.Dataset
             Unwrapped phase grid(s) in radar coordinates.
 
         Returns
@@ -429,16 +429,16 @@ def elevation_m(self, unwrap, baseline):
         import xarray as xr
         import numpy as np
 
-        assert self.is_ra(unwrap), 'ERROR: unwrapped phase needs to be defined in radar coordinates'
+        assert self.is_ra(data), 'ERROR: unwrapped phase needs to be defined in radar coordinates'
 
         # 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, unwrap.shape[1]),
-                                   coords=[unwrap.coords['x']])
-        incidence = self.incidence_angle().reindex_like(unwrap, method='nearest')
+        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')
         # sign corresponding to baseline and phase signs
-        return -(wavelength*unwrap*slant_range*np.cos(incidence)/(4*np.pi*baseline)).rename('ele')
+        return -(wavelength*data*slant_range*np.cos(incidence)/(4*np.pi*baseline)).rename('ele')
 
     def compute_satellite_look_vector(self, interactive=False):
         #import dask