diff --git a/mintpy/ifgram_inversion.py b/mintpy/ifgram_inversion.py index db7d912a6..f0c64a001 100755 --- a/mintpy/ifgram_inversion.py +++ b/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] diff --git a/mintpy/objects/ramp.py b/mintpy/objects/ramp.py index dce861421..fc91f02f9 100644 --- a/mintpy/objects/ramp.py +++ b/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: diff --git a/mintpy/objects/stack.py b/mintpy/objects/stack.py index 3f7a92ab8..535ef30e7 100644 --- a/mintpy/objects/stack.py +++ b/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' diff --git a/mintpy/utils/readfile.py b/mintpy/utils/readfile.py index 677b8403a..963dfd099 100644 --- a/mintpy/utils/readfile.py +++ b/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 diff --git a/mintpy/utils/utils1.py b/mintpy/utils/utils1.py index 477a81b75..862c3cd72 100644 --- a/mintpy/utils/utils1.py +++ b/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) diff --git a/mintpy/utils/writefile.py b/mintpy/utils/writefile.py index 45e7337d9..1b36132fa 100644 --- a/mintpy/utils/writefile.py +++ b/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