From a62275732301e34880613832e1d9068854fbfc6e Mon Sep 17 00:00:00 2001 From: Yuankailiu Date: Tue, 16 May 2023 09:52:29 -0700 Subject: [PATCH 1/4] [dev.loadion] load_data: load iono timesereis from isce2 topsStack 1. new template keywords for the ion timeseries: + mintpy.load.ionFile: ../ion_dates/*.ion #[path pattern of ionosphere timeseries files] + mintpy.load.ionBurstRampFile = ../ion_burst_ramp_merged_dates/*.float #[path pattern of ionosphere burst ramp timeseries files] 2. load_data: more flexible loading datasets + allow to specify which dataset to load with -l option. choice from {ifg/offset, geom, ion} 3. stackDict: for ion timeseries files from topsStack + class timeseriesAcqDict(): read(): data_unit & phase2range --- src/mintpy/cli/load_data.py | 7 +- src/mintpy/defaults/auto_path.py | 5 +- src/mintpy/defaults/smallbaselineApp.cfg | 5 +- src/mintpy/load_data.py | 189 ++++++++++- src/mintpy/objects/stack.py | 2 + src/mintpy/objects/stackDict.py | 392 +++++++++++++++++++++++ src/mintpy/prep_isce.py | 27 +- src/mintpy/utils/readfile.py | 4 +- 8 files changed, 596 insertions(+), 35 deletions(-) diff --git a/src/mintpy/cli/load_data.py b/src/mintpy/cli/load_data.py index e649914f6..7faf24474 100755 --- a/src/mintpy/cli/load_data.py +++ b/src/mintpy/cli/load_data.py @@ -39,10 +39,11 @@ # load & write the following HDF5 files: # ./inputs/ifgramStack.h5 for interferogram stack - # ./inputs/ionStack.h5 for ionosphere stack # ./inputs/offsetStack.h5 for range/azimuth offset stack # ./inputs/geometryRadar.h5 for geometry in radar coordinates # ./inputs/geometryGeo.h5 for geometry in geo coordinates + # ./inputs/ion.h5 for smooth ionosphere time-series (from ISCE topStack) + # ./inputs/ionBurstRamp.h5 for ionosphere burst ramp time-series (from ISCE topStack) load_data.py -t smallbaselineApp.cfg load_data.py -t smallbaselineApp.cfg GalapagosSenDT128.txt --project GalapagosSenDT128 @@ -67,8 +68,8 @@ def create_parser(subparsers=None): # input files parser.add_argument('-t', '--template', dest='template_file', type=str, nargs='+', help='template file(s) with path info.') - parser.add_argument('--geom','--geometry', dest='only_load_geometry', action='store_true', - help='Load the geometry file(s) ONLY.') + parser.add_argument('-l','--listDset', nargs='+', help='a list of datasets to be loadded (default: %(default)s)', + default=['ifg','geom','ion'], choices=['ifg','geom','ion']) # options from template file name & content parser.add_argument('--project', type=str, dest='PROJECT_NAME', diff --git a/src/mintpy/defaults/auto_path.py b/src/mintpy/defaults/auto_path.py index 14e0bbdbd..8ca8ec28d 100644 --- a/src/mintpy/defaults/auto_path.py +++ b/src/mintpy/defaults/auto_path.py @@ -23,9 +23,8 @@ mintpy.load.connCompFile = ../merged/interferograms/*/filt*.unw.conncomp mintpy.load.intFile = None -mintpy.load.ionUnwFile = ../ion/*/ion_cal/filt.ion -mintpy.load.ionCorFile = ../ion/*/ion_cal/raw_no_projection.cor -mintpy.load.ionConnCompFile = None +mintpy.load.ionFile = ../ion_dates/*.ion +mintpy.load.ionBurstRampFile = ../ion_burst_ramp_merged_dates/*.float mintpy.load.demFile = ../merged/geom_reference/hgt.rdr mintpy.load.lookupYFile = ../merged/geom_reference/lat.rdr diff --git a/src/mintpy/defaults/smallbaselineApp.cfg b/src/mintpy/defaults/smallbaselineApp.cfg index 82f56464a..f812b4e40 100644 --- a/src/mintpy/defaults/smallbaselineApp.cfg +++ b/src/mintpy/defaults/smallbaselineApp.cfg @@ -40,9 +40,8 @@ mintpy.load.connCompFile = auto #[path pattern of connected components fi mintpy.load.intFile = auto #[path pattern of wrapped interferogram files], optional mintpy.load.magFile = auto #[path pattern of interferogram magnitude files], optional ##---------ionosphere stack (optional): -mintpy.load.ionUnwFile = auto #[path pattern of unwrapped interferogram files] -mintpy.load.ionCorFile = auto #[path pattern of spatial coherence files] -mintpy.load.ionConnCompFile = auto #[path pattern of connected components files], optional but recommended +mintpy.load.ionFile = ../ion_dates/*.ion #[path pattern of ionosphere timeseries files] +mintpy.load.ionBurstRampFile = ../ion_burst_ramp_merged_dates/*.float #[path pattern of ionosphere burst ramp timeseries files] ##---------offset stack (optional): mintpy.load.azOffFile = auto #[path pattern of azimuth offset file] mintpy.load.rgOffFile = auto #[path pattern of range offset file] diff --git a/src/mintpy/load_data.py b/src/mintpy/load_data.py index 0452f1637..2d48724ef 100644 --- a/src/mintpy/load_data.py +++ b/src/mintpy/load_data.py @@ -16,11 +16,18 @@ from mintpy.objects import ( GEOMETRY_DSET_NAMES, IFGRAM_DSET_NAMES, + TIMESERIES_DSET_NAMES, geometry, ifgramStack, sensor, ) -from mintpy.objects.stackDict import geometryDict, ifgramDict, ifgramStackDict +from mintpy.objects.stackDict import ( + geometryDict, + ifgramDict, + ifgramStackDict, + timeseriesAcqDict, + timeseriesDict, +) from mintpy.utils import ptime, readfile, utils as ut ################################################################# @@ -38,9 +45,8 @@ } ION_DSET_NAME2TEMPLATE_KEY = { - 'unwrapPhase' : 'mintpy.load.ionUnwFile', - 'coherence' : 'mintpy.load.ionCorFile', - 'connectComponent': 'mintpy.load.ionConnCompFile', + 'ionosphericDelay' : 'mintpy.load.ionFile', + 'ionosphericBurstRamp' : 'mintpy.load.ionBurstRampFile', } OFF_DSET_NAME2TEMPLATE_KEY = { @@ -314,13 +320,11 @@ def read_inps_dict2ifgram_stack_dict_object(iDict, ds_name2template_key): ds_name2template_key - dict, to relate the HDF5 dataset name to the template key Returns: stackObj - ifgramStackDict object or None """ - if iDict['only_load_geometry']: + if not iDict['load_ifg']: return None if 'mintpy.load.unwFile' in ds_name2template_key.values(): obs_type = 'interferogram' - elif 'mintpy.load.ionUnwFile' in ds_name2template_key.values(): - obs_type = 'ionosphere' elif 'mintpy.load.azOffFile' in ds_name2template_key.values(): obs_type = 'offset' @@ -435,6 +439,121 @@ def read_inps_dict2ifgram_stack_dict_object(iDict, ds_name2template_key): return stackObj +def read_inps_dict2timeseries_dict_object(iDict, ds_name2tmpl_opt): + """Read input arguments into timeseriesDict object(s). + + Parameters: iDict - dict, input arguments from command line & template file + ds_name2template_key - dict, to relate the HDF5 dataset name to the template key + Returns: timeseriesObj - timeseriesDict object or None + """ + if not iDict['load_ion']: + return None + + if 'mintpy.load.unwFile' in ds_name2tmpl_opt.values(): + obs_type = 'interferogram' + elif 'mintpy.load.ionFile' in ds_name2tmpl_opt.values(): + obs_type = 'ionosphere' + elif 'mintpy.load.ionBurstRampFile' in ds_name2tmpl_opt.values(): + obs_type = 'ionosphere burst ramp' + elif 'mintpy.load.azOffFile' in ds_name2tmpl_opt.values(): + obs_type = 'offset' + + # iDict --> dsPathDict + print('-'*50) + print(f'searching {obs_type} pairs info') + print('input data files:') + max_digit = max(len(i) for i in list(ds_name2tmpl_opt.keys())) + dsPathDict = {} + for dsName in [i for i in TIMESERIES_DSET_NAMES if i in ds_name2tmpl_opt.keys()]: + key = ds_name2tmpl_opt[dsName] + if key in iDict.keys(): + files = sorted(glob.glob(str(iDict[key]))) + if len(files) > 0: + dsPathDict[dsName] = files + print(f'{dsName:<{max_digit}}: {iDict[key]}') + + # Check 1: required dataset + dsName0s = [x for x in TIMESERIES_DSET_NAMES if x in ds_name2tmpl_opt.keys()] + dsName0 = [i for i in dsName0s if i in dsPathDict.keys()] + if len(dsName0) == 0: + print(f'WARNING: No data files found for the required dataset: {dsName0s}! Skip loading for {obs_type} stack.') + return None + else: + dsName0 = dsName0[0] + + # Check 2: data dimension for unwrapPhase files + dsPathDict = skip_files_with_inconsistent_size( + dsPathDict=dsPathDict, + pix_box=iDict['box'], + dsName=dsName0) + + + + # extra metadata from observations + # e.g. EARTH_RADIUS, HEIGHT, etc. + obsMetaGeo = None + obsMetaRadar = None + + for obsName in TIMESERIES_DSET_NAMES: + if obsName in ds_name2tmpl_opt.keys(): + print(obsName) + obsFiles = sorted(glob.glob(iDict[ds_name2tmpl_opt[obsName]])) + if len(obsFiles) > 0: + atr = readfile.read_attribute(obsFiles[0]) + if 'Y_FIRST' in atr.keys(): + obsMetaGeo = atr.copy() + else: + obsMetaRadar = atr.copy() + break + + + # Check 3: number of files for all dataset types + # dsPathDict --> dsNumDict + dsNumDict = {} + for key in dsPathDict.keys(): + num_file = len(dsPathDict[key]) + dsNumDict[key] = num_file + print(f'number of {key:<{max_digit}}: {num_file}') + + dsNumList = list(dsNumDict.values()) + if any(i != dsNumList[0] for i in dsNumList): + msg = 'WARNING: NOT all types of dataset have the same number of files.' + msg += ' -> skip interferograms with missing files and continue.' + print(msg) + raise Exception(msg) + + # dsPathDict --> pairsDict --> tsObj + dsNameList = list(dsPathDict.keys()) + + + + ##################################### + + datesDict = {} + for i, dsPath0 in enumerate(dsPathDict[dsName0]): + date8s = os.path.basename(dsPath0).split('.')[0] + + tsPathDict = {} + for dsName in dsNameList: + # search the matching data file for the given date12 + dsPath1 = dsPathDict[dsName][i] + if date8s in dsPath1: + tsPathDict[dsName] = dsPath1 + + + # initiate timeseriesAcqDict object + acqObj = timeseriesAcqDict(datasetDict=tsPathDict, metadata={'data_unit': 'radian'}) + + # update datesDict object + datesDict[date8s] = acqObj + + if len(datesDict) > 0: + tsObj = timeseriesDict(datesDict=datesDict, dsName0=dsName0) + else: + tsObj = None + return tsObj + + def read_inps_dict2geometry_dict_object(iDict, dset_name2template_key): """Read input arguments into geometryDict object(s). @@ -442,6 +561,8 @@ def read_inps_dict2geometry_dict_object(iDict, dset_name2template_key): Returns: geomGeoObj - geometryDict object in geo coordinates or None geomRadarObj - geometryDict object in radar coordinates or None """ + if not iDict['load_geom']: + return None # eliminate lookup table dsName for input files in radar-coordinates if iDict['processor'] in ['isce', 'doris']: @@ -681,7 +802,8 @@ def prepare_metadata(iDict): # --dset-dir / --file-pattern obs_keys = [ 'mintpy.load.unwFile', - 'mintpy.load.ionUnwFile', + 'mintpy.load.ionFile', + 'mintpy.load.ionBurstRampFile', 'mintpy.load.rgOffFile', 'mintpy.load.azOffFile', ] @@ -827,7 +949,12 @@ def load_data(inps): # read subset info [need the metadata from above] iDict = read_subset_box(iDict) - # geometry in geo / radar coordinates + # read specific datasets + iDict['load_ifg']=True if 'ifg' in iDict['listDset'] else False + iDict['load_geom']=True if 'geom' in iDict['listDset'] else False + iDict['load_ion']=True if 'ion' in iDict['listDset'] else False + + # 3. geometry in geo / radar coordinates geom_dset_name2template_key = { **GEOM_DSET_NAME2TEMPLATE_KEY, **IFG_DSET_NAME2TEMPLATE_KEY, @@ -856,24 +983,20 @@ def load_data(inps): compression='lzf', extra_metadata=extraDict) - # observations: ifgram, ion or offset + # 4. observations: ifgram or offset # loop over obs stacks stack_ds_name2tmpl_key_list = [ IFG_DSET_NAME2TEMPLATE_KEY, - ION_DSET_NAME2TEMPLATE_KEY, OFF_DSET_NAME2TEMPLATE_KEY, ] - stack_files = ['ifgramStack.h5', 'ionStack.h5', 'offsetStack.h5'] + stack_files = ['ifgramStack.h5', 'offsetStack.h5'] stack_files = [os.path.abspath(os.path.join('./inputs', x)) for x in stack_files] for ds_name2tmpl_opt, stack_file in zip(stack_ds_name2tmpl_key_list, stack_files): # initiate dict objects stack_obj = read_inps_dict2ifgram_stack_dict_object(iDict, ds_name2tmpl_opt) - # use geom_obj as size reference while loading ionosphere geom_obj = None - if os.path.basename(stack_file).startswith('ion'): - geom_obj = geom_geo_obj if iDict['geocoded'] else geom_radar_obj # write dict objects to HDF5 files if run_or_skip(stack_file, stack_obj, iDict['box'], geom_obj=geom_obj, **kwargs) == 'run': @@ -888,6 +1011,42 @@ def load_data(inps): extra_metadata=extraDict, geom_obj=geom_obj) + + # 5. observations: ionosphericDelay, ionosphericBurstRamp + # loop over obs of timeseries + stack_ds_name2tmpl_key_list = [ + {k: v} for k, v in ION_DSET_NAME2TEMPLATE_KEY.items() + ] + ts_files = ['ion.h5', 'ionBurstRamp.h5'] + ts_files = [os.path.abspath(os.path.join('./inputs', x)) for x in ts_files] + for (ds_name2tmpl_opt, ts_file) in zip(stack_ds_name2tmpl_key_list, ts_files): + # rename the object key to timeseries + for k in ION_DSET_NAME2TEMPLATE_KEY.keys(): + if k in ds_name2tmpl_opt: + ds_name2tmpl_opt[TIMESERIES_DSET_NAMES[0]] = ds_name2tmpl_opt[k] + del ds_name2tmpl_opt[k] + + # initiate dict objects (a new func getting timeseries_obj) + ts_obj = read_inps_dict2timeseries_dict_object(iDict, ds_name2tmpl_opt) + + # use geom_obj as size reference while loading ionosphere + geom_obj = None + if os.path.basename(ts_file).startswith('ion'): + geom_obj = geom_geo_obj if iDict['geocoded'] else geom_radar_obj + + # write dict objects to HDF5 files + if run_or_skip(ts_file, ts_obj, iDict['box'], geom_obj=geom_obj, **kwargs) == 'run': + ts_obj.write2hdf5( + outputFile=ts_file, + access_mode='w', + box=iDict['box'], + xstep=iDict['xstep'], + ystep=iDict['ystep'], + mli_method=iDict['method'], + compression=iDict['compression'], + extra_metadata=extraDict, + geom_obj=geom_obj) + # used time m, s = divmod(time.time()-start_time, 60) print(f'time used: {m:02.0f} mins {s:02.1f} secs.\n') diff --git a/src/mintpy/objects/stack.py b/src/mintpy/objects/stack.py index 56bcd012f..ac6919c40 100644 --- a/src/mintpy/objects/stack.py +++ b/src/mintpy/objects/stack.py @@ -39,6 +39,8 @@ 'raw', 'troposphericDelay', 'topographicResidual', + 'ionosphericDelay', + 'ionosphericBurstRamp', 'ramp', 'displacement', ] diff --git a/src/mintpy/objects/stackDict.py b/src/mintpy/objects/stackDict.py index 54c340f1d..8b9ecf0b2 100644 --- a/src/mintpy/objects/stackDict.py +++ b/src/mintpy/objects/stackDict.py @@ -25,6 +25,7 @@ DATA_TYPE_DICT, GEOMETRY_DSET_NAMES, IFGRAM_DSET_NAMES, + TIMESERIES_DSET_NAMES, ) from mintpy.utils import attribute as attr, ptime, readfile, utils0 as ut @@ -411,6 +412,397 @@ def get_metadata(self, family=IFGRAM_DSET_NAMES[0]): return self.metadata +######################################################################################## +class timeseriesDict: + """ + timeseriesStack object for a set of InSAR acquisitions from the same platform and track. + + Example: + from mintpy.objects.insarobj import timeseriesDict + datesDict = {('20160524'):acqObj1, + ('20160524'):acqObj2, + ('20160524'):acqObj3, + ('20160530'):acqObj4, + ... + } + tsObj = timeseriesDict(datesDict=datesDict) + tsObj.write2hdf5(outputFile='timeseries.h5', box=(200,500,300,600)) + """ + + def __init__(self, name='timeseries', datesDict=None, dsName0=TIMESERIES_DSET_NAMES[0]): + self.name = name + self.datesDict = datesDict + self.dsName0 = dsName0 #reference dataset name, unwrapPhase OR azimuthOffset OR rangeOffset + + def get_size(self, box=None, xstep=1, ystep=1, geom_obj=None): + """Get size in 3D""" + num_date = len(self.datesDict) + acqObj = [v for v in self.datesDict.values()][0] + length, width = acqObj.get_size(family=self.dsName0) + + # use the reference geometry obj size + # for low-reso ionosphere from isce2/topsStack + if geom_obj: + length, width = geom_obj.get_size() + + # update due to subset + if box: + length, width = box[3] - box[1], box[2] - box[0] + + # update due to multilook + length = length // ystep + width = width // xstep + + return num_date, length, width + + def get_date_list(self): + self.dateList = list(self.datesDict.keys()) + return self.dateList + + def get_dataset_list(self): + acqObj = [x for x in self.datesDict.values()][0] + dsetList = [x for x in acqObj.datasetDict.keys()] + return dsetList + + def get_metadata(self): + acqObj = [v for v in self.datesDict.values()][0] + self.metadata = acqObj.get_metadata(family=self.dsName0) + return self.metadata + + def get_dataset_data_type(self, dsName): + acqObj = [v for v in self.datesDict.values()][0] + dsFile = acqObj.datasetDict[dsName] + metadata = readfile.read_attribute(dsFile) + dataType = DATA_TYPE_DICT[metadata.get('DATA_TYPE', 'float32').lower()] + return dataType + + def write2hdf5(self, outputFile='timeseries.h5', access_mode='w', box=None, xstep=1, ystep=1, mli_method='nearest', + compression=None, extra_metadata=None, geom_obj=None): + """Save/write an timeseriesDict object into an HDF5 file with the structure defined in: + + https://mintpy.readthedocs.io/en/latest/api/data_structure/#ifgramstack (Kai need to update the doc?) + + Parameters: outputFile - str, Name of the HDF5 file for the InSAR stack + access_mode - str, access mode of output File, e.g. w, r+ + box - tuple, subset range in (x0, y0, x1, y1) + x/ystep - int, multilook number in x/y direction + mli_method - str, multilook method, nearest, mean or median + compression - str, HDF5 dataset compression method, None, lzf or gzip + extra_metadata - dict, extra metadata to be added into output file + geom_obj - geometryDict object, size reference to determine the resizing operation. + Returns: outputFile - str, Name of the HDF5 file for the InSAR stack + """ + print('-'*50) + + # output directory + output_dir = os.path.dirname(outputFile) + if not os.path.isdir(output_dir): + os.makedirs(output_dir) + print(f'create directory: {output_dir}') + + self.dates = sorted(self.get_date_list()) # Kai + self.dsNames = list(self.datesDict[self.dates[0]].datasetDict.keys()) # Kai + self.dsNames = [i for i in TIMESERIES_DSET_NAMES if i in self.dsNames] + maxDigit = max(len(i) for i in self.dsNames) + num_date, length, width = self.get_size( + box=box, + xstep=xstep, + ystep=ystep) + + # check if resize is needed for a lower resolution stack, e.g. ionosphere from isce2/topsStack + resize2shape = None + if geom_obj and os.path.basename(outputFile).startswith('ion'): + # compare the original data size between ionosphere and geometry w/o subset/multilook + ion_size = self.get_size()[1:] + geom_size = geom_obj.get_size() + if ion_size != geom_size: + msg = 'lower resolution ionosphere file detected' + msg += f' --> resize from {ion_size} to {geom_size} via skimage.transform.resize ...' + print(msg) + + # matrix shape for the original geometry size w/o subset/multilook + resize2shape = geom_size + # data size of the output HDF5 file w/ resize/subset/multilook + length, width = self.get_size( + box=box, + xstep=xstep, + ystep=ystep, + geom_obj=geom_obj)[1:] + + # write HDF5 file + with h5py.File(outputFile, access_mode) as f: + print(f'create HDF5 file {outputFile} with {access_mode} mode') + + ############################### + # 3D datasets containing unwrapPhase, magnitude, coherence, connectComponent, wrapPhase, etc. + for dsName in self.dsNames: + dsShape = (num_date, length, width) + dsDataType = np.float32 + dsCompression = compression + if dsName in ['connectComponent']: + dsDataType = np.int16 + dsCompression = 'lzf' + mli_method = 'nearest' + + print(('create dataset /{d:<{w}} of {t:<25} in size of {s}' + ' with compression = {c}').format(d=dsName, + w=maxDigit, + t=str(dsDataType), + s=dsShape, + c=dsCompression)) + ds = f.create_dataset(dsName, + shape=dsShape, + maxshape=(None, dsShape[1], dsShape[2]), + dtype=dsDataType, + chunks=True, + compression=dsCompression) + + # set no-data value - printout msg + if dsName.endswith('OffsetVar'): + print(f'set no-data value for {dsName} from 99 to NaN.') + dsFile = self.datesDict[self.dates[0]].datasetDict[dsName] + if dsFile.endswith('cov.bip'): + print('convert variance to standard deviation.') + + # msg + if xstep * ystep > 1: + print(f'apply {xstep} x {ystep} multilooking/downsampling via {mli_method} ...') + + prog_bar = ptime.progressBar(maxValue=num_date) + for i, date in enumerate(self.dates): + prog_bar.update(i+1, suffix=f'{date[0]}') + + # read and/or resize + acqObj = self.datesDict[date] + data = acqObj.read(dsName, + box=box, + xstep=xstep, + ystep=ystep, + mli_method=mli_method, + resize2shape=resize2shape)[0] + + # write + ds[i, :, :] = data + + ds.attrs['MODIFICATION_TIME'] = str(time.time()) + prog_bar.close() + + ############################### + # 2D dataset containing reference and secondary dates of all dates + dsName = 'date' + dsDataType = np.string_ + dsShape = (num_date, 2) + print('create dataset /{d:<{w}} of {t:<25} in size of {s}'.format(d=dsName, + w=maxDigit, + t=str(dsDataType), + s=dsShape)) + data = np.array(self.dates, dtype=dsDataType) + f.create_dataset(dsName, data=data) + + ############################### + # 1D dataset containing perpendicular baseline of all dates + if False: # Mute this since we don't have it when reading timeseries dataset + dsName = 'bperp' + dsDataType = np.float32 + dsShape = (num_date,) + print('create dataset /{d:<{w}} of {t:<25} in size of {s}'.format(d=dsName, + w=maxDigit, + t=str(dsDataType), + s=dsShape)) + # get bperp + data = np.zeros(num_date, dtype=dsDataType) + for i in range(num_date): + acqObj = self.datesDict[self.dates[i]] + data[i] = acqObj.get_perp_baseline(family=self.dsName0) + # write + f.create_dataset(dsName, data=data) + + ############################### + # 1D dataset containing bool value of dropping the interferograms or not + dsName = 'dropIfgram' # Kai: need to delete this? + dsDataType = np.bool_ + dsShape = (num_date,) + print('create dataset /{d:<{w}} of {t:<25} in size of {s}'.format(d=dsName, + w=maxDigit, + t=str(dsDataType), + s=dsShape)) + data = np.ones(dsShape, dtype=dsDataType) + f.create_dataset(dsName, data=data) + + ############################### + # Attributes + # read metadata from original data file w/o resize/subset/multilook + meta = self.get_metadata() + if extra_metadata: + meta.update(extra_metadata) + print(f'add extra metadata: {extra_metadata}') + + # update metadata due to resize + # for low resolution ionosphere from isce2/topsStack + if resize2shape: + print('update metadata due to resize') + meta = attr.update_attribute4resize(meta, resize2shape) + + # update metadata due to subset + if box: + print('update metadata due to subset') + meta = attr.update_attribute4subset(meta, box) + + # update metadata due to multilook + if xstep * ystep > 1: + print('update metadata due to multilook') + meta = attr.update_attribute4multilook(meta, ystep, xstep) + + # write metadata to HDF5 file at the root level + meta['FILE_TYPE'] = self.name + for key, value in meta.items(): + f.attrs[key] = value + + print(f'Finished writing to {outputFile}') + return outputFile + + +######################################################################################## +class timeseriesAcqDict: + """ + Timeseries object for timeseries, date, bperp, ... from the same platform and track. + + Example: + from mintpy.utils import readfile + from mintpy.utils.insarobj import timeseriesAcqDict + datasetDict = {'timeseries' :'$PROJECT_DIR/ion/*.rdr', + 'date' :'$PROJECT_DIR/', + 'bperp' :bperpDict + ... + } + bperpDict = {'20160406':'$PROJECT_DIR/merged/baselines/20160406/bperp', + '20160418':'$PROJECT_DIR/merged/baselines/20160418/bperp', + ... + } + metadata = readfile.read_attribute('$PROJECT_DIR/merged/interferograms/20160629_20160723/filt_fine.unw') + acqObj = timeseriesAcqDict(processor='isce', datasetDict=datasetDict, extraMetadata=metadata) + """ + + def __init__(self, name='timeseries', datasetDict={}, metadata=None): + self.name = name + self.datasetDict = datasetDict + + self.platform = None + self.track = None + self.processor = None + # platform, track and processor can get values from metadata if they exist + if metadata is not None: + for key, value in metadata.items(): + setattr(self, key, value) + + def read(self, family, box=None, xstep=1, ystep=1, mli_method='nearest', resize2shape=None): + """Read data for the given dataset name. + + Parameters: self - ifgramDict object + family - str, dataset name + box - tuple of 4 int, in (x0, y0, x1, y1) with respect to the full resolution + x/ystep - int, number of pixels to skip, with respect to the full resolution + mli_method - str, interpolation method, nearest, mean, median + resize2shape - tuple of 2 int, resize the native matrix to the given shape + Set to None for not resizing + Returns: data - 2D np.ndarray + meta - dict, metadata + """ + self.file = self.datasetDict[family] + dirname = os.path.dirname(self.file).split('/')[-1] + box2read = None if resize2shape else box + + # 1. read input file + data, meta = readfile.read(self.file, + datasetName=family, + box=box2read, + xstep=1, + ystep=1) + + # 2. resize + if resize2shape: + # link: https://scikit-image.org/docs/dev/api/skimage.transform.html#skimage.transform.resize + data = resize(data, + output_shape=resize2shape, + order=1, + mode='constant', + anti_aliasing=True, + preserve_range=True) + + # 3. subset by box + if box: + data = data[box[1]:box[3], + box[0]:box[2]] + + # 4. multilook + if xstep * ystep > 1: + if mli_method == 'nearest': + # multilook - nearest resampling + # output data size + xsize = int(data.shape[1] / xstep) + ysize = int(data.shape[0] / ystep) + # sampling + data = data[int(ystep/2)::ystep, + int(xstep/2)::xstep] + data = data[:ysize, :xsize] + + else: + # multilook - mean or median resampling + data = multilook_data(data, + lks_y=ystep, + lks_x=xstep, + method=mli_method) + + # 5. check the input unit + if self.data_unit: + if self.data_unit == 'm': data *= 1e0 + elif self.data_unit == 'cm': data *= 1e-2 + elif self.data_unit == 'mm': data *= 1e-3 + elif self.data_unit == 'radian': + if dirname in ['ion_dates', 'ion_burst_ramp_merged_dates']: + phase2range = 1 * float(meta['WAVELENGTH']) / (4.*np.pi) + else: + phase2range = -1 * float(meta['WAVELENGTH']) / (4.*np.pi) + data *= phase2range + + return data, meta + + def get_size(self, family=TIMESERIES_DSET_NAMES[0]): + self.file = self.datasetDict[family] + metadata = readfile.read_attribute(self.file) + length = int(metadata['LENGTH']) + width = int(metadata['WIDTH']) + return length, width + + def get_perp_baseline(self, family=TIMESERIES_DSET_NAMES[0]): + self.file = self.datasetDict[family] + metadata = readfile.read_attribute(self.file) + self.bperp_top = float(metadata['P_BASELINE_TOP_HDR']) + self.bperp_bottom = float(metadata['P_BASELINE_BOTTOM_HDR']) + self.bperp = (self.bperp_top + self.bperp_bottom) / 2.0 + return self.bperp + + def get_metadata(self, family=TIMESERIES_DSET_NAMES[0]): + self.file = self.datasetDict[family] + self.metadata = readfile.read_attribute(self.file) + self.length = int(self.metadata['LENGTH']) + self.width = int(self.metadata['WIDTH']) + self.metadata['UNIT'] = 'm' + + if self.track: + self.metadata['TRACK'] = self.track + + if self.platform: + self.metadata['PLATFORM'] = self.platform + + return self.metadata + + + def get_dataset_list(self): + self.datasetList = list(self.datasetDict.keys()) + return self.datasetList + + ######################################################################################## class geometryDict: """ diff --git a/src/mintpy/prep_isce.py b/src/mintpy/prep_isce.py index 3159a48a3..424323708 100644 --- a/src/mintpy/prep_isce.py +++ b/src/mintpy/prep_isce.py @@ -153,19 +153,28 @@ def prepare_stack(obs_file, metadata=dict(), baseline_dict=dict(), update_mode=T print_msg = True if num_file > 5 else False # do not print progress bar for <=5 files prog_bar = ptime.progressBar(maxValue=num_file, print_msg=print_msg) for i, isce_file in enumerate(isce_files): - # get date1/2 - date12 = ptime.get_date12_from_path(isce_file) - dates = ptime.yyyymmdd(date12.replace('-','_').split('_')) - prog_bar.update(i+1, suffix=f'{dates[0]}_{dates[1]} {i+1}/{num_file}') + if obs_file.endswith('ion_dates/*.ion') or obs_file.endswith('ion_burst_ramp_merged_dates/*.float'): + # get date for timeseries + date = ptime.get_date_str_format(isce_file) + prog_bar.update(i+1, suffix=f'{date} {i+1}/{num_file}') + ts_meta = {**meta} + ts_meta.update(readfile.read_attribute(isce_file, metafile_ext='.xml')) + out_meta = ts_meta + else: + # get date1/2 + date12 = ptime.get_date12_from_path(isce_file) + dates = ptime.yyyymmdd(date12.replace('-','_').split('_')) + prog_bar.update(i+1, suffix=f'{dates[0]}_{dates[1]} {i+1}/{num_file}') - # merge metadata from: data.rsc, *.unw.xml and DATE12/P_BASELINE_TOP/BOTTOM_HDR - ifg_meta = {**meta} - ifg_meta.update(readfile.read_attribute(isce_file, metafile_ext='.xml')) - ifg_meta = add_ifgram_metadata(ifg_meta, dates, baseline_dict) + # merge metadata from: data.rsc, *.unw.xml and DATE12/P_BASELINE_TOP/BOTTOM_HDR + ifg_meta = {**meta} + ifg_meta.update(readfile.read_attribute(isce_file, metafile_ext='.xml')) + ifg_meta = add_ifgram_metadata(ifg_meta, dates, baseline_dict) + out_meta = ifg_meta # write .rsc file rsc_file = isce_file+'.rsc' - writefile.write_roipac_rsc(ifg_meta, rsc_file, + writefile.write_roipac_rsc(out_meta, rsc_file, update_mode=update_mode, print_msg=False) diff --git a/src/mintpy/utils/readfile.py b/src/mintpy/utils/readfile.py index 24c1c029c..5f467cfee 100644 --- a/src/mintpy/utils/readfile.py +++ b/src/mintpy/utils/readfile.py @@ -585,7 +585,7 @@ def read_binary_file(fname, datasetName=None, box=None, xstep=1, ystep=1): if datasetName: if datasetName.startswith(('mag', 'amp')): cpx_band = 'magnitude' - elif datasetName in ['phase', 'angle']: + elif any(x in datasetName.lower() for x in ['phase', 'angle']): cpx_band = 'phase' elif datasetName.lower() == 'real': cpx_band = 'real' @@ -771,7 +771,7 @@ def get_slice_list(fname, no_complex=False): with h5py.File(fname, 'r') as f: d1_list = [i for i in f.keys() if isinstance(f[i], h5py.Dataset)] - if ftype == 'timeseries' and ftype in d1_list: + if ftype == 'timeseries' and any(x.startswith(ftype) for x in d1_list): obj = timeseries(fname) obj.open(print_msg=False) slice_list = obj.sliceList From cbda18c5cffc499203f7bd82377263299cbb4521 Mon Sep 17 00:00:00 2001 From: Yuankailiu Date: Thu, 8 Jun 2023 17:20:38 -0700 Subject: [PATCH 2/4] [dev.loadion] cli/load_data examples for -l option --- src/mintpy/cli/load_data.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/mintpy/cli/load_data.py b/src/mintpy/cli/load_data.py index 7faf24474..654be7929 100755 --- a/src/mintpy/cli/load_data.py +++ b/src/mintpy/cli/load_data.py @@ -49,7 +49,11 @@ # load geometry ONLY smallbaselineApp.py SaltonSeaSenDT173.txt -g - load_data.py -t smallbaselineApp.cfg --geom + load_data.py -t smallbaselineApp.cfg -l geom + + # load ionoshpere time series ONLY + smallbaselineApp.py SaltonSeaSenDT173.txt -g + load_data.py -t smallbaselineApp.cfg -l ion """ From 7510d076bc9bdf19517d4eea0cb2252cdf5fe3e6 Mon Sep 17 00:00:00 2001 From: Yuankailiu Date: Wed, 20 Sep 2023 16:52:08 -0700 Subject: [PATCH 3/4] [dev.loadion] rebase on main branch, fix Codacy issues --- src/mintpy/load_data.py | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/src/mintpy/load_data.py b/src/mintpy/load_data.py index 2d48724ef..2e100940a 100644 --- a/src/mintpy/load_data.py +++ b/src/mintpy/load_data.py @@ -487,26 +487,6 @@ def read_inps_dict2timeseries_dict_object(iDict, ds_name2tmpl_opt): pix_box=iDict['box'], dsName=dsName0) - - - # extra metadata from observations - # e.g. EARTH_RADIUS, HEIGHT, etc. - obsMetaGeo = None - obsMetaRadar = None - - for obsName in TIMESERIES_DSET_NAMES: - if obsName in ds_name2tmpl_opt.keys(): - print(obsName) - obsFiles = sorted(glob.glob(iDict[ds_name2tmpl_opt[obsName]])) - if len(obsFiles) > 0: - atr = readfile.read_attribute(obsFiles[0]) - if 'Y_FIRST' in atr.keys(): - obsMetaGeo = atr.copy() - else: - obsMetaRadar = atr.copy() - break - - # Check 3: number of files for all dataset types # dsPathDict --> dsNumDict dsNumDict = {} From dd6dbe99cb43981b679be1990dc9171edb75399b Mon Sep 17 00:00:00 2001 From: Yuankailiu Date: Tue, 26 Sep 2023 10:05:19 -0700 Subject: [PATCH 4/4] [dev.loadion] load_data: always read geom file --- src/mintpy/load_data.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mintpy/load_data.py b/src/mintpy/load_data.py index 2e100940a..2bb752cc3 100644 --- a/src/mintpy/load_data.py +++ b/src/mintpy/load_data.py @@ -930,8 +930,8 @@ def load_data(inps): iDict = read_subset_box(iDict) # read specific datasets + iDict['load_geom']=True iDict['load_ifg']=True if 'ifg' in iDict['listDset'] else False - iDict['load_geom']=True if 'geom' in iDict['listDset'] else False iDict['load_ion']=True if 'ion' in iDict['listDset'] else False # 3. geometry in geo / radar coordinates