From 4d92f0c2c602cd3f3c3371572f9cd26bb1135d53 Mon Sep 17 00:00:00 2001
From: Alexey Pechnikov <pechnikov@mobigroup.ru>
Date: Sat, 21 Sep 2024 00:16:48 +0700
Subject: [PATCH] Apply on-the-fly subswaths aligning and stacking (incomplete)

---
 pygmtsar/pygmtsar/IO.py          | 192 ++++++++++++++++++--
 pygmtsar/pygmtsar/Stack_align.py | 302 -------------------------------
 pygmtsar/pygmtsar/Stack_trans.py |   9 +-
 3 files changed, 185 insertions(+), 318 deletions(-)

diff --git a/pygmtsar/pygmtsar/IO.py b/pygmtsar/pygmtsar/IO.py
index b23d4d3..78ee1f8 100644
--- a/pygmtsar/pygmtsar/IO.py
+++ b/pygmtsar/pygmtsar/IO.py
@@ -248,38 +248,204 @@ def get_filenames(self, pairs, name, add_subswath=False):
                 filenames.append(filename)
         return filenames
 
+#     # 2.5e-07 is Sentinel-1 scale factor
+#     def open_data(self, dates=None, scale=2.5e-07, debug=False):
+#         import xarray as xr
+#         import pandas as pd
+#         import numpy as np
+#         import os
+# 
+#         if debug:
+#             print ('DEBUG: open_data: apply scale:', scale)
+# 
+#         if dates is None:
+#             dates = self.df.index.values
+#         #print ('dates', dates)
+# 
+#         filenames = [self.PRM(date).filename[:-4] + '.grd' for date in dates]
+#         #print ('filenames', filenames)
+#         ds = xr.open_mfdataset(
+#             filenames,
+#             engine=self.netcdf_engine,
+#             chunks=self.chunksize,
+#             parallel=True,
+#             concat_dim='date',
+#             combine='nested'
+#         ).assign(date=pd.to_datetime(dates)).rename({'a': 'y', 'r': 'x'})
+#         if scale is None:
+#             # there is no complex int16 datatype, so return two variables for real and imag parts
+#             return ds
+#         # scale and return as complex values
+#         ds_scaled = (scale*(ds.re.astype(np.float32) + 1j*ds.im.astype(np.float32))).assign_attrs(ds.attrs).rename('data')
+#         del ds
+#         # zero in np.int16 type means NODATA
+#         return ds_scaled.where(ds_scaled != 0)
+
     # 2.5e-07 is Sentinel-1 scale factor
     def open_data(self, dates=None, scale=2.5e-07, debug=False):
         import xarray as xr
         import pandas as pd
         import numpy as np
+        from scipy import constants
         import os
 
         if debug:
             print ('DEBUG: open_data: apply scale:', scale)
 
         if dates is None:
-            dates = self.df.index.values
+            #dates = self.df.index.values
+            dates = np.unique(self.df.index.values)
         #print ('dates', dates)
 
-        filenames = [self.PRM(date).filename[:-4] + '.grd' for date in dates]
-        #print ('filenames', filenames)
-        ds = xr.open_mfdataset(
-            filenames,
-            engine=self.netcdf_engine,
-            chunks=self.chunksize,
-            parallel=True,
-            concat_dim='date',
-            combine='nested'
-        ).assign(date=pd.to_datetime(dates)).rename({'a': 'y', 'r': 'x'})
+        subswaths = self.get_subswaths()
+        if not isinstance(subswaths, (str, int)):
+            subswaths = ''.join(map(str, subswaths))
+
+        # DEM extent in radar coordinates, merged reference PRM required
+        #print ('minx, miny, maxx, maxy', minx, miny, maxx, maxy)
+        extent_ra = np.round(self.get_extent_ra(subswath=subswaths[0]).bounds).astype(int)
+
+        if len(subswaths) == 1:
+            # stack single subswath
+            stack = []
+            shape = None
+            for date in dates:
+                prm = self.PRM(date, subswath=int(subswaths))
+                # disable scaling
+                slc = prm.read_SLC_int(scale=scale, shape=shape)
+                stack.append(slc.assign_coords(date=date))
+                if shape is None:
+                    shape = (slc.y.size, slc.x.size)
+                del slc, prm
+        else:
+            # calculate the offsets to merge subswaths
+            prms = []
+            ylims = []
+            xlims = []
+            for subswath in subswaths:
+                #print (subswath)
+                prm = self.PRM(subswath=subswath)
+                prms.append(prm)
+                ylims.append(prm.get('num_valid_az'))
+                xlims.append(prm.get('num_rng_bins'))
+    
+            assert len(np.unique([prm.get('PRF') for prm in prms])), 'Image PRFs are not consistent'
+            assert len(np.unique([prm.get('rng_samp_rate') for prm in prms])), 'Image range sampling rates are not consistent'
+    
+            bottoms = [0] + [int(np.round(((prm.get('clock_start') - prms[0].get('clock_start')) * 86400 * prms[0].get('PRF')))) for prm in prms[1:]]
+            # head123: 0, 466, -408
+            if debug:
+                print ('bottoms init', bottoms)
+            # minh: -408
+            minh = min(bottoms)
+            if debug:
+                print ('minh', minh)
+            #head123: 408, 874, 0
+            bottoms = np.asarray(bottoms) - minh
+            if debug:
+                print ('bottoms', bottoms)
+    
+            #ovl12,23: 2690, 2558
+            ovls = [prm1.get('num_rng_bins') - \
+                    int(np.round(((prm2.get('near_range') - prm1.get('near_range')) / (constants.speed_of_light/ prm1.get('rng_samp_rate') / 2)))) \
+                    for (prm1, prm2) in zip(prms[:-1], prms[1:])]
+            if debug:
+                print ('ovls', ovls)
+    
+            #Writing the grid files..Size(69158x13075)...
+            #maxy: 13075
+            # for SLC
+            maxy = max([prm.get('num_valid_az') + bottom for prm, bottom in zip(prms, bottoms)])
+            if debug:
+                print ('maxy', maxy)
+            maxx = sum([prm.get('num_rng_bins') - ovl - 1 for prm, ovl in zip(prms, [-1] + ovls)])
+            if debug:
+                print ('maxx', maxx)
+    
+            #Stitching location n1 = 1045
+            #Stitching location n2 = 935
+            ns = [np.ceil(-prm.get('rshift') + prm.get('first_sample') + 150.0).astype(int) for prm in prms[1:]]
+            ns = [10 if n < 10 else n for n in ns]
+            if debug:
+                print ('ns', ns)
+    
+            # left and right coordinates for every subswath valid area
+            lefts = []
+            rights = []
+    
+            # 1st
+            xlim = prms[0].get('num_rng_bins') - ovls[0] + ns[0]
+            lefts.append(0)
+            rights.append(xlim)
+    
+            # 2nd
+            if len(prms) == 2:
+                xlim = prms[1].get('num_rng_bins') - 1
+            else:
+                # for 3 subswaths
+                xlim = prms[1].get('num_rng_bins') - ovls[1] + ns[1]
+            lefts.append(ns[0])
+            rights.append(xlim)
+    
+            # 3rd
+            if len(prms) == 3:
+                xlim = prms[2].get('num_rng_bins') - 2
+                lefts.append(ns[1])
+                rights.append(xlim)
+    
+            # check and merge SLCs
+            sumx = sum([right-left for right, left in zip(rights, lefts)])
+            if debug:
+                print ('assert maxx == sum(...)', maxx, sumx)
+            assert maxx == sumx, 'Incorrect output grid range dimension size'
+    
+            offsets = {'bottoms': bottoms, 'lefts': lefts, 'rights': rights, 'bottom': minh, 'extent': [maxy, maxx], 'ylims': ylims, 'xlims': xlims}
+            if debug:
+                print ('offsets', offsets)
+    
+            # merge subswaths
+            stack = []
+            for date in dates:
+                slcs = []
+                prms = []
+                
+                for subswath, bottom, left, right, ylim, xlim in zip(subswaths, bottoms, lefts, rights, ylims, xlims):
+                    print (date, subswath)
+                    prm = self.PRM(date, subswath=int(subswath))
+                    # disable scaling
+                    slc = prm.read_SLC_int(scale=scale, shape=(ylim, xlim))
+                    slc = slc.isel(x=slice(left, right)).assign_coords(y=slc.y + bottom)
+                    slcs.append(slc)
+                    prms.append(prm)
+        
+                # check and merge SLCs, use zero fill for np.int16 datatype
+                slc = xr.concat(slcs, dim='x', fill_value=0).assign_coords(x=0.5 + np.arange(maxx))
+        
+                if debug:
+                    print ('assert slc.y.size == maxy', slc.y.size, maxy)
+                assert slc.y.size == maxy, 'Incorrect output grid azimuth dimension size'
+                if debug:
+                    print ('assert slc.x.size == maxx', slc.x.size, maxx)
+                assert slc.x.size == maxx, 'Incorrect output grid range dimension sizes'
+                del slcs
+    
+                stack.append(slc.assign_coords(date=date))
+                del slc
+
+        # minx, miny, maxx, maxy = extent_ra
+        ds = xr.concat(stack, dim='date').assign(date=pd.to_datetime(dates))\
+            .sel(y=slice(extent_ra[1], extent_ra[3]), x=slice(extent_ra[0], extent_ra[2])) \
+            .chunk({'y': self.chunksize, 'x': self.chunksize})
+        del stack
+
         if scale is None:
             # there is no complex int16 datatype, so return two variables for real and imag parts
             return ds
         # scale and return as complex values
-        ds_scaled = (scale*(ds.re.astype(np.float32) + 1j*ds.im.astype(np.float32))).assign_attrs(ds.attrs).rename('data')
+        ds_complex = (ds.re.astype(np.float32) + 1j*ds.im.astype(np.float32))
         del ds
         # zero in np.int16 type means NODATA
-        return ds_scaled.where(ds_scaled != 0)
+        return ds_complex.where(ds_complex != 0).rename('data')
 
 #     def open_geotif(self, dates=None, subswath=None, intensity=False, chunksize=None):
 #         """
diff --git a/pygmtsar/pygmtsar/Stack_align.py b/pygmtsar/pygmtsar/Stack_align.py
index b30ec37..cda4f61 100644
--- a/pygmtsar/pygmtsar/Stack_align.py
+++ b/pygmtsar/pygmtsar/Stack_align.py
@@ -260,13 +260,6 @@ def _align_rep_subswath(self, subswath, date=None, degrees=12.0/3600, debug=Fals
         # Restoring $tmp_da lines shift to the image... 
         PRM.from_file(rep_prm).set(ashift=0 if abs(tmp_da) < 1000 else tmp_da, rshift=0).update()
 
-        # that is safe to rewrite source files
-        prm1 = PRM.from_file(ref_prm)
-        prm1.resamp(PRM.from_file(rep_prm),
-                    repeatSLC_tofile=rep_prm[:-4]+'.SLC',
-                    interp=1, debug=debug
-        ).to_file(rep_prm)
-
         PRM.from_file(rep_prm).set(PRM.fitoffset(3, 3, par_tmp)).update()
 
         PRM.from_file(rep_prm).calc_dop_orb(earth_radius, 0, inplace=True, debug=debug).update()
@@ -276,281 +269,6 @@ def _align_rep_subswath(self, subswath, date=None, degrees=12.0/3600, debug=Fals
             #if os.path.exists(filename):
             os.remove(filename)
 
-    # define bottoms from reference scene and apply for all scenes
-    def get_subswaths_offsets(self, date, offsets=None, debug=False):
-        import xarray as xr
-        import numpy as np
-        from scipy import constants
-
-        subswaths = self.get_subswaths()
-
-        # define offset parameters to merge subswaths
-        prms = []
-        for subswath in subswaths:
-            #print (subswath)
-            prm = self.PRM(date, subswath=subswath)
-            prms.append(prm)
-
-        assert len(np.unique([prm.get('PRF') for prm in prms])), 'Image PRFs are not consistent'
-        assert len(np.unique([prm.get('rng_samp_rate') for prm in prms])), 'Image range sampling rates are not consistent'
-
-        if offsets is None:
-            bottoms = [0] + [((prm.get('clock_start') - prms[0].get('clock_start')) * 86400 * prms[0].get('PRF')).round().astype(int) for prm in prms[1:]]
-            # head123: 0, 466, -408
-            if debug:
-                print ('bottoms init', bottoms)
-            # minh: -408
-            minh = min(bottoms)
-            if debug:
-                print ('minh', minh)
-            #head123: 408, 874, 0
-            bottoms = np.asarray(bottoms) - minh
-        else:
-            bottoms = offsets['bottoms']
-            minh = offsets['bottom']
-        if debug:
-            print ('bottoms', bottoms)
-
-        #ovl12,23: 2690, 2558
-        ovls = [prm1.get('num_rng_bins') - ((prm2.get('near_range') - prm1.get('near_range')) \
-                    / (constants.speed_of_light/ prm1.get('rng_samp_rate') / 2)).round().astype(int) \
-                for (prm1, prm2) in zip(prms[:-1], prms[1:])]
-        if debug:
-            print ('ovls', ovls)
-
-        #Writing the grid files..Size(69158x13075)...
-        #maxy: 13075
-        # for SLC
-        maxy = max([prm.get('num_valid_az') + bottom for prm, bottom in zip(prms, bottoms)])
-        if debug:
-            print ('maxy', maxy)
-        maxx = sum([prm.get('num_rng_bins') - ovl -1 for prm, ovl in zip(prms, [-1] + ovls)])
-        if debug:
-            print ('maxx', maxx)
-
-        #Stitching location n1 = 1045
-        #Stitching location n2 = 935
-        ns = [np.ceil(-prm.get('rshift') + prm.get('first_sample') + 150.0).astype(int) for prm in prms[1:]]
-        ns = [10 if n < 10 else n for n in ns]
-        if debug:
-            print ('ns', ns)
-        #ns [1070, 963]
-
-        # left and right coordinates for every subswath valid area
-        x1s = []
-        x2s = []
-
-        # 1st
-        xlim = prms[0].get('num_rng_bins') - ovls[0] + ns[0]
-        x1s.append(0)
-        x2s.append(xlim)
-
-        # 2nd
-        if len(prms) == 2:
-            xlim = prms[1].get('num_rng_bins') - 1
-        else:
-            # for 3 subswaths
-            xlim = prms[1].get('num_rng_bins') - ovls[1] + ns[1]
-        x1s.append(ns[0])
-        x2s.append(xlim)
-
-        # 3rd
-        if len(prms) == 3:
-            xlim = prms[2].get('num_rng_bins') - 2
-            x1s.append(ns[1])
-            x2s.append(xlim)
-
-        # check and merge SLCs
-        sumx = sum([right-left for right, left in zip(x2s, x1s)])
-        if debug:
-            print ('assert maxx == sum(...)', maxx, sumx)
-        assert maxx == sumx, 'Incorrect output grid range dimension size'
-
-        # create reference merged PRM for geocoding to extent radar coordinates calculation
-        filename = prms[0].filename[:-5] + ''.join(map(str, subswaths))
-        prm_filename = filename + '.PRM'
-        #print ('prm_filename', prm_filename)
-        prm = PRM(prms[0])
-        dt = -minh / prm.get('PRF') / 86400
-        prm = prm.set(SLC_file=None,
-                      num_lines=maxy, nrows=maxy, num_valid_az=maxy,
-                      num_rng_bins=maxx, bytes_per_line=4*maxx, good_bytes=4*maxx,
-                      SC_clock_start=prm.get('SC_clock_start') - dt,
-                      clock_start=prm.get('clock_start') - dt,
-                      SC_clock_stop=prm.get('SC_clock_start') + maxy / prm.get('PRF') / 86400,
-                      clock_stop=prm.get('clock_start') + maxy / prm.get('PRF') / 86400)\
-            .to_file(prm_filename)
-
-        return {'bottoms': bottoms, 'lefts': x1s, 'rights': x2s, 'bottom': minh, 'extent': [maxy, maxx]}
-
-    # merge_swath.c modified for SLCs
-    # use reference scene vertical subswath aligments
-    def _merge_subswaths(self, date, offsets, xlim1=None, ylim1=None, xlim2=None, ylim2=None, debug=False):
-        import xarray as xr
-        import numpy as np
-        import os
-
-        subswaths = self.get_subswaths()
-
-        #if debug:
-        #    print ('offsets', offsets)
-        # use reference scene offsets
-        maxy = offsets['extent'][0]
-        minh = offsets['bottom']
-        bottoms = offsets['bottoms']
-        maxx = offsets['extent'][1]
-        lefts = offsets['lefts']
-        rights = offsets['rights']
-
-        slcs = []
-        prms = []
-        for subswath, bottom, left, right in zip(subswaths, bottoms, lefts, rights):
-            prm = self.PRM(date, subswath=subswath)
-            # disable scaling
-            slc = prm.read_SLC_int(scale=None)
-            slc = slc.isel(x=slice(left, right)).assign_coords(y=slc.y + bottom)
-            slcs.append(slc)
-            prms.append(prm)
-
-        # check and merge SLCs, use zero fill for np.int16 datatype
-        slc = xr.concat(slcs, dim='x', fill_value=0).assign_coords(x=0.5 + np.arange(maxx))
-
-        if debug:
-            print ('assert slc.y.size == maxy', slc.y.size, maxy)
-        assert slc.y.size == maxy, 'Incorrect output grid azimuth dimension size'
-        if debug:
-            print ('assert slc.x.size == maxx', slc.x.size, maxx)
-        assert slc.x.size == maxx, 'Incorrect output grid range dimension sizes'
-        del slcs
-
-        # define merge filenames  
-        filename = prms[0].filename[:-5] + ''.join(map(str, subswaths))
-        prm_filename = filename + '.PRM'
-        netcdf_filename = os.path.basename(filename + '.grd')
-        if debug:
-            print ('prm_filename', prm_filename, 'netcdf_filename', netcdf_filename)
-
-        # merge PRM
-        prm = PRM(prms[0])
-        dt = -minh / prm.get('PRF') / 86400
-        prm = prm.set(SLC_file=None, netcdf_filename=netcdf_filename,
-                      num_lines=maxy, nrows=maxy, num_valid_az=maxy,
-                      num_rng_bins=maxx, bytes_per_line=4*maxx, good_bytes=4*maxx,
-                      SC_clock_start=prm.get('SC_clock_start') - dt,
-                      clock_start=prm.get('clock_start') - dt,
-                      SC_clock_stop=prm.get('SC_clock_start') + maxy / prm.get('PRF') / 86400,
-                      clock_stop=prm.get('clock_start') + maxy / prm.get('PRF') / 86400)\
-                 .to_file(prm_filename)
-        #.calc_dop_orb(prm.get('earth_radius'), 0, inplace=True, debug=debug)\
-
-        # add PRM to grid
-        #slc.attrs['prm'] = str(prm)
-
-        # save merged SLC
-        #prm.write_SLC_int(slc, chunksize=chunksize)
-        # encoding = {vn: self._compression(slc[vn].shape, complevel=0,
-#                     chunksize=(self.chunksize, 4*self.chunksize)) for vn in slc.data_vars}
-        encoding = {vn: self._compression(slc[vn].shape) for vn in slc.data_vars}
-        # add directory name
-        netcdf_filename = os.path.join(self.basedir, netcdf_filename)
-        # rename dimensions to prevent issue with square output
-        slc.rename({'y': 'a', 'x': 'r'})\
-            .sel(a=slice(ylim1, ylim2), r=slice(xlim1, xlim2))\
-            .to_netcdf(netcdf_filename, encoding=encoding, engine=self.netcdf_engine)
-
-        # add calculated offsets to single subswaths
-    #         for idx, prm in enumerate(prms):
-    #             prm.set(smath_maxy=maxy, swath_maxx=maxx,
-    #                     swath_bottom=bottoms[idx],
-    #                     swath_left=lefts[idx], swath_right=rights[idx]).update()
-        # for single SLCs add crop coordinates
-        #slc.attrs['bottom'] = bottom[idx]
-        #slc.attrs['left'] = left[idx]
-        #slc.attrs['right'] = rights[idx]
-
-        # cleanup
-        # [os.path.join(self.basedir, prm.get('led_file')) for prm in prms]
-        # [prm.filename for prm in prms]
-        cleanup = [os.path.join(self.basedir, prm.get('SLC_file')) for prm in prms]
-        for filename in cleanup:
-            if debug:
-                print ('DEBUG: remove', filename)
-            os.remove(filename)
-
-#     # merge_swath.c modified for SLCs
-#     # use reference scene vertical subswath aligments
-#     def aligh_subswaths(self, date, offsets, debug=False):
-#         import numpy as np
-#         import os
-# 
-#         subswaths = self.get_subswaths()
-# 
-#         #offsets = self.get_subswaths_offsets(date, offsets=offsets, debug=debug)
-#         #if debug:
-#         #    print ('offsets', offsets)
-#         # use reference scene vertical offsets
-#         maxy = offsets['extent'][0]
-#         minh = offsets['bottom']
-#         bottoms = offsets['bottoms']
-#         maxx = offsets['extent'][1]
-#         lefts = offsets['lefts']
-#         rights = offsets['rights']
-# 
-#         prms = []
-#         for subswath, bottom, left, right in zip(subswaths, bottoms, lefts, rights):
-#             prm = self.PRM(date, subswath=subswath)
-#             prms.append(prm)
-# 
-#         # define merge filenames  
-#         prm_filename = prms[0].filename[:-5] + ''.join(map(str, subswaths)) + '.PRM'
-#         if debug:
-#             print ('prm_filename', prm_filename)
-# 
-#         # merge PRM
-#         prm = PRM(prms[0])
-#         dt = -minh / prm.get('PRF') / 86400
-#         prm = prm.set(SLC_file=None,
-#                       swath_bottom=';'.join(map(str, bottoms)),
-#                       swath_left=';'.join(map(str, lefts)),
-#                       swath_right=';'.join(map(str, rights)),
-#                       num_lines=maxy, nrows=maxy, num_valid_az=maxy,
-#                       num_rng_bins=maxx, bytes_per_line=4*maxx, good_bytes=4*maxx,
-#                       SC_clock_start=prm.get('SC_clock_start') - dt,
-#                       clock_start=prm.get('clock_start') - dt,
-#                       SC_clock_stop=prm.get('SC_clock_start') + maxy / prm.get('PRF') / 86400,
-#                       clock_stop=prm.get('clock_start') + maxy / prm.get('PRF') / 86400)\
-#                  .to_file(prm_filename)
-#         #.calc_dop_orb(prm.get('earth_radius'), 0, inplace=True, debug=debug)\
-
-    def _convert_subswath(self, subswath, date, xlim1=None, ylim1=None, xlim2=None, ylim2=None, debug=False):
-        import xarray as xr
-        #import numpy as np
-        import os
-
-        prm = self.PRM(date, subswath=subswath)
-        # crop the full grid if needed
-        slc = prm.read_SLC_int(scale=None).sel(y=slice(ylim1, ylim2), x=slice(xlim1, xlim2))
-        # add PRM to grid
-        #slc.attrs['prm'] = str(prm)
-
-        netcdf_filename = prm.filename[:-4] + '.grd'
-        #print ('netcdf_filename', netcdf_filename)
-        # cleanup before saving
-        if os.path.exists(netcdf_filename):
-            os.remove(netcdf_filename)
-        # complevel=0 means disabled compression and the fastest saving
-        encoding = {vn: self._compression(slc[vn].shape) for vn in slc.data_vars}
-        # rename dimensions to prevent issue with square output
-        slc.rename({'y': 'a', 'x': 'r'})\
-            .to_netcdf(netcdf_filename, encoding=encoding, engine=self.netcdf_engine)
-
-        # cleanup
-        slc_filename = os.path.join(self.basedir, prm.get('SLC_file'))
-        if debug:
-            print ('DEBUG: remove', slc_filename)
-        # file always exists, no check required
-        os.remove(slc_filename)
-
     def baseline_table(self, n_jobs=-1, debug=False):
         """
         Generates a baseline table for Sentinel-1 data, containing dates and baseline components.
@@ -667,26 +385,6 @@ def compute_align(self, geometry='auto', dates=None, n_jobs=-1, degrees=12.0/360
             joblib.Parallel(n_jobs=n_jobs, backend=joblib_backend)(joblib.delayed(self._align_rep_subswath)(subswath, date, degrees=degrees, debug=debug) \
                                            for date in dates_rep for subswath in subswaths)
 
-        if len(subswaths) > 1:
-            # calculate the offsets and also create merged reference PRM
-            offsets = self.get_subswaths_offsets(self.reference, debug=debug)
-            # DEM extent in radar coordinates, merged reference PRM required
-            extent_ra = self.get_extent_ra()
-            minx, miny, maxx, maxy = np.round(extent_ra.bounds).astype(int)
-            #print ('minx, miny, maxx, maxy', minx, miny, maxx, maxy)
-            with self.tqdm_joblib(tqdm(desc=f'Merging Subswaths', total=len(dates))) as progress_bar:
-                joblib.Parallel(n_jobs=n_jobs, backend=joblib_backend)(joblib.delayed(self._merge_subswaths)(date, offsets, minx, miny, maxx, maxy, debug=debug) \
-                                               for date in dates)
-        else:
-            # DEM extent in radar coordinates, merged reference PRM required
-            extent_ra = self.get_extent_ra()
-            minx, miny, maxx, maxy = np.round(extent_ra.bounds).astype(int)
-            #print ('minx, miny, maxx, maxy', minx, miny, maxx, maxy)
-            # in case of a single subswath only convert SLC to NetCDF grid
-            with self.tqdm_joblib(tqdm(desc='Convert Subswath', total=len(dates))) as progress_bar:
-                joblib.Parallel(n_jobs=n_jobs, backend=joblib_backend)(joblib.delayed(self._convert_subswath)(subswaths[0], date, minx, miny, maxx, maxy, debug=debug) \
-                                               for date in dates)
-
         # merge subswaths, datapath and metapath converted to lists even for a single subswath, geometry merges bursts
         df = self.df.groupby(self.df.index).agg({'datetime': 'min', 'orbit': 'min', 'mission': 'min', 'polarization': 'min',
                                             'subswath': lambda s: int(''.join(map(str,list(s)))),
diff --git a/pygmtsar/pygmtsar/Stack_trans.py b/pygmtsar/pygmtsar/Stack_trans.py
index 4321e85..e751815 100644
--- a/pygmtsar/pygmtsar/Stack_trans.py
+++ b/pygmtsar/pygmtsar/Stack_trans.py
@@ -11,12 +11,15 @@
 from .tqdm_dask import tqdm_dask
 
 class Stack_trans(Stack_align):
-    
+
     def define_trans_grid(self, coarsen):
         import numpy as np
         # select radar coordinates extent
-        rng_max, yvalid, num_patch = self.PRM().get('num_rng_bins', 'num_valid_az', 'num_patches')
-        azi_max = yvalid * num_patch
+        #rng_max, yvalid, num_patch = self.PRM().get('num_rng_bins', 'num_valid_az', 'num_patches')
+        #azi_max = yvalid * num_patch
+        data = self.open_data(dates=[self.reference])
+        rng_max = float(data.x[-1])
+        azi_max = float(data.y[-1])
         #print ('azi_max', azi_max, 'rng_max', rng_max)
         # this grid covers the full interferogram area
         # common single pixel resolution