From 39621983f0a30724e618980fc186df73616b858f Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Fri, 17 Feb 2023 13:37:37 -0800 Subject: [PATCH 01/13] added loading of correction layers --- src/mintpy/cli/prep_aria.py | 17 ++- src/mintpy/prep_aria.py | 241 ++++++++++++++++++++++++++++++++++++ 2 files changed, 257 insertions(+), 1 deletion(-) diff --git a/src/mintpy/cli/prep_aria.py b/src/mintpy/cli/prep_aria.py index 720cea0e1..367c723db 100755 --- a/src/mintpy/cli/prep_aria.py +++ b/src/mintpy/cli/prep_aria.py @@ -49,7 +49,8 @@ prep_aria.py -t SanFranSenDT42.txt prep_aria.py -s ../stack/ -d ../DEM/SRTM_3arcsec.dem -i '../incidenceAngle/*.vrt' prep_aria.py -s ../stack/ -d ../DEM/SRTM_3arcsec.dem -i '../incidenceAngle/*.vrt' -a '../azimuthAngle/*.vrt' -w ../mask/watermask.msk - + prep_aria.py -s ../stack/ -d ../DEM/SRTM_3arcsec.dem -i '../incidenceAngle/*.vrt' -tropo '../stack/troposphereTotal/GMAO_stack.vrt' -ion '../stack/ionStack.vrt' + # download / extract / prepare inteferograms stack from ARIA using ARIA-tools: # reference: https://github.com/aria-tools/ARIA-tools ariaDownload.py -b '37.25 38.1 -122.6 -121.75' --track 42 @@ -106,6 +107,20 @@ def create_parser(subparsers=None): help='Name of the azimuth angle file.') geom.add_argument('-w','--water-mask', dest='waterMaskFile', type=str, help='Name of the water mask file') + + # correction layers: troposphereTotal, ionosphere, solidEarthTides + corr = parser.add_argument_group('corrections') + corr.add_argument('-ct', '--tropo', dest='tropoFile', type=str, + help='Name of the Troposhere Delay stack file', default=None) + corr.add_argument('-ci', '--ion', dest='ionoFile', type=str, + help='Name of the Ionosphere Delay stack file', default=None) + corr.add_argument('-cs', '--set', dest='setFile', type=str, + help='Name of the Solid Earth Tides stack file', default=None) + corr.add_argument('--cluster', dest='cluster', type=str, choices={'local', 'pbs', None}, + help='Parallelize inversion of correction layers w Dask', default=None) + corr.add_argument('-n', '--num-workers', dest='num_workers', type=str, + help='Dask number of workers', default='2') + return parser diff --git a/src/mintpy/prep_aria.py b/src/mintpy/prep_aria.py index 8246976e7..359c396f6 100644 --- a/src/mintpy/prep_aria.py +++ b/src/mintpy/prep_aria.py @@ -424,6 +424,200 @@ def write_ifgram_stack(outfile, unwStack, cohStack, connCompStack, ampStack=None dsAmp = None return outfile +# OPTIONAL - ARIA corrections troposphereTotal, ionosphere, solidearthtides +def write_correction(outfile, corrStack, box=None, + xstep=1, ystep=1, mli_method='nearest'): + """Write corrections to object ifgramStack HDF5 file from stack VRT files + Correction layers are in a form of differential observations for each interferometric pair + + ARIA correction layers: + troposhereTotal : models GMAO, HRRR, HRES, ERA5 '/science/grids/corrections/external/troposphere/' + ionosphere '/science/grids/corrections/derived/ionosphere/ionosphere' + solidEarthTides '/science/grids/corrections/derived/solidearthtides/' + """ + + print('-'*50) + max_digit = len(os.path.basename(str(corrStack))) + if corrStack is not None: + print('open {f:<{w}} with gdal ...'.format(f=os.path.basename(corrStack), w=max_digit)) + + dsCor = gdal.Open(corrStack, gdal.GA_ReadOnly) + # get the layer name (for tropo this will get the model name) + layer = dsCor.GetRasterBand(1).GetMetadataDomainList()[0] + + # extract NoDataValue (from the last */date2_date1.vrt file for example) + ds = gdal.Open(dsCor.GetFileList()[-1], gdal.GA_ReadOnly) + noDataValue = ds.GetRasterBand(1).GetNoDataValue() + print(f'grab NoDataValue for {layer}: {noDataValue:<5} and convert to 0.') + + # sort the order of correction layer pairs based on date1_date2 with date1 < date2 + nPairs = dsCor.RasterCount + d12BandDict = {} + for ii in range(nPairs): + bnd = dsCor.GetRasterBand(ii+1) + d12 = bnd.GetMetadata(layer)["Dates"] + d12 = sorted(d12.split("_")) + d12 = f'{d12[0]}_{d12[1]}' + d12BandDict[d12] = ii+1 + d12List = sorted(d12BandDict.keys()) + print(f'number of {layer} pairs: {len(d12List)}') + + # box to gdal arguments + # link: https://gdal.org/python/osgeo.gdal.Band-class.html#ReadAsArray + if box is not None: + kwargs = dict( + xoff=box[0], + yoff=box[1], + win_xsize=box[2]-box[0], + win_ysize=box[3]-box[1]) + else: + kwargs = dict() + + if xstep * ystep > 1: + msg = f'apply {xstep} x {ystep} multilooking/downsampling via {mli_method} to {layer}' + print(msg) + print(f'writing data to HDF5 file {outfile} with a mode ...') + with h5py.File(outfile, "a") as f: + + prog_bar = ptime.progressBar(maxValue=nPairs) + for ii in range(nPairs): + d12 = d12List[ii] + bndIdx = d12BandDict[d12] + prog_bar.update(ii+1, suffix=f'{d12} {ii+1}/{nPairs}') + + f["date"][ii,0] = d12.split("_")[0].encode("utf-8") + f["date"][ii,1] = d12.split("_")[1].encode("utf-8") + f["dropIfgram"][ii] = True + + bnd = dsCor.GetRasterBand(bndIdx) + data = bnd.ReadAsArray(**kwargs) + data = multilook_data(data, ystep, xstep, method=mli_method) + data[data == noDataValue] = 0 #assign pixel with no-data to 0 + data[np.isnan(data)] = 0 #assign nan pixel to 0 + data *= -1.0 #date2_date1 -> date1_date2 + f["unwrapPhase"][ii,:,:] = data + # NOTE: ifgramStack accepts only 'unwrapPhase', 'rangeOffset', 'azimuthOffset' + # if changed to different layer name, changes need be be done in + # stack.py line: 782 + # reference_point.py line: 45 + # ifgram_inversion.py line 651 + # add option to recognize the layer + + bperp = float(bnd.GetMetadata(str(layer))["perpendicularBaseline"]) + bperp *= -1.0 #date2_date1 -> date1_date2 + f["bperp"][ii] = bperp + + prog_bar.close() + + # add MODIFICATION_TIME metadata to each 3D dataset + for dsName in ['unwrapPhase']: + f[dsName].attrs['MODIFICATION_TIME'] = str(time.time()) + + print(f'finished writing to HD5 file: {outfile}') + ds = None + dsCor = None + + return outfile + +def invert_diff_corrections(input_filename, output_filename, dataset, + cluster = None, num_workers = '4', waterMask=None, + maskDataset = None, maskThreshold = 0.4): + ''' + Invert differential ARIA correction layers to get correction for SAR acquistion dates + - the inversion gives corrections relative to the REF_DATE (typically the first acquisition) + + NOTE: 1.the inversion network needs to be connected, otherwise inversion will give wrong estimates + for the isolated clusters + 2. each SAR date needs to have min two datasets with that date to have min_degree of freedom >=2 + otherwise inversion will give wrong estimates for those dates + ''' + + # create inps dummy + class dummy(): + pass + + # PREPARE INPUT OBJECT + inps = dummy() + # input + inps.ifgramStackFile = input_filename + inps.obsDatasetName = 'unwrapPhase' #only this available for use at the moment + inps.skip_ref = False + + # solver + inps.minNormVelocity = False + inps.minRedundancy = 1.0 + # Note: minRedun set to 2.0 gives wrong estim, and dask has some errors + inps.weightFunc = 'no' + inps.calcCov = False + + # mask + inps.waterMaskFile = waterMask + inps.maskDataset = maskDataset + inps.maskThreshold = maskThreshold + + # cluster - Expose / leave None for now + inps.cluster = cluster + inps.maxMemory = 2 + inps.numWorker = num_workers + inps.config = 'local' #not sure what to put here + + # outputs + #Avoid setting 'no' for invQualityFile + #self.invQualityFile = 'modelTempCoh.h5' # not needed, maybe leave it to avoid dask outout issues + inps.invQualityFile = 'modelTempCoh.h' + inps.numInvFile = 'numInvModel.h5' + inps.tsFile = output_filename + + ### INVERSION + from mintpy.utils import readfile + from mintpy.cli import reference_point + from mintpy.ifgram_inversion import run_ifgram_inversion + + # remove if exists + try: + os.remove(inps.tsFile) + print('Delete existing file') + except FileNotFoundError: + print("File is not present in the system.") + + # 1. get reference pixel + + reference_point.main([inps.ifgramStackFile]) + + # 2. invert differential model obs + run_ifgram_inversion(inps) + + # 3. compensate for range2phase conversion as not needed for models + if dataset.startswith(('tropo', 'set')): + print('Return back units to original') + data, metadata = readfile.read(inps.tsFile, datasetName='timeseries') + phase2range = -1 * float(metadata['WAVELENGTH']) / (4.*np.pi) + + #Replace values + with h5py.File(inps.tsFile, 'r+') as f: + f['timeseries'][:] = data / phase2range + + # 4. clean up uncessary files + os.remove(inps.invQualityFile) + os.remove(inps.numInvFile) + + +def get_correction_layer(correction_filename): + ds = gdal.Open(correction_filename, gdal.GA_ReadOnly) + # get the layer name (for tropo this will get the model name) + layer_name = ds.GetRasterBand(1).GetMetadataDomainList()[0] + + # Get type of correction + if layer_name in ['GMAO', 'HRES', 'HRRR', "ERA5"]: + layer_type = 'tropo' + else: + # ionosphere, solid earth tides + layer_type = layer_name + + #close + ds = None + + return layer_name, layer_type #################################################################################### def load_aria(inps): @@ -520,6 +714,53 @@ def load_aria(inps): ystep=inps.ystep, ) + ########## output file 3 - correction layers + correction_layers = [inps.tropoFile, inps.ionoFile, inps.setFile] + + # define correction dataset structure for ifgramStack + # is it already resampled?? + ds_name_dict = { + "date" : (np.dtype('S8'), (num_pair, 2)), + "dropIfgram" : (np.bool_, (num_pair,)), + "bperp" : (np.float32, (num_pair,)), + "unwrapPhase" : (np.float32, (num_pair, length, width)), + } + meta['FILE_TYPE'] = 'ifgramStack' + + # Loop thourgh defined correction layers + for layer in correction_layers: + if layer: + layer_name, layer_type = get_correction_layer(layer) + writefile.layout_hdf5( + f'./inputs/d{layer_name}.h5', + ds_name_dict, + metadata=meta, + compression=inps.compression, + ) + + # write data to disk + write_correction( + f'./inputs/d{layer_name}.h5', + corrStack=layer, + box=box, + xstep=inps.xstep, + ystep=inps.ystep, + ) + + + # invert layer to get correction + # for SAR acquistion dates + invert_diff_corrections(f'./inputs/d{layer_name}.h5', + f'./inputs/{layer_name}.h5', + layer_type, + cluster=inps.cluster, + num_workers=inps.num_workers, + waterMask = None, + maskDataset = None, + maskThreshold = 0.4) + + + print('-'*50) # used time From 76880063563ec42008e3747e5b5fbe88ec7cfb27 Mon Sep 17 00:00:00 2001 From: Simran S Sangha Date: Wed, 15 Mar 2023 19:33:32 -0700 Subject: [PATCH 02/13] Add support for updated epoch stacks --- src/mintpy/prep_aria.py | 130 ++++++++++++++++++++++++++++------------ 1 file changed, 92 insertions(+), 38 deletions(-) diff --git a/src/mintpy/prep_aria.py b/src/mintpy/prep_aria.py index 359c396f6..685802164 100644 --- a/src/mintpy/prep_aria.py +++ b/src/mintpy/prep_aria.py @@ -5,6 +5,7 @@ ############################################################ +import datetime as dt import os import time @@ -715,50 +716,103 @@ def load_aria(inps): ) ########## output file 3 - correction layers - correction_layers = [inps.tropoFile, inps.ionoFile, inps.setFile] - # define correction dataset structure for ifgramStack - # is it already resampled?? - ds_name_dict = { - "date" : (np.dtype('S8'), (num_pair, 2)), - "dropIfgram" : (np.bool_, (num_pair,)), - "bperp" : (np.float32, (num_pair,)), - "unwrapPhase" : (np.float32, (num_pair, length, width)), - } - meta['FILE_TYPE'] = 'ifgramStack' + # Invert Iono stack and write out cube + if inps.ionoFile: + # define correction dataset structure for ifgramStack + ds_name_dict = { + 'date' : (np.dtype('S8'), (num_pair, 2)), + 'dropIfgram' : (np.bool_, (num_pair,)), + 'bperp' : (np.float32, (num_pair,)), + 'unwrapPhase' : (np.float32, (num_pair, length, width)), + } + meta['FILE_TYPE'] = 'ifgramStack' - # Loop thourgh defined correction layers + layer_name, layer_type = get_correction_layer(inps.ionoFile) + writefile.layout_hdf5( + f'./inputs/d{layer_name}.h5', + ds_name_dict, + metadata=meta, + compression=inps.compression, + ) + + # write data to disk + write_correction( + f'./inputs/d{layer_name}.h5', + corrStack=inps.ionoFile, + box=box, + xstep=inps.xstep, + ystep=inps.ystep, + ) + + + # invert layer to get correction + # for SAR acquistion dates + invert_diff_corrections(f'./inputs/d{layer_name}.h5', + f'./inputs/{layer_name}.h5', + layer_type, + cluster=inps.cluster, + num_workers=inps.num_workers, + waterMask = None, + maskDataset = None, + maskThreshold = 0.4) + + # Loop through other correction layers also provided as epochs + correction_layers = [inps.tropoFile, inps.setFile] for layer in correction_layers: if layer: + # get name and type layer_name, layer_type = get_correction_layer(layer) - writefile.layout_hdf5( - f'./inputs/d{layer_name}.h5', - ds_name_dict, - metadata=meta, - compression=inps.compression, - ) - - # write data to disk - write_correction( - f'./inputs/d{layer_name}.h5', - corrStack=layer, - box=box, - xstep=inps.xstep, - ystep=inps.ystep, - ) - - - # invert layer to get correction - # for SAR acquistion dates - invert_diff_corrections(f'./inputs/d{layer_name}.h5', - f'./inputs/{layer_name}.h5', - layer_type, - cluster=inps.cluster, - num_workers=inps.num_workers, - waterMask = None, - maskDataset = None, - maskThreshold = 0.4) + # open raster + ds = gdal.Open(layer, gdal.GA_ReadOnly) + + # get times and arrays + date_list = [] + dt_objs = [] + num_dates = ds.RasterCount + ts_arr = np.zeros((num_dates, length, width), dtype=np.float32) + for ii in range(num_dates): + # read band + bnd = ds.GetRasterBand(ii+1) + bnd_arr = bnd.ReadAsArray() + bnd_name = bnd.GetMetadataDomainList()[0] + d12 = bnd.GetMetadata(bnd_name)['Dates'] + utc = bnd.GetMetadata(bnd_name)['UTCTime (HH:MM:SS.ss)'] + # set time + utc = time.strptime(utc, "%H:%M:%S.%f") + center_line_utc = utc.tm_hour*3600.0 + \ + utc.tm_min*60.0 + utc.tm_sec + utc_sec = dt.timedelta(seconds=float(center_line_utc)) + dt_obj = dt.datetime.strptime(d12, '%Y%m%d') + utc_sec + date_list.append(d12) + dt_objs.append(dt_obj) + # set array + phase2range = float(meta['WAVELENGTH']) / (4.*np.pi) + ts_arr[ii,:,:] = bnd_arr * phase2range + ds = None + + # adjust attribute + meta['FILE_TYPE'] = 'timeseries' + meta['DATA_TYPE'] = 'float32' + meta['UNIT'] = 'm' + for key in ['REF_Y', 'REF_X', 'REF_DATE']: + if key in meta.keys(): + meta.pop(key) + + # define correction dataset structure for timeseries + ds_name_dict = { + 'timeseries' : ts_arr, + 'date' : np.array(date_list, dtype=np.string_), + 'sensingMid' : np.array( + [i.strftime('%Y%m%dT%H%M%S') for i in dt_objs], + dtype=np.string_), + } + + # write + writefile.write(ds_name_dict, + out_file=f'./inputs/{layer_name}.h5', + metadata=meta) print('-'*50) From 338ec1d404b45ab20acf5747e7496f82dd43342f Mon Sep 17 00:00:00 2001 From: Simran S Sangha Date: Thu, 16 Mar 2023 13:36:54 -0700 Subject: [PATCH 03/13] Mask out 0s --- src/mintpy/prep_aria.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mintpy/prep_aria.py b/src/mintpy/prep_aria.py index 685802164..53576cfc8 100644 --- a/src/mintpy/prep_aria.py +++ b/src/mintpy/prep_aria.py @@ -791,6 +791,8 @@ def load_aria(inps): phase2range = float(meta['WAVELENGTH']) / (4.*np.pi) ts_arr[ii,:,:] = bnd_arr * phase2range ds = None + # mask nodata + ts_arr[ts_arr == 0] = np.nan # adjust attribute meta['FILE_TYPE'] = 'timeseries' From 43869a092b67daa1a5378b6236ab5c31befadad5 Mon Sep 17 00:00:00 2001 From: mgovorcin Date: Wed, 22 Mar 2023 14:28:59 -0700 Subject: [PATCH 04/13] updating correction func to be consistent w others --- src/mintpy/cli/prep_aria.py | 4 +- src/mintpy/prep_aria.py | 253 +++++++++++++++++++++++------------- 2 files changed, 166 insertions(+), 91 deletions(-) diff --git a/src/mintpy/cli/prep_aria.py b/src/mintpy/cli/prep_aria.py index 367c723db..65f3b8ebe 100755 --- a/src/mintpy/cli/prep_aria.py +++ b/src/mintpy/cli/prep_aria.py @@ -49,7 +49,7 @@ prep_aria.py -t SanFranSenDT42.txt prep_aria.py -s ../stack/ -d ../DEM/SRTM_3arcsec.dem -i '../incidenceAngle/*.vrt' prep_aria.py -s ../stack/ -d ../DEM/SRTM_3arcsec.dem -i '../incidenceAngle/*.vrt' -a '../azimuthAngle/*.vrt' -w ../mask/watermask.msk - prep_aria.py -s ../stack/ -d ../DEM/SRTM_3arcsec.dem -i '../incidenceAngle/*.vrt' -tropo '../stack/troposphereTotal/GMAO_stack.vrt' -ion '../stack/ionStack.vrt' + prep_aria.py -s ../stack/ -d ../DEM/SRTM_3arcsec.dem -i '../incidenceAngle/*.vrt' --tropo '../stack/troposphereTotal/GMAO_stack.vrt' --iono '../stack/ionStack.vrt' # download / extract / prepare inteferograms stack from ARIA using ARIA-tools: # reference: https://github.com/aria-tools/ARIA-tools @@ -112,7 +112,7 @@ def create_parser(subparsers=None): corr = parser.add_argument_group('corrections') corr.add_argument('-ct', '--tropo', dest='tropoFile', type=str, help='Name of the Troposhere Delay stack file', default=None) - corr.add_argument('-ci', '--ion', dest='ionoFile', type=str, + corr.add_argument('-ci', '--iono', dest='ionoFile', type=str, help='Name of the Ionosphere Delay stack file', default=None) corr.add_argument('-cs', '--set', dest='setFile', type=str, help='Name of the Solid Earth Tides stack file', default=None) diff --git a/src/mintpy/prep_aria.py b/src/mintpy/prep_aria.py index 53576cfc8..bf0e94114 100644 --- a/src/mintpy/prep_aria.py +++ b/src/mintpy/prep_aria.py @@ -5,12 +5,12 @@ ############################################################ -import datetime as dt import os import time import h5py import numpy as np +import datetime as dt try: from osgeo import gdal @@ -297,7 +297,7 @@ def write_geometry(outfile, demFile, incAngleFile, azAngleFile=None, waterMaskFi # write f['waterMask'][:,:] = water_mask - print(f'finished writing to HD5 file: {outfile}') + print(f'finished writing to HD5 file: {outfile}\n') return outfile @@ -418,32 +418,28 @@ def write_ifgram_stack(outfile, unwStack, cohStack, connCompStack, ampStack=None for dsName in ['unwrapPhase','coherence','connectComponent']: f[dsName].attrs['MODIFICATION_TIME'] = str(time.time()) - print(f'finished writing to HD5 file: {outfile}') + print(f'finished writing to HD5 file: {outfile}\n') dsUnw = None dsCoh = None dsComp = None dsAmp = None return outfile + # OPTIONAL - ARIA corrections troposphereTotal, ionosphere, solidearthtides -def write_correction(outfile, corrStack, box=None, +def write_iono_stack(outfile, ionStack, box=None, xstep=1, ystep=1, mli_method='nearest'): - """Write corrections to object ifgramStack HDF5 file from stack VRT files - Correction layers are in a form of differential observations for each interferometric pair - - ARIA correction layers: - troposhereTotal : models GMAO, HRRR, HRES, ERA5 '/science/grids/corrections/external/troposphere/' - ionosphere '/science/grids/corrections/derived/ionosphere/ionosphere' - solidEarthTides '/science/grids/corrections/derived/solidearthtides/' + """Write ionospheric corrections to object ifgramStack HDF5 file from stack VRT files + Ionospheric layer is in a form of differential observations for each interferometric pair + ARIA_GUNW_NC_PATH : ionosphere '/science/grids/corrections/derived/ionosphere/ionosphere' """ - print('-'*50) - max_digit = len(os.path.basename(str(corrStack))) - if corrStack is not None: - print('open {f:<{w}} with gdal ...'.format(f=os.path.basename(corrStack), w=max_digit)) + max_digit = len(os.path.basename(str(ionStack))) + if ionStack is not None: + print('open {f:<{w}} with gdal ...'.format(f=os.path.basename(ionStack), w=max_digit)) - dsCor = gdal.Open(corrStack, gdal.GA_ReadOnly) - # get the layer name (for tropo this will get the model name) + dsCor = gdal.Open(ionStack, gdal.GA_ReadOnly) + # get the layer name layer = dsCor.GetRasterBand(1).GetMetadataDomainList()[0] # extract NoDataValue (from the last */date2_date1.vrt file for example) @@ -514,12 +510,109 @@ def write_correction(outfile, corrStack, box=None, for dsName in ['unwrapPhase']: f[dsName].attrs['MODIFICATION_TIME'] = str(time.time()) - print(f'finished writing to HD5 file: {outfile}') + print(f'finished writing to HD5 file: {outfile}\n') ds = None dsCor = None return outfile + +def write_model_stack(outfile, corrStack, box=None, + xstep=1, ystep=1, mli_method='nearest'): + """Write SET and TropsphericDelay corrections to HDF5 file from stack VRT files + Correction layers are stored for each SAR acquisition date + + ARIA_GUNW_NC_PATH: + troposhereTotal : models GMAO, HRRR, HRES, ERA5 '/science/grids/corrections/external/troposphere/' + solidEarthTides '/science/grids/corrections/derived/solidearthtides/' + """ + + print('-'*50) + max_digit = len(os.path.basename(str(corrStack))) + + if corrStack is not None: + print('open {f:<{w}} with gdal ...'.format(f=os.path.basename(corrStack), w=max_digit)) + + # open raster + dsCor = gdal.Open(corrStack, gdal.GA_ReadOnly) + # extract NoDataValue (from the last date.vrt file for example) + ds = gdal.Open(dsCor.GetFileList()[-1], gdal.GA_ReadOnly) + noDataValue = ds.GetRasterBand(1).GetNoDataValue() + ds = None + + # get the layer name (for tropo this will get the model name) + layer = dsCor.GetRasterBand(1).GetMetadataDomainList()[0] + + # Get the wavelength. need to convert radians to meters + wavelength = np.float64(dsCor.GetRasterBand(1).GetMetadata(layer)["Wavelength (m)"]) + phase2range = -1 * wavelength / (4.*np.pi) + + # get model dates and time + nDate = dsCor.RasterCount + dateDict = {} + sensingDict = {} + for ii in range(nDate): + bnd = dsCor.GetRasterBand(ii+1) + date = bnd.GetMetadata(layer)["Dates"] + utc = dt.datetime.strptime(date + ',' + \ + bnd.GetMetadata(layer)["UTCTime (HH:MM:SS.ss)"], + "%Y%m%d,%H:%M:%S.%f") + dateDict[date] = ii+1 + sensingDict[ii+1] = str(utc) + dateList = sorted(dateDict.keys()) + print(f'number of {layer} datasets: {len(dateList)}') + + # box to gdal arguments + # link: https://gdal.org/python/osgeo.gdal.Band-class.html#ReadAsArray + if box is not None: + kwargs = dict( + xoff=box[0], + yoff=box[1], + win_xsize=box[2]-box[0], + win_ysize=box[3]-box[1]) + else: + kwargs = dict() + + if xstep * ystep > 1: + msg = f'apply {xstep} x {ystep} multilooking/downsampling via {mli_method} to {layer}' + print(msg) + + print(f'writing data to HDF5 file {outfile} with a mode ...') + with h5py.File(outfile, "a") as f: + prog_bar = ptime.progressBar(maxValue=nDate) + for ii in range(nDate): + date = dateList[ii] + bndIdx = dateDict[date] + utc = sensingDict[bndIdx] + prog_bar.update(ii+1, suffix=f'{date} {ii+1}/{nDate}') + + f["date"][ii] = date.encode("utf-8") + f["sensingMid"][ii] = utc.encode("utf-8") + + bnd = dsCor.GetRasterBand(bndIdx) + data = bnd.ReadAsArray(**kwargs) + data = multilook_data(data, ystep, xstep, method=mli_method) + data[data == noDataValue] = 0 #assign pixel with no-data to 0 + data[np.isnan(data)] = 0 #assign nan pixel to 0 + f["timeseries"][ii,:,:] = data * phase2range + + prog_bar.close() + + # add MODIFICATION_TIME metadata to each 3D dataset + for dsName in ['timeseries']: + f[dsName].attrs['MODIFICATION_TIME'] = str(time.time()) + + print(f'finished writing to HD5 file: {outfile}\n') + dsCor = None + + return outfile + + +def get_number_of_epochs(vrtfile): + ds = gdal.Open(vrtfile, gdal.GA_ReadOnly) + + return ds.RasterCount + def invert_diff_corrections(input_filename, output_filename, dataset, cluster = None, num_workers = '4', waterMask=None, maskDataset = None, maskThreshold = 0.4): @@ -717,6 +810,7 @@ def load_aria(inps): ########## output file 3 - correction layers + # 3.1 - ionosphere # Invert Iono stack and write out cube if inps.ionoFile: # define correction dataset structure for ifgramStack @@ -729,98 +823,79 @@ def load_aria(inps): meta['FILE_TYPE'] = 'ifgramStack' layer_name, layer_type = get_correction_layer(inps.ionoFile) - writefile.layout_hdf5( - f'./inputs/d{layer_name}.h5', - ds_name_dict, - metadata=meta, - compression=inps.compression, - ) - - # write data to disk - write_correction( - f'./inputs/d{layer_name}.h5', - corrStack=inps.ionoFile, - box=box, - xstep=inps.xstep, - ystep=inps.ystep, - ) - - - # invert layer to get correction - # for SAR acquistion dates - invert_diff_corrections(f'./inputs/d{layer_name}.h5', - f'./inputs/{layer_name}.h5', - layer_type, - cluster=inps.cluster, - num_workers=inps.num_workers, - waterMask = None, - maskDataset = None, - maskThreshold = 0.4) + if run_or_skip(inps, ds_name_dict, out_file=inps.outfile[0]) == 'run': + writefile.layout_hdf5( + f'{out_dir}/d{layer_name}.h5', + ds_name_dict, + metadata=meta, + compression=inps.compression, + ) + + # write data to disk + write_iono_stack( + f'{out_dir}/d{layer_name}.h5', + ionStack=inps.ionoFile, + box=box, + xstep=inps.xstep, + ystep=inps.ystep, + ) + + # invert layer to get ionosphere phase delay + # on SAR acquistion dates + invert_diff_corrections(f'{out_dir}/d{layer_name}.h5', + f'{out_dir}/{layer_name}.h5', + layer_type, + cluster=inps.cluster, + num_workers=inps.num_workers, + waterMask = None, + maskDataset = None, + maskThreshold = 0.4) + + # 3.2 - model based corrections: SolidEarthTides and Troposphere # Loop through other correction layers also provided as epochs + correction_layers = [inps.tropoFile, inps.setFile] for layer in correction_layers: if layer: # get name and type layer_name, layer_type = get_correction_layer(layer) + num_dates = get_number_of_epochs(layer) - # open raster - ds = gdal.Open(layer, gdal.GA_ReadOnly) - - # get times and arrays - date_list = [] - dt_objs = [] - num_dates = ds.RasterCount - ts_arr = np.zeros((num_dates, length, width), dtype=np.float32) - for ii in range(num_dates): - # read band - bnd = ds.GetRasterBand(ii+1) - bnd_arr = bnd.ReadAsArray() - bnd_name = bnd.GetMetadataDomainList()[0] - d12 = bnd.GetMetadata(bnd_name)['Dates'] - utc = bnd.GetMetadata(bnd_name)['UTCTime (HH:MM:SS.ss)'] - # set time - utc = time.strptime(utc, "%H:%M:%S.%f") - center_line_utc = utc.tm_hour*3600.0 + \ - utc.tm_min*60.0 + utc.tm_sec - utc_sec = dt.timedelta(seconds=float(center_line_utc)) - dt_obj = dt.datetime.strptime(d12, '%Y%m%d') + utc_sec - date_list.append(d12) - dt_objs.append(dt_obj) - # set array - phase2range = float(meta['WAVELENGTH']) / (4.*np.pi) - ts_arr[ii,:,:] = bnd_arr * phase2range - ds = None - # mask nodata - ts_arr[ts_arr == 0] = np.nan - - # adjust attribute meta['FILE_TYPE'] = 'timeseries' - meta['DATA_TYPE'] = 'float32' meta['UNIT'] = 'm' + meta['DATA_TYPE'] = 'float32' for key in ['REF_Y', 'REF_X', 'REF_DATE']: if key in meta.keys(): meta.pop(key) # define correction dataset structure for timeseries ds_name_dict = { - 'timeseries' : ts_arr, - 'date' : np.array(date_list, dtype=np.string_), - 'sensingMid' : np.array( - [i.strftime('%Y%m%dT%H%M%S') for i in dt_objs], - dtype=np.string_), + 'date' : (np.dtype('S8'), (num_dates, )), + 'sensingMid' : (np.dtype('S15'), (num_dates, )), + 'timeseries' : (np.float32, (num_dates, length, width)), } - # write - writefile.write(ds_name_dict, - out_file=f'./inputs/{layer_name}.h5', - metadata=meta) - - + if run_or_skip(inps, ds_name_dict, out_file=inps.outfile[0]) == 'run': + writefile.layout_hdf5( + f'{out_dir}/{layer_name}.h5', + ds_name_dict, + metadata=meta, + compression=inps.compression, + ) + + # write data to disk + write_model_stack( + f'{out_dir}/{layer_name}.h5', + corrStack=layer, + box=box, + xstep=inps.xstep, + ystep=inps.ystep, + ) print('-'*50) # used time m, s = divmod(time.time() - start_time, 60) print(f'time used: {m:02.0f} mins {s:02.1f} secs.') - return + From 90bebdbed16cf7be7a798204ef7b40067a57eb4c Mon Sep 17 00:00:00 2001 From: Simran S Sangha Date: Tue, 28 Mar 2023 10:54:05 -0700 Subject: [PATCH 05/13] Adjust sign error and handle multiple tropo files --- src/mintpy/cli/prep_aria.py | 10 +++++++--- src/mintpy/prep_aria.py | 12 +++++++----- 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/src/mintpy/cli/prep_aria.py b/src/mintpy/cli/prep_aria.py index 65f3b8ebe..65ab4c249 100755 --- a/src/mintpy/cli/prep_aria.py +++ b/src/mintpy/cli/prep_aria.py @@ -161,10 +161,14 @@ def cmd_line_parse(iargs = None): # search for wildcard pattern fnames = glob.glob(iDict[key]) if iDict[key] else [] - # user the first element if more than one exist + # return the first element if more than one exist + # except for tropo, for which multiple inputs could be passed if len(fnames) > 0: - iDict[key] = fnames[0] - print('{k:<{w}} : {f}'.format(k=key, w=max_digit, f=fnames[0])) + if 'tropo' not in key: + iDict[key] = fnames[0] + else: + iDict[key] = fnames + print('{k:<{w}} : {f}'.format(k=key, w=max_digit, f=iDict[key])) elif key in required_ds_keys: # raise exception if any required DS is missing diff --git a/src/mintpy/prep_aria.py b/src/mintpy/prep_aria.py index bf0e94114..e2cd17845 100644 --- a/src/mintpy/prep_aria.py +++ b/src/mintpy/prep_aria.py @@ -545,7 +545,7 @@ def write_model_stack(outfile, corrStack, box=None, # Get the wavelength. need to convert radians to meters wavelength = np.float64(dsCor.GetRasterBand(1).GetMetadata(layer)["Wavelength (m)"]) - phase2range = -1 * wavelength / (4.*np.pi) + phase2range = wavelength / (4.*np.pi) # get model dates and time nDate = dsCor.RasterCount @@ -854,8 +854,10 @@ def load_aria(inps): # 3.2 - model based corrections: SolidEarthTides and Troposphere # Loop through other correction layers also provided as epochs - - correction_layers = [inps.tropoFile, inps.setFile] + # handle multiple tropo stacks (if specified) + if inps.tropoFile is None: + inps.tropoFile = [None] + correction_layers = inps.tropoFile + [inps.setFile] for layer in correction_layers: if layer: # get name and type @@ -878,7 +880,7 @@ def load_aria(inps): if run_or_skip(inps, ds_name_dict, out_file=inps.outfile[0]) == 'run': writefile.layout_hdf5( - f'{out_dir}/{layer_name}.h5', + f'{out_dir}/{layer_name}_ARIA.h5', ds_name_dict, metadata=meta, compression=inps.compression, @@ -886,7 +888,7 @@ def load_aria(inps): # write data to disk write_model_stack( - f'{out_dir}/{layer_name}.h5', + f'{out_dir}/{layer_name}_ARIA.h5', corrStack=layer, box=box, xstep=inps.xstep, From 92bc5a3bde7182dd909752a1c7e659bdda0e1f0d Mon Sep 17 00:00:00 2001 From: Simran S Sangha Date: Tue, 13 Aug 2024 13:13:26 -0700 Subject: [PATCH 06/13] Update cli examples in prep_aria --- src/mintpy/cli/prep_aria.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mintpy/cli/prep_aria.py b/src/mintpy/cli/prep_aria.py index 65ab4c249..c0ae14bd7 100755 --- a/src/mintpy/cli/prep_aria.py +++ b/src/mintpy/cli/prep_aria.py @@ -49,7 +49,7 @@ prep_aria.py -t SanFranSenDT42.txt prep_aria.py -s ../stack/ -d ../DEM/SRTM_3arcsec.dem -i '../incidenceAngle/*.vrt' prep_aria.py -s ../stack/ -d ../DEM/SRTM_3arcsec.dem -i '../incidenceAngle/*.vrt' -a '../azimuthAngle/*.vrt' -w ../mask/watermask.msk - prep_aria.py -s ../stack/ -d ../DEM/SRTM_3arcsec.dem -i '../incidenceAngle/*.vrt' --tropo '../stack/troposphereTotal/GMAO_stack.vrt' --iono '../stack/ionStack.vrt' + prep_aria.py -s ../stack/ -d ../DEM/SRTM_3arcsec.dem -i '../incidenceAngle/*.vrt' --set '../stack/setStack.vrt' --tropo '../stack/troposphereTotal/HRRR_stack.vrt' --iono '../stack/ionStack.vrt' # download / extract / prepare inteferograms stack from ARIA using ARIA-tools: # reference: https://github.com/aria-tools/ARIA-tools From 66f486fc14aeb6b0b680e86108930616fe4add5d Mon Sep 17 00:00:00 2001 From: rzinke Date: Sun, 27 Oct 2024 20:15:11 -0400 Subject: [PATCH 07/13] Reuse write_ifgram_stack and iono_split_spectrum to reduce maintainance burden --- src/mintpy/prep_aria.py | 343 ++++++++++------------------------------ 1 file changed, 83 insertions(+), 260 deletions(-) diff --git a/src/mintpy/prep_aria.py b/src/mintpy/prep_aria.py index b4c1260e0..d7b838786 100644 --- a/src/mintpy/prep_aria.py +++ b/src/mintpy/prep_aria.py @@ -301,52 +301,43 @@ def write_geometry(outfile, demFile, incAngleFile, azAngleFile=None, waterMaskFi return outfile -def write_ifgram_stack(outfile, unwStack, cohStack, connCompStack, ampStack=None, - box=None, xstep=1, ystep=1, mli_method='nearest'): - """Write ifgramStack HDF5 file from stack VRT files +def write_ifgram_stack(outfile, stackFiles, box=None, xstep=1, ystep=1, mli_method='nearest'): + """Write stacks to HDF5 files from stack VRT files """ print('-'*50) - stackFiles = [unwStack, cohStack, connCompStack, ampStack] - max_digit = max(len(os.path.basename(str(i))) for i in stackFiles) - for stackFile in stackFiles: - if stackFile is not None: - print('open {f:<{w}} with gdal ...'.format(f=os.path.basename(stackFile), w=max_digit)) - dsUnw = gdal.Open(unwStack, gdal.GA_ReadOnly) - dsCoh = gdal.Open(cohStack, gdal.GA_ReadOnly) - dsComp = gdal.Open(connCompStack, gdal.GA_ReadOnly) - if ampStack is not None: - dsAmp = gdal.Open(ampStack, gdal.GA_ReadOnly) - else: - dsAmp = None + # remove None entries + stackFiles = {key:val for key, val in stackFiles.items() if val is not None} - # extract NoDataValue (from the last */date2_date1.vrt file for example) - ds = gdal.Open(dsUnw.GetFileList()[-1], gdal.GA_ReadOnly) - noDataValueUnw = ds.GetRasterBand(1).GetNoDataValue() - print(f'grab NoDataValue for unwrapPhase : {noDataValueUnw:<5} and convert to 0.') + # check all files exist + for dsName, fname in stackFiles.items(): + if not os.path.exists(fname): + raise Exception("%s does not exist" % fname) - ds = gdal.Open(dsCoh.GetFileList()[-1], gdal.GA_ReadOnly) - noDataValueCoh = ds.GetRasterBand(1).GetNoDataValue() - print(f'grab NoDataValue for coherence : {noDataValueCoh:<5} and convert to 0.') + # determine field length for printing + max_digit = max(len(os.path.basename(str(i))) for i in stackFiles.values()) + for stackFile in stackFiles.values(): + if stackFile is not None: + print('open {f:<{w}} with gdal ...'.format(f=os.path.basename(stackFile), w=max_digit)) - ds = gdal.Open(dsComp.GetFileList()[-1], gdal.GA_ReadOnly) - noDataValueComp = ds.GetRasterBand(1).GetNoDataValue() - print(f'grab NoDataValue for connectComponent: {noDataValueComp:<5} and convert to 0.') - ds = None + # extract NoDataValue for each stack (from the last */date2_date1.vrt file for example) + noDataValues = {} + for dsName in stackFiles.keys(): + dsStack = gdal.Open(stackFiles[dsName], gdal.GA_ReadOnly) + ds = gdal.Open(dsStack.GetFileList()[-1], gdal.GA_ReadOnly) + noDataValues[dsName] = ds.GetRasterBand(1).GetNoDataValue() - if dsAmp is not None: - ds = gdal.Open(dsAmp.GetFileList()[-1], gdal.GA_ReadOnly) - noDataValueAmp = ds.GetRasterBand(1).GetNoDataValue() - print(f'grab NoDataValue for magnitude : {noDataValueAmp:<5} and convert to 0.') + fileName = os.path.basename(stackFiles[dsName]) + print(f'grab NoDataValue for {fileName:<{max_digit}}: {noDataValues[dsName]:<5} and convert to 0.') ds = None # sort the order of interferograms based on date1_date2 with date1 < date2 - nPairs = dsUnw.RasterCount + nPairs = dsStack.RasterCount d12BandDict = {} for ii in range(nPairs): - bnd = dsUnw.GetRasterBand(ii+1) - d12 = bnd.GetMetadata("unwrappedPhase")["Dates"] + bnd = dsStack.GetRasterBand(ii+1) + d12 = bnd.GetMetadata(bnd.GetMetadataDomainList()[0])["Dates"] d12 = sorted(d12.split("_")) d12 = f'{d12[0]}_{d12[1]}' d12BandDict[d12] = ii+1 @@ -364,14 +355,9 @@ def write_ifgram_stack(outfile, unwStack, cohStack, connCompStack, ampStack=None else: kwargs = dict() - if xstep * ystep > 1: - msg = f'apply {xstep} x {ystep} multilooking/downsampling via {mli_method} to: unwrapPhase, coherence' - msg += ', magnitude' if dsAmp is not None else '' - msg += f'\napply {xstep} x {ystep} multilooking/downsampling via nearest to: connectComponent' - print(msg) + # write to HDF5 file print(f'writing data to HDF5 file {outfile} with a mode ...') with h5py.File(outfile, "a") as f: - prog_bar = ptime.progressBar(maxValue=nPairs) for ii in range(nPairs): d12 = d12List[ii] @@ -382,141 +368,61 @@ def write_ifgram_stack(outfile, unwStack, cohStack, connCompStack, ampStack=None f["date"][ii,1] = d12.split("_")[1].encode("utf-8") f["dropIfgram"][ii] = True - bnd = dsUnw.GetRasterBand(bndIdx) - data = bnd.ReadAsArray(**kwargs) - data = multilook_data(data, ystep, xstep, method=mli_method) - data[data == noDataValueUnw] = 0 #assign pixel with no-data to 0 - data *= -1.0 #date2_date1 -> date1_date2 - f["unwrapPhase"][ii,:,:] = data - - bperp = float(bnd.GetMetadata("unwrappedPhase")["perpendicularBaseline"]) - bperp *= -1.0 #date2_date1 -> date1_date2 - f["bperp"][ii] = bperp - - bnd = dsCoh.GetRasterBand(bndIdx) - data = bnd.ReadAsArray(**kwargs) - data = multilook_data(data, ystep, xstep, method=mli_method) - data[data == noDataValueCoh] = 0 #assign pixel with no-data to 0 - f["coherence"][ii,:,:] = data - - bnd = dsComp.GetRasterBand(bndIdx) - data = bnd.ReadAsArray(**kwargs) - data = multilook_data(data, ystep, xstep, method='nearest') - data[data == noDataValueComp] = 0 #assign pixel with no-data to 0 - f["connectComponent"][ii,:,:] = data - - if dsAmp is not None: - bnd = dsAmp.GetRasterBand(bndIdx) + # loop through stacks + print(stackFiles.keys()) + for dsName in stackFiles.keys(): + dsStack = gdal.Open(stackFiles[dsName], gdal.GA_ReadOnly) + bnd = dsStack.GetRasterBand(bndIdx) data = bnd.ReadAsArray(**kwargs) - data = multilook_data(data, ystep, xstep, method=mli_method) - data[data == noDataValueAmp] = 0 #assign pixel with no-data to 0 - f["magnitude"][ii,:,:] = data + if xstep * ystep > 1: + mli_method_spec = mli_method if stackName not in \ + ['connCompStack'] else 'nearest' + print(f'apply {xstep} x {ystep} multilooking/downsampling via ' + f'{mli_method_spec} to: {stackName}') + data = multilook_data(data, ystep, xstep, method=mli_method_spec) + data[data == noDataValues[dsName]] = 0 #assign pixel with no-data to 0 - prog_bar.close() - - # add MODIFICATION_TIME metadata to each 3D dataset - for dsName in ['unwrapPhase','coherence','connectComponent']: - f[dsName].attrs['MODIFICATION_TIME'] = str(time.time()) + if dsName == 'unwrapPhase': + data *= -1 # date2_date1 -> date1_date2 + f['unwrapPhase'][ii,:,:] = data - print(f'finished writing to HD5 file: {outfile}\n') - dsUnw = None - dsCoh = None - dsComp = None - dsAmp = None - return outfile + bperp = float(bnd.GetMetadata("unwrappedPhase")["perpendicularBaseline"]) + bperp *= -1.0 # date2_date1 -> date1_date2 + f["bperp"][ii] = bperp + elif dsName == 'coherence': + f["coherence"][ii,:,:] = data -# OPTIONAL - ARIA corrections troposphereTotal, ionosphere, solidearthtides -def write_iono_stack(outfile, ionStack, box=None, - xstep=1, ystep=1, mli_method='nearest'): - """Write ionospheric corrections to object ifgramStack HDF5 file from stack VRT files - Ionospheric layer is in a form of differential observations for each interferometric pair - ARIA_GUNW_NC_PATH : ionosphere '/science/grids/corrections/derived/ionosphere/ionosphere' - """ - print('-'*50) - max_digit = len(os.path.basename(str(ionStack))) - if ionStack is not None: - print('open {f:<{w}} with gdal ...'.format(f=os.path.basename(ionStack), w=max_digit)) + elif dsName == 'connectComponent': + f["connectComponent"][ii,:,:] = data - dsCor = gdal.Open(ionStack, gdal.GA_ReadOnly) - # get the layer name - layer = dsCor.GetRasterBand(1).GetMetadataDomainList()[0] + elif dsName == 'magnitude': + f["magnitude"][ii,:,:] = data - # extract NoDataValue (from the last */date2_date1.vrt file for example) - ds = gdal.Open(dsCor.GetFileList()[-1], gdal.GA_ReadOnly) - noDataValue = ds.GetRasterBand(1).GetNoDataValue() - print(f'grab NoDataValue for {layer}: {noDataValue:<5} and convert to 0.') + elif dsName == 'ionosphere': + data *= -1.0 #date2_date1 -> date1_date2 + f["unwrapPhase"][ii,:,:] = data - # sort the order of correction layer pairs based on date1_date2 with date1 < date2 - nPairs = dsCor.RasterCount - d12BandDict = {} - for ii in range(nPairs): - bnd = dsCor.GetRasterBand(ii+1) - d12 = bnd.GetMetadata(layer)["Dates"] - d12 = sorted(d12.split("_")) - d12 = f'{d12[0]}_{d12[1]}' - d12BandDict[d12] = ii+1 - d12List = sorted(d12BandDict.keys()) - print(f'number of {layer} pairs: {len(d12List)}') - - # box to gdal arguments - # link: https://gdal.org/python/osgeo.gdal.Band-class.html#ReadAsArray - if box is not None: - kwargs = dict( - xoff=box[0], - yoff=box[1], - win_xsize=box[2]-box[0], - win_ysize=box[3]-box[1]) - else: - kwargs = dict() - - if xstep * ystep > 1: - msg = f'apply {xstep} x {ystep} multilooking/downsampling via {mli_method} to {layer}' - print(msg) - print(f'writing data to HDF5 file {outfile} with a mode ...') - with h5py.File(outfile, "a") as f: - - prog_bar = ptime.progressBar(maxValue=nPairs) - for ii in range(nPairs): - d12 = d12List[ii] - bndIdx = d12BandDict[d12] - prog_bar.update(ii+1, suffix=f'{d12} {ii+1}/{nPairs}') - - f["date"][ii,0] = d12.split("_")[0].encode("utf-8") - f["date"][ii,1] = d12.split("_")[1].encode("utf-8") - f["dropIfgram"][ii] = True - - bnd = dsCor.GetRasterBand(bndIdx) - data = bnd.ReadAsArray(**kwargs) - data = multilook_data(data, ystep, xstep, method=mli_method) - data[data == noDataValue] = 0 #assign pixel with no-data to 0 - data[np.isnan(data)] = 0 #assign nan pixel to 0 - data *= -1.0 #date2_date1 -> date1_date2 - f["unwrapPhase"][ii,:,:] = data - # NOTE: ifgramStack accepts only 'unwrapPhase', 'rangeOffset', 'azimuthOffset' - # if changed to different layer name, changes need be be done in - # stack.py line: 782 - # reference_point.py line: 45 - # ifgram_inversion.py line 651 - # add option to recognize the layer - - bperp = float(bnd.GetMetadata(str(layer))["perpendicularBaseline"]) - bperp *= -1.0 #date2_date1 -> date1_date2 - f["bperp"][ii] = bperp + bperp = float(bnd.GetMetadata("ionosphere")["perpendicularBaseline"]) + bperp *= -1.0 #date2_date1 -> date1_date2 + f["bperp"][ii] = bperp prog_bar.close() # add MODIFICATION_TIME metadata to each 3D dataset - for dsName in ['unwrapPhase']: + for dsName in stackFiles.keys(): + dsName = 'unwrapPhase' if dsName == 'ionosphere' else dsName f[dsName].attrs['MODIFICATION_TIME'] = str(time.time()) print(f'finished writing to HD5 file: {outfile}\n') - ds = None - dsCor = None - + dsUnw = None + dsCoh = None + dsComp = None + dsAmp = None return outfile +# OPTIONAL - ARIA model-based corrections troposphereTotal, solidearthtides def write_model_stack(outfile, corrStack, box=None, xstep=1, ystep=1, mli_method='nearest'): """Write SET and TropsphericDelay corrections to HDF5 file from stack VRT files @@ -528,11 +434,17 @@ def write_model_stack(outfile, corrStack, box=None, """ print('-'*50) + + # determine field length for printing max_digit = len(os.path.basename(str(corrStack))) if corrStack is not None: print('open {f:<{w}} with gdal ...'.format(f=os.path.basename(corrStack), w=max_digit)) + # check all files exist + if not os.path.exists(corrStack): + raise Exception("%s does not exist" % corrStack) + # open raster dsCor = gdal.Open(corrStack, gdal.GA_ReadOnly) # extract NoDataValue (from the last date.vrt file for example) @@ -613,88 +525,6 @@ def get_number_of_epochs(vrtfile): return ds.RasterCount -def invert_diff_corrections(input_filename, output_filename, dataset, - cluster = None, num_workers = '4', waterMask=None, - maskDataset = None, maskThreshold = 0.4): - ''' - Invert differential ARIA correction layers to get correction for SAR acquistion dates - - the inversion gives corrections relative to the REF_DATE (typically the first acquisition) - - NOTE: 1.the inversion network needs to be connected, otherwise inversion will give wrong estimates - for the isolated clusters - 2. each SAR date needs to have min two datasets with that date to have min_degree of freedom >=2 - otherwise inversion will give wrong estimates for those dates - ''' - - # create inps dummy - class dummy(): - pass - - # PREPARE INPUT OBJECT - inps = dummy() - # input - inps.ifgramStackFile = input_filename - inps.obsDatasetName = 'unwrapPhase' #only this available for use at the moment - inps.skip_ref = False - - # solver - inps.minNormVelocity = False - inps.minRedundancy = 1.0 - # Note: minRedun set to 2.0 gives wrong estim, and dask has some errors - inps.weightFunc = 'no' - inps.calcCov = False - - # mask - inps.waterMaskFile = waterMask - inps.maskDataset = maskDataset - inps.maskThreshold = maskThreshold - - # cluster - Expose / leave None for now - inps.cluster = cluster - inps.maxMemory = 2 - inps.numWorker = num_workers - inps.config = 'local' #not sure what to put here - - # outputs - #Avoid setting 'no' for invQualityFile - #self.invQualityFile = 'modelTempCoh.h5' # not needed, maybe leave it to avoid dask outout issues - inps.invQualityFile = 'modelTempCoh.h' - inps.numInvFile = 'numInvModel.h5' - inps.tsFile = output_filename - - ### INVERSION - from mintpy.utils import readfile - from mintpy.cli import reference_point - from mintpy.ifgram_inversion import run_ifgram_inversion - - # remove if exists - try: - os.remove(inps.tsFile) - print('Delete existing file') - except FileNotFoundError: - print("File is not present in the system.") - - # 1. get reference pixel - - reference_point.main([inps.ifgramStackFile]) - - # 2. invert differential model obs - run_ifgram_inversion(inps) - - # 3. compensate for range2phase conversion as not needed for models - if dataset.startswith(('tropo', 'set')): - print('Return back units to original') - data, metadata = readfile.read(inps.tsFile, datasetName='timeseries') - phase2range = -1 * float(metadata['WAVELENGTH']) / (4.*np.pi) - - #Replace values - with h5py.File(inps.tsFile, 'r+') as f: - f['timeseries'][:] = data / phase2range - - # 4. clean up uncessary files - os.remove(inps.invQualityFile) - os.remove(inps.numInvFile) - def get_correction_layer(correction_filename): ds = gdal.Open(correction_filename, gdal.GA_ReadOnly) @@ -761,13 +591,12 @@ def load_aria(inps): compression=inps.compression, ) - # write data to h5 file in disk write_ifgram_stack( inps.outfile[0], - unwStack=inps.unwFile, - cohStack=inps.corFile, - connCompStack=inps.connCompFile, - ampStack=inps.magFile, + stackFiles={'unwrapPhase': inps.unwFile, + 'coherence': inps.corFile, + 'connectComponent': inps.connCompFile, + 'magnitude': inps.magFile}, box=box, xstep=inps.xstep, ystep=inps.ystep, @@ -819,39 +648,33 @@ def load_aria(inps): 'dropIfgram' : (np.bool_, (num_pair,)), 'bperp' : (np.float32, (num_pair,)), 'unwrapPhase' : (np.float32, (num_pair, length, width)), + "coherence" : (np.float32, (num_pair, length, width)), } meta['FILE_TYPE'] = 'ifgramStack' layer_name, layer_type = get_correction_layer(inps.ionoFile) if run_or_skip(inps, ds_name_dict, out_file=inps.outfile[0]) == 'run': + outname = f'{out_dir}/ionStack.h5' + writefile.layout_hdf5( - f'{out_dir}/d{layer_name}.h5', + outname, ds_name_dict, metadata=meta, compression=inps.compression, ) - + # write data to disk - write_iono_stack( - f'{out_dir}/d{layer_name}.h5', - ionStack=inps.ionoFile, + write_ifgram_stack( + outname, + stackFiles={'ionosphere': inps.ionoFile, + 'coherence': inps.corFile,}, box=box, xstep=inps.xstep, ystep=inps.ystep, - ) + mli_method=inps.method, + ) - # invert layer to get ionosphere phase delay - # on SAR acquistion dates - invert_diff_corrections(f'{out_dir}/d{layer_name}.h5', - f'{out_dir}/{layer_name}.h5', - layer_type, - cluster=inps.cluster, - num_workers=inps.num_workers, - waterMask = None, - maskDataset = None, - maskThreshold = 0.4) - # 3.2 - model based corrections: SolidEarthTides and Troposphere # Loop through other correction layers also provided as epochs # handle multiple tropo stacks (if specified) From 4d7d2bed1732399fa5882a15cfca3729d4a05748 Mon Sep 17 00:00:00 2001 From: ehavazli Date: Wed, 30 Oct 2024 11:52:14 -0700 Subject: [PATCH 08/13] add ARIA-tools example --- src/mintpy/cli/prep_aria.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mintpy/cli/prep_aria.py b/src/mintpy/cli/prep_aria.py index c0ae14bd7..05cf1635e 100755 --- a/src/mintpy/cli/prep_aria.py +++ b/src/mintpy/cli/prep_aria.py @@ -52,9 +52,9 @@ prep_aria.py -s ../stack/ -d ../DEM/SRTM_3arcsec.dem -i '../incidenceAngle/*.vrt' --set '../stack/setStack.vrt' --tropo '../stack/troposphereTotal/HRRR_stack.vrt' --iono '../stack/ionStack.vrt' # download / extract / prepare inteferograms stack from ARIA using ARIA-tools: - # reference: https://github.com/aria-tools/ARIA-tools - ariaDownload.py -b '37.25 38.1 -122.6 -121.75' --track 42 - ariaTSsetup.py -f 'products/*.nc' -b '37.25 38.1 -122.6 -121.75' --mask Download --num_threads 4 --verbose + ariaDownload.py --track 71 --bbox "34.21 34.31 -118.60 -118.43" --start 20240101 --end 20240301 + ariaTSsetup.py -f "products/*.nc" --bbox "34.21 34.31 -118.60 -118.43" --mask Download --num_threads 4 --layers "solidEarthTide, troposphereTotal, ionosphere" --tropo_models HRRR + prep_aria.py -s stack -d DEM/glo_90.dem -i incidenceAngle/*.vrt -a azimuthAngle/*.vrt -w mask/esa_world_cover_2021.msk --set stack/setStack.vrt --tropo stack/troposphereTotal/HRRRStack.vrt --iono stack/ionoStack.vrt """ def create_parser(subparsers=None): From 512fb7b20faa122f2e02fd84fef34ba143eb2eb7 Mon Sep 17 00:00:00 2001 From: ehavazli Date: Wed, 30 Oct 2024 14:45:35 -0700 Subject: [PATCH 09/13] change 'write_model_stack' to 'write_timeseries' and correct indentations --- src/mintpy/prep_aria.py | 187 ++++++++++++++++++++-------------------- 1 file changed, 93 insertions(+), 94 deletions(-) diff --git a/src/mintpy/prep_aria.py b/src/mintpy/prep_aria.py index d7b838786..6d7678999 100644 --- a/src/mintpy/prep_aria.py +++ b/src/mintpy/prep_aria.py @@ -329,7 +329,8 @@ def write_ifgram_stack(outfile, stackFiles, box=None, xstep=1, ystep=1, mli_meth noDataValues[dsName] = ds.GetRasterBand(1).GetNoDataValue() fileName = os.path.basename(stackFiles[dsName]) - print(f'grab NoDataValue for {fileName:<{max_digit}}: {noDataValues[dsName]:<5} and convert to 0.') + print(f'grab NoDataValue for {fileName:<{max_digit}}: ' + f'{noDataValues[dsName]:<5} and convert to 0.') ds = None # sort the order of interferograms based on date1_date2 with date1 < date2 @@ -375,10 +376,10 @@ def write_ifgram_stack(outfile, stackFiles, box=None, xstep=1, ystep=1, mli_meth bnd = dsStack.GetRasterBand(bndIdx) data = bnd.ReadAsArray(**kwargs) if xstep * ystep > 1: - mli_method_spec = mli_method if stackName not in \ + mli_method_spec = mli_method if dsName not in \ ['connCompStack'] else 'nearest' print(f'apply {xstep} x {ystep} multilooking/downsampling via ' - f'{mli_method_spec} to: {stackName}') + f'{mli_method_spec} to: {dsName}') data = multilook_data(data, ystep, xstep, method=mli_method_spec) data[data == noDataValues[dsName]] = 0 #assign pixel with no-data to 0 @@ -423,101 +424,101 @@ def write_ifgram_stack(outfile, stackFiles, box=None, xstep=1, ystep=1, mli_meth # OPTIONAL - ARIA model-based corrections troposphereTotal, solidearthtides -def write_model_stack(outfile, corrStack, box=None, +def write_timeseries(outfile, corrStack, box=None, xstep=1, ystep=1, mli_method='nearest'): - """Write SET and TropsphericDelay corrections to HDF5 file from stack VRT files - Correction layers are stored for each SAR acquisition date + """Write SET and TropsphericDelay corrections to HDF5 file from stack VRT files + Correction layers are stored for each SAR acquisition date - ARIA_GUNW_NC_PATH: - troposhereTotal : models GMAO, HRRR, HRES, ERA5 '/science/grids/corrections/external/troposphere/' - solidEarthTides '/science/grids/corrections/derived/solidearthtides/' - """ + ARIA_GUNW_NC_PATH: + troposhereTotal : models GMAO, HRRR, HRES, ERA5 '/science/grids/corrections/external/troposphere/' + solidEarthTides '/science/grids/corrections/derived/solidearthtides/' + """ - print('-'*50) + print('-'*50) - # determine field length for printing - max_digit = len(os.path.basename(str(corrStack))) + # determine field length for printing + max_digit = len(os.path.basename(str(corrStack))) - if corrStack is not None: - print('open {f:<{w}} with gdal ...'.format(f=os.path.basename(corrStack), w=max_digit)) + if corrStack is not None: + print('open {f:<{w}} with gdal ...'.format(f=os.path.basename(corrStack), w=max_digit)) - # check all files exist - if not os.path.exists(corrStack): - raise Exception("%s does not exist" % corrStack) + # check all files exist + if not os.path.exists(corrStack): + raise Exception("%s does not exist" % corrStack) - # open raster - dsCor = gdal.Open(corrStack, gdal.GA_ReadOnly) - # extract NoDataValue (from the last date.vrt file for example) - ds = gdal.Open(dsCor.GetFileList()[-1], gdal.GA_ReadOnly) - noDataValue = ds.GetRasterBand(1).GetNoDataValue() - ds = None + # open raster + dsCor = gdal.Open(corrStack, gdal.GA_ReadOnly) + # extract NoDataValue (from the last date.vrt file for example) + ds = gdal.Open(dsCor.GetFileList()[-1], gdal.GA_ReadOnly) + noDataValue = ds.GetRasterBand(1).GetNoDataValue() + ds = None - # get the layer name (for tropo this will get the model name) - layer = dsCor.GetRasterBand(1).GetMetadataDomainList()[0] + # get the layer name (for tropo this will get the model name) + layer = dsCor.GetRasterBand(1).GetMetadataDomainList()[0] + + # Get the wavelength. need to convert radians to meters + wavelength = np.float64(dsCor.GetRasterBand(1).GetMetadata(layer)["Wavelength (m)"]) + phase2range = -wavelength / (4.*np.pi) + + # get model dates and time + nDate = dsCor.RasterCount + dateDict = {} + sensingDict = {} + for ii in range(nDate): + bnd = dsCor.GetRasterBand(ii+1) + date = bnd.GetMetadata(layer)["Dates"] + utc = dt.datetime.strptime(date + ',' + \ + bnd.GetMetadata(layer)["UTCTime (HH:MM:SS.ss)"], + "%Y%m%d,%H:%M:%S.%f") + dateDict[date] = ii+1 + sensingDict[ii+1] = str(utc) + dateList = sorted(dateDict.keys()) + print(f'number of {layer} datasets: {len(dateList)}') + + # box to gdal arguments + # link: https://gdal.org/python/osgeo.gdal.Band-class.html#ReadAsArray + if box is not None: + kwargs = dict( + xoff=box[0], + yoff=box[1], + win_xsize=box[2]-box[0], + win_ysize=box[3]-box[1]) + else: + kwargs = dict() - # Get the wavelength. need to convert radians to meters - wavelength = np.float64(dsCor.GetRasterBand(1).GetMetadata(layer)["Wavelength (m)"]) - phase2range = wavelength / (4.*np.pi) + if xstep * ystep > 1: + msg = f'apply {xstep} x {ystep} multilooking/downsampling via {mli_method} to {layer}' + print(msg) - # get model dates and time - nDate = dsCor.RasterCount - dateDict = {} - sensingDict = {} + print(f'writing data to HDF5 file {outfile} with a mode ...') + with h5py.File(outfile, "a") as f: + prog_bar = ptime.progressBar(maxValue=nDate) for ii in range(nDate): - bnd = dsCor.GetRasterBand(ii+1) - date = bnd.GetMetadata(layer)["Dates"] - utc = dt.datetime.strptime(date + ',' + \ - bnd.GetMetadata(layer)["UTCTime (HH:MM:SS.ss)"], - "%Y%m%d,%H:%M:%S.%f") - dateDict[date] = ii+1 - sensingDict[ii+1] = str(utc) - dateList = sorted(dateDict.keys()) - print(f'number of {layer} datasets: {len(dateList)}') - - # box to gdal arguments - # link: https://gdal.org/python/osgeo.gdal.Band-class.html#ReadAsArray - if box is not None: - kwargs = dict( - xoff=box[0], - yoff=box[1], - win_xsize=box[2]-box[0], - win_ysize=box[3]-box[1]) - else: - kwargs = dict() - - if xstep * ystep > 1: - msg = f'apply {xstep} x {ystep} multilooking/downsampling via {mli_method} to {layer}' - print(msg) - - print(f'writing data to HDF5 file {outfile} with a mode ...') - with h5py.File(outfile, "a") as f: - prog_bar = ptime.progressBar(maxValue=nDate) - for ii in range(nDate): - date = dateList[ii] - bndIdx = dateDict[date] - utc = sensingDict[bndIdx] - prog_bar.update(ii+1, suffix=f'{date} {ii+1}/{nDate}') - - f["date"][ii] = date.encode("utf-8") - f["sensingMid"][ii] = utc.encode("utf-8") - - bnd = dsCor.GetRasterBand(bndIdx) - data = bnd.ReadAsArray(**kwargs) - data = multilook_data(data, ystep, xstep, method=mli_method) - data[data == noDataValue] = 0 #assign pixel with no-data to 0 - data[np.isnan(data)] = 0 #assign nan pixel to 0 - f["timeseries"][ii,:,:] = data * phase2range - - prog_bar.close() - - # add MODIFICATION_TIME metadata to each 3D dataset - for dsName in ['timeseries']: - f[dsName].attrs['MODIFICATION_TIME'] = str(time.time()) - - print(f'finished writing to HD5 file: {outfile}\n') - dsCor = None - - return outfile + date = dateList[ii] + bndIdx = dateDict[date] + utc = sensingDict[bndIdx] + prog_bar.update(ii+1, suffix=f'{date} {ii+1}/{nDate}') + + f["date"][ii] = date.encode("utf-8") + f["sensingMid"][ii] = utc.encode("utf-8") + + bnd = dsCor.GetRasterBand(bndIdx) + data = bnd.ReadAsArray(**kwargs) + data = multilook_data(data, ystep, xstep, method=mli_method) + data[data == noDataValue] = 0 #assign pixel with no-data to 0 + data[np.isnan(data)] = 0 #assign nan pixel to 0 + f["timeseries"][ii,:,:] = data * phase2range + + prog_bar.close() + + # add MODIFICATION_TIME metadata to each 3D dataset + for dsName in ['timeseries']: + f[dsName].attrs['MODIFICATION_TIME'] = str(time.time()) + + print(f'finished writing to HD5 file: {outfile}\n') + dsCor = None + + return outfile def get_number_of_epochs(vrtfile): @@ -537,10 +538,10 @@ def get_correction_layer(correction_filename): else: # ionosphere, solid earth tides layer_type = layer_name - + #close ds = None - + return layer_name, layer_type #################################################################################### @@ -679,7 +680,7 @@ def load_aria(inps): # Loop through other correction layers also provided as epochs # handle multiple tropo stacks (if specified) if inps.tropoFile is None: - inps.tropoFile = [None] + inps.tropoFile = [None] correction_layers = inps.tropoFile + [inps.setFile] for layer in correction_layers: if layer: @@ -710,9 +711,9 @@ def load_aria(inps): ) # write data to disk - write_model_stack( + write_timeseries( f'{out_dir}/{layer_name}_ARIA.h5', - corrStack=layer, + corrStack=layer, box=box, xstep=inps.xstep, ystep=inps.ystep, @@ -722,5 +723,3 @@ def load_aria(inps): # used time m, s = divmod(time.time() - start_time, 60) print(f'time used: {m:02.0f} mins {s:02.1f} secs.') - - From 36d606a4dcd9af610f9e154029e82e8f2ccae01e Mon Sep 17 00:00:00 2001 From: ehavazli Date: Wed, 30 Oct 2024 14:59:58 -0700 Subject: [PATCH 10/13] fix pre-commit hooks --- src/mintpy/prep_aria.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/mintpy/prep_aria.py b/src/mintpy/prep_aria.py index 6d7678999..acfb0f2b6 100644 --- a/src/mintpy/prep_aria.py +++ b/src/mintpy/prep_aria.py @@ -5,12 +5,12 @@ ############################################################ +import datetime as dt import os import time import h5py import numpy as np -import datetime as dt try: from osgeo import gdal @@ -425,7 +425,7 @@ def write_ifgram_stack(outfile, stackFiles, box=None, xstep=1, ystep=1, mli_meth # OPTIONAL - ARIA model-based corrections troposphereTotal, solidearthtides def write_timeseries(outfile, corrStack, box=None, - xstep=1, ystep=1, mli_method='nearest'): + xstep=1, ystep=1, mli_method='nearest'): """Write SET and TropsphericDelay corrections to HDF5 file from stack VRT files Correction layers are stored for each SAR acquisition date @@ -496,7 +496,7 @@ def write_timeseries(outfile, corrStack, box=None, for ii in range(nDate): date = dateList[ii] bndIdx = dateDict[date] - utc = sensingDict[bndIdx] + utc = sensingDict[bndIdx] prog_bar.update(ii+1, suffix=f'{date} {ii+1}/{nDate}') f["date"][ii] = date.encode("utf-8") @@ -676,7 +676,7 @@ def load_aria(inps): mli_method=inps.method, ) - # 3.2 - model based corrections: SolidEarthTides and Troposphere + # 3.2 - model based corrections: SolidEarthTides and Troposphere # Loop through other correction layers also provided as epochs # handle multiple tropo stacks (if specified) if inps.tropoFile is None: @@ -717,7 +717,7 @@ def load_aria(inps): box=box, xstep=inps.xstep, ystep=inps.ystep, - ) + ) print('-'*50) # used time From cde446bf726fc58121b01f36ff3ac8345fba730b Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Sun, 3 Nov 2024 12:23:39 +0800 Subject: [PATCH 11/13] Update prep_aria.py --- src/mintpy/cli/prep_aria.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mintpy/cli/prep_aria.py b/src/mintpy/cli/prep_aria.py index 05cf1635e..cd28df1a8 100755 --- a/src/mintpy/cli/prep_aria.py +++ b/src/mintpy/cli/prep_aria.py @@ -50,7 +50,7 @@ prep_aria.py -s ../stack/ -d ../DEM/SRTM_3arcsec.dem -i '../incidenceAngle/*.vrt' prep_aria.py -s ../stack/ -d ../DEM/SRTM_3arcsec.dem -i '../incidenceAngle/*.vrt' -a '../azimuthAngle/*.vrt' -w ../mask/watermask.msk prep_aria.py -s ../stack/ -d ../DEM/SRTM_3arcsec.dem -i '../incidenceAngle/*.vrt' --set '../stack/setStack.vrt' --tropo '../stack/troposphereTotal/HRRR_stack.vrt' --iono '../stack/ionStack.vrt' - + # download / extract / prepare inteferograms stack from ARIA using ARIA-tools: ariaDownload.py --track 71 --bbox "34.21 34.31 -118.60 -118.43" --start 20240101 --end 20240301 ariaTSsetup.py -f "products/*.nc" --bbox "34.21 34.31 -118.60 -118.43" --mask Download --num_threads 4 --layers "solidEarthTide, troposphereTotal, ionosphere" --tropo_models HRRR @@ -107,7 +107,7 @@ def create_parser(subparsers=None): help='Name of the azimuth angle file.') geom.add_argument('-w','--water-mask', dest='waterMaskFile', type=str, help='Name of the water mask file') - + # correction layers: troposphereTotal, ionosphere, solidEarthTides corr = parser.add_argument_group('corrections') corr.add_argument('-ct', '--tropo', dest='tropoFile', type=str, From 87151529789dd8881d1fa5e1b9383566a7e1a3ed Mon Sep 17 00:00:00 2001 From: ehavazli Date: Tue, 5 Nov 2024 11:43:10 -0800 Subject: [PATCH 12/13] remove unused variable --- src/mintpy/prep_aria.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mintpy/prep_aria.py b/src/mintpy/prep_aria.py index acfb0f2b6..e0be72554 100644 --- a/src/mintpy/prep_aria.py +++ b/src/mintpy/prep_aria.py @@ -653,7 +653,7 @@ def load_aria(inps): } meta['FILE_TYPE'] = 'ifgramStack' - layer_name, layer_type = get_correction_layer(inps.ionoFile) + layer_name, _ = get_correction_layer(inps.ionoFile) if run_or_skip(inps, ds_name_dict, out_file=inps.outfile[0]) == 'run': outname = f'{out_dir}/ionStack.h5' @@ -685,7 +685,7 @@ def load_aria(inps): for layer in correction_layers: if layer: # get name and type - layer_name, layer_type = get_correction_layer(layer) + layer_name, _ = get_correction_layer(layer) num_dates = get_number_of_epochs(layer) meta['FILE_TYPE'] = 'timeseries' From 981d3ce608964df8d71b6c2fdbf03a7c5de55c03 Mon Sep 17 00:00:00 2001 From: ehavazli Date: Thu, 14 Nov 2024 10:53:06 -0800 Subject: [PATCH 13/13] remove unused arguments and comments --- src/mintpy/cli/prep_aria.py | 4 ---- src/mintpy/prep_aria.py | 1 - 2 files changed, 5 deletions(-) diff --git a/src/mintpy/cli/prep_aria.py b/src/mintpy/cli/prep_aria.py index cd28df1a8..d55a16db7 100755 --- a/src/mintpy/cli/prep_aria.py +++ b/src/mintpy/cli/prep_aria.py @@ -116,10 +116,6 @@ def create_parser(subparsers=None): help='Name of the Ionosphere Delay stack file', default=None) corr.add_argument('-cs', '--set', dest='setFile', type=str, help='Name of the Solid Earth Tides stack file', default=None) - corr.add_argument('--cluster', dest='cluster', type=str, choices={'local', 'pbs', None}, - help='Parallelize inversion of correction layers w Dask', default=None) - corr.add_argument('-n', '--num-workers', dest='num_workers', type=str, - help='Dask number of workers', default='2') return parser diff --git a/src/mintpy/prep_aria.py b/src/mintpy/prep_aria.py index e0be72554..b21fe6824 100644 --- a/src/mintpy/prep_aria.py +++ b/src/mintpy/prep_aria.py @@ -641,7 +641,6 @@ def load_aria(inps): ########## output file 3 - correction layers # 3.1 - ionosphere - # Invert Iono stack and write out cube if inps.ionoFile: # define correction dataset structure for ifgramStack ds_name_dict = {