Improve memory efficiency in the steps of calculating ramp, deramp, residual_RMS

mintpy/ifgram_inversion.py

@@ -1002,8 +1002,7 @@ def ifgram_inversion(inps=None):
     meta['UNIT'] = 'm'
     meta['REF_DATE'] = date_list[0]
 
-    ts_obj = timeseries(inps.tsFile)
-    ts_obj.layout_hdf5(dsNameDict, meta)
+    writefile.layout_hdf5(inps.tsFile, dsNameDict, meta)
 
     # write date time-series
     date_list_utf8 = [dt.encode('utf-8') for dt in date_list]

mintpy/objects/ramp.py

@@ -75,6 +75,7 @@ def deramp(data, mask_in, ramp_type='linear', metadata=None):
     # estimate ramp
     X = np.dot(np.linalg.pinv(G[mask, :], rcond=1e-15), data[mask, :])
     ramp = np.dot(G, X)
+    ramp = np.array(ramp, dtype=data.dtype)
     del X
 
     # reference in space if metadata
@@ -85,7 +86,6 @@ def deramp(data, mask_in, ramp_type='linear', metadata=None):
 
     # do not change pixel with original zero value
     ramp[data == 0] = 0
-    ramp = np.array(ramp, dtype=data.dtype)
 
     data_out = data - ramp
     if len(dshape) == 3:

mintpy/objects/stack.py

@@ -398,13 +398,26 @@ def timeseries_rms(self, maskFile=None, outFile=None):
         """Calculate the Root Mean Square for each acquisition of time-series
             and output result to a text file.
         """
-        # Calculate RMS
-        data = self.read()
+        # Get date list
+        date_list = self.get_date_list()
+        num_date  = len(date_list)
+
+        # Get mask
         if maskFile and os.path.isfile(maskFile):
             print('read mask from file: '+maskFile)
             mask = singleDataset(maskFile).read()
-            data[:, mask == 0] = np.nan
-        self.rms = np.sqrt(np.nanmean(np.square(data), axis=(1, 2)))
+
+        # Calculate RMS one date at a time
+        self.rms = np.zeros(num_date) * np.nan
+        print('reading {} data from file: {} ...'.format(self.name, self.file))
+        prog_bar = ptime.progressBar(maxValue=num_date)
+        for i in range(num_date):
+            data = self.read(datasetName='{}'.format(date_list[i]), print_msg=False)
+            if maskFile and os.path.isfile(maskFile):
+                data[mask == 0] = np.nan
+            self.rms[i] = np.sqrt(np.nanmean(np.square(data), axis=(0, 1)))
+            prog_bar.update(i+1, suffix='{}/{}'.format(i+1, num_date))
+        prog_bar.close()
 
         # Write text file
         header = 'Root Mean Square in space for each acquisition of time-series\n'

mintpy/utils/readfile.py

@@ -796,7 +796,7 @@ def get_hdf5_dataset(name, obj):
             'cfloat': 'complex64',
         }
         data_type = atr.get('DATA_TYPE', 'none').lower()
-        if data_type is not 'none' and data_type in dataTypeDict.keys():
+        if data_type != 'none' and data_type in dataTypeDict.keys():
             atr['DATA_TYPE'] = dataTypeDict[data_type]
 
     # UNIT

mintpy/utils/utils1.py

@@ -716,11 +716,45 @@ def run_deramp(fname, ramp_type, mask_file=None, out_file=None, datasetName=None
 
     # deramping
     if k == 'timeseries':
-        print('reading data ...')
-        data = readfile.read(fname)[0]
-        print('estimating phase ramp ...')
-        data = deramp(data, mask, ramp_type=ramp_type, metadata=atr)[0]
-        writefile.write(data, out_file, ref_file=fname)
+        # initialize dataset structure
+        ts_obj = timeseries(fname)
+        ts_obj.open(print_msg=False)
+        date_list  = ts_obj.dateList
+        num_date   = len(date_list)
+        date_dtype = np.dtype('S{}'.format(len(date_list[0])))
+        dsNameDict = {
+            "date"       : (date_dtype, (num_date,)),
+            "bperp"      : (np.float32, (num_date,)),
+            "timeseries" : (np.float32, (num_date, ts_obj.length, ts_obj.width)),
+        }
+
+        # write HDF5 file with defined metadata and (empty) dataset structure
+        writefile.layout_hdf5(out_file, dsNameDict, atr, print_msg=False)
+
+        # write date time-series
+        date_list_utf8 = [dt.encode('utf-8') for dt in date_list]
+        writefile.write_hdf5_block(out_file, date_list_utf8, datasetName='date', print_msg=False)
+
+        # write bperp time-series
+        bperp = ts_obj.pbase
+        writefile.write_hdf5_block(out_file, bperp, datasetName='bperp', print_msg=False)
+
+        print('estimating phase ramp one date at a time ...')
+        prog_bar = ptime.progressBar(maxValue=num_date)
+        for i in range(num_date):
+            # read
+            dsName = 'timeseries-{}'.format(date_list[i])
+            data = readfile.read(fname, datasetName=dsName)[0]
+            # deramp
+            data = deramp(data, mask, ramp_type=ramp_type, metadata=atr)[0]
+            # write
+            writefile.write_hdf5_block(out_file, data,
+                                       datasetName='timeseries',
+                                       block=[i, i+1, 0, ts_obj.length, 0, ts_obj.width],
+                                       print_msg=False)
+            prog_bar.update(i+1, suffix='{}/{}'.format(i+1, num_date))
+        prog_bar.close()
+        print('finished writing to file: {}'.format(out_file))
 
     elif k == 'ifgramStack':
         obj = ifgramStack(fname)
@@ -734,8 +768,11 @@ def run_deramp(fname, ramp_type, mask_file=None, out_file=None, datasetName=None
                 dsOut = f[dsNameOut]
                 print('access HDF5 dataset /{}'.format(dsNameOut))
             else:
-                dsOut = f.create_dataset(dsNameOut, shape=(obj.numIfgram, obj.length, obj.width),
-                                         dtype=np.float32, chunks=True, compression=None)
+                dsOut = f.create_dataset(dsNameOut,
+                                         shape=(obj.numIfgram, obj.length, obj.width),
+                                         dtype=np.float32,
+                                         chunks=True,
+                                         compression=None)
                 print('create HDF5 dataset /{}'.format(dsNameOut))
 
             prog_bar = ptime.progressBar(maxValue=obj.numIfgram)

mintpy/utils/writefile.py

@@ -191,7 +191,7 @@ def write(datasetDict, out_file, metadata=None, ref_file=None, compression=None)
 
 #########################################################################
 
-def layout_hdf5(fname, dsNameDict, metadata):
+def layout_hdf5(fname, dsNameDict, metadata, print_msg=True):
     """Create HDF5 file with defined metadata and (empty) dataset structure
 
     Parameters: fname      - str, HDF5 file path
@@ -226,9 +226,10 @@ def layout_hdf5(fname, dsNameDict, metadata):
     }
     """
 
-    print('-'*50)
-    print('create HDF5 file {} with w mode'.format(fname))
     h5 = h5py.File(fname, "w")
+    if print_msg:
+        print('-'*50)
+        print('create HDF5 file {} with w mode'.format(fname))
 
     # initiate dataset
     for key in dsNameDict.keys():
@@ -247,9 +248,10 @@ def layout_hdf5(fname, dsNameDict, metadata):
         else:
             maxShape = data_shape
 
-        print("create dataset: {d:<25} of {t:<25} in size of {s}".format(d=key,
-                                                                         t=str(data_type),
-                                                                         s=data_shape))
+        if print_msg:
+            print("create dataset: {d:<25} of {t:<25} in size of {s}".format(d=key,
+                                                                             t=str(data_type),
+                                                                             s=data_shape))
         h5.create_dataset(key,
                           shape=data_shape,
                           maxshape=maxShape,
@@ -262,11 +264,12 @@ def layout_hdf5(fname, dsNameDict, metadata):
         h5.attrs[key] = metadata[key]
 
     h5.close()
-    print('close  HDF5 file {}'.format(fname))
+    if print_msg:
+        print('close  HDF5 file {}'.format(fname))
     return fname
 
 
-def write_hdf5_block(fname, data, datasetName, block=None, mode='a'):
+def write_hdf5_block(fname, data, datasetName, block=None, mode='a', print_msg=True):
     """Write data to existing HDF5 dataset in disk block by block.
     Parameters: data        - np.ndarray 1/2/3D matrix
                 datasetName - str, dataset name
@@ -299,11 +302,11 @@ def write_hdf5_block(fname, data, datasetName, block=None, mode='a'):
                      0, shape[2]]
 
     # write
-    print('-'*50)
-    print('open  HDF5 file {} in {} mode'.format(fname, mode))
-    with h5py.File(fname, mode) as f:
+    if print_msg:
+        print('-'*50)
+        print('open  HDF5 file {} in {} mode'.format(fname, mode))
         print("writing dataset /{:<25} block: {}".format(datasetName, block))
-
+    with h5py.File(fname, mode) as f:
         if len(block) == 6:
             f[datasetName][block[0]:block[1],
                            block[2]:block[3],
@@ -316,7 +319,8 @@ def write_hdf5_block(fname, data, datasetName, block=None, mode='a'):
         elif len(block) == 2:
             f[datasetName][block[0]:block[1]] = data
 
-    print('close HDF5 file {}.'.format(fname))
+    if print_msg:
+        print('close HDF5 file {}.'.format(fname))
     return fname