Revision f24076bbf4a4a4af8e90100f344e10fb237e9b3b authored by Zhang Yunjun on 08 November 2022, 07:17:03 UTC, committed by GitHub on 08 November 2022, 07:17:03 UTC
+ .circleci/config.yml update: - use $BASH_ENV to set and share environment variable (PATH) among multiple run steps - rename to workflow/job to "unit-n-workflow-tests" + cli/save_gmt.py: fix a typo in the module import + cli/load_data: check `-t smallbaselineApp.cfg` existence and print out error msg + tsview: show the reference index/date info on the slider as the title
1 parent 64ffdef
smallbaselineApp.py
############################################################
# Program is part of MintPy #
# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi #
# Author: Zhang Yunjun, Heresh Fattahi, Jul 2013 #
############################################################
import glob
import os
import shutil
import time
import numpy as np
import mintpy
from mintpy.objects import RAMP_LIST, cluster, sensor
from mintpy.utils import readfile, utils as ut, writefile
# comment out the workflow import in favor of individual imports for speed and robustness
# import mintpy.workflow # dynamic import of modules for smallbaselineApp
##########################################################################
def get_the_latest_default_template_file(work_dir):
"""Get the latest version of default template file.
If an obsolete file exists in the working directory, the existing option values are kept.
"""
# get the latest and current versions of the default template files
lfile = os.path.join(os.path.dirname(mintpy.__file__), 'defaults/smallbaselineApp.cfg')
cfile = os.path.join(work_dir, 'smallbaselineApp.cfg')
if not os.path.isfile(cfile):
print(f'copy default template file {lfile} to work directory')
shutil.copy2(lfile, work_dir)
else:
#cfile is obsolete if any key is missing
ldict = readfile.read_template(lfile)
cdict = readfile.read_template(cfile)
if any([key not in cdict.keys() for key in ldict.keys()]):
print('obsolete default template detected, update to the latest version.')
shutil.copy2(lfile, work_dir)
#keep the existing option value from obsolete template file
ut.update_template_file(cfile, cdict)
return cfile
##########################################################################
class TimeSeriesAnalysis:
""" Routine processing workflow for time series analysis of small baseline InSAR stacks
"""
def __init__(self, customTemplateFile=None, workDir=None):
self.customTemplateFile = customTemplateFile
self.workDir = os.path.abspath(workDir)
self.cwd = os.path.abspath(os.getcwd())
def open(self):
"""The starting point of the workflow. It runs everytime.
It 1) grab project name if given
2) go to work directory
3) get and read template(s) options
"""
#1. Get projectName
self.projectName = None
if self.customTemplateFile:
self.projectName = os.path.splitext(os.path.basename(self.customTemplateFile))[0]
print('Project name:', self.projectName)
#2. Go to work directory
os.makedirs(self.workDir, exist_ok=True)
os.chdir(self.workDir)
print('Go to work directory:', self.workDir)
#3. Read templates
self.templateFile = get_the_latest_default_template_file(self.workDir)
self._read_template()
def _read_template(self):
# 1) update default template
self.customTemplate = None
if self.customTemplateFile:
# customTemplateFile --> customTemplate
print('read custom template file:', self.customTemplateFile)
cdict = readfile.read_template(self.customTemplateFile)
# correct some loose type errors
standardValues = {
'def':'auto', 'default':'auto',
'y':'yes', 'on':'yes', 'true':'yes',
'n':'no', 'off':'no', 'false':'no',
}
for key, value in cdict.items():
if value in standardValues.keys():
cdict[key] = standardValues[value]
for key in ['mintpy.deramp', 'mintpy.troposphericDelay.method']:
if key in cdict.keys():
cdict[key] = cdict[key].lower().replace('-', '_')
if 'processor' in cdict.keys():
cdict['mintpy.load.processor'] = cdict['processor']
# these metadata are used in load_data.py only, not needed afterwards
# (in order to manually add extra offset when the lookup table is shifted)
# (seen in ROI_PAC product sometimes)
for key in ['SUBSET_XMIN', 'SUBSET_YMIN']:
if key in cdict.keys():
cdict.pop(key)
self.customTemplate = dict(cdict)
# customTemplate --> templateFile
print('update default template based on input custom template')
self.templateFile = ut.update_template_file(self.templateFile, self.customTemplate)
# 2) backup custome/default template file in inputs/pic folder
flen = len(os.path.basename(self.templateFile))
if self.customTemplateFile:
flen = max(flen, len(os.path.basename(self.customTemplateFile)))
for backup_dirname in ['inputs', 'pic']:
backup_dir = os.path.join(self.workDir, backup_dirname)
# create directory
os.makedirs(backup_dir, exist_ok=True)
# back up to the directory
tfiles = [x for x in [self.customTemplateFile, self.templateFile] if x]
for tfile in tfiles:
out_file = os.path.join(backup_dir, os.path.basename(tfile))
if ut.run_or_skip(out_file, in_file=tfile, readable=False, print_msg=False) == 'run':
shutil.copy2(tfile, backup_dir)
print('copy {f:<{l}} to {d:<8} directory for backup.'.format(
f=os.path.basename(tfile), l=flen, d=os.path.basename(backup_dir)))
# 3) read default template file
print('read default template file:', self.templateFile)
self.template = readfile.read_template(self.templateFile)
self.template = ut.check_template_auto_value(self.template)
# correct some loose setup conflicts
if self.template['mintpy.geocode'] is False:
for key in ['mintpy.save.hdfEos5', 'mintpy.save.kmz']:
if self.template[key] is True:
self.template['mintpy.geocode'] = True
print(f'Turn ON mintpy.geocode in order to run {key}.')
break
def run_load_data(self, step_name):
"""Load InSAR stacks into HDF5 files in ./inputs folder.
It 1) copy auxiliary files into work directory (for Unvi of Miami only)
2) load all interferograms stack files into mintpy/inputs directory.
3) check loading result
4) add custom metadata (optional, for HDF-EOS5 format only)
"""
# 1) copy aux files (optional)
self._copy_aux_file()
# 2) loading data
# compose list of input arguments
# instead of using command line then split
# to support path with whitespace
iargs = ['--template', self.templateFile]
if self.customTemplateFile:
iargs += [self.customTemplateFile]
if self.projectName:
iargs += ['--project', self.projectName]
# run command line
print('\nload_data.py', ' '.join(iargs))
import mintpy.cli.load_data
mintpy.cli.load_data.main(iargs)
# come back to working directory
os.chdir(self.workDir)
# 3) check loading result
stack_file, geom_file, _, ion_file = ut.check_loaded_dataset(self.workDir, print_msg=True)[:4]
# 4) add custom metadata (optional)
if self.customTemplateFile:
# use ut.add_attribute() instead of add_attribute.py because of
# better control of special metadata, such as SUBSET_X/YMIN
msg = f'updating metadata based on custom template file {os.path.basename(self.customTemplateFile)}'
for fname in [stack_file, ion_file, geom_file]:
if fname:
print(f'{msg} for file: {os.path.basename(fname)}')
ut.add_attribute(fname, self.customTemplate)
def _copy_aux_file(self):
# for Univ of Miami
if os.getenv('SCRATCHDIR') and self.projectName:
proj_dir = os.path.join(os.getenv('SCRATCHDIR'), self.projectName)
flist = ['PROCESS/unavco_attributes.txt',
'PROCESS/bl_list.txt',
'SLC/summary*slc.jpg']
flist = ut.get_file_list([os.path.join(proj_dir, i) for i in flist], abspath=True)
for fname in flist:
if ut.run_or_skip(os.path.basename(fname), in_file=fname, readable=False) == 'run':
shutil.copy2(fname, self.workDir)
print(f'copy {os.path.basename(fname)} to work directory')
def run_network_modification(self, step_name):
"""Modify network of interferograms before the network inversion."""
# check the existence of ifgramStack.h5
stack_file, geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[:2]
coh_txt = os.path.join(self.workDir, 'coherenceSpatialAvg.txt')
net_fig = [os.path.join(self.workDir, i, 'network.pdf') for i in ['', 'pic']]
net_fig = [x for x in net_fig if os.path.isfile(x)]
# 1) output waterMask.h5 to simplify the detection/use of waterMask
water_mask_file = os.path.join(self.workDir, 'waterMask.h5')
if 'waterMask' in readfile.get_dataset_list(geom_file):
print(f'generate {water_mask_file} from {geom_file} for conveniency')
if ut.run_or_skip(out_file=water_mask_file, in_file=geom_file) == 'run':
water_mask, atr = readfile.read(geom_file, datasetName='waterMask')
# ignore no-data pixels in geometry files
ds_name_list = readfile.get_dataset_list(geom_file)
for ds_name in ['latitude','longitude']:
if ds_name in ds_name_list:
print(f'set pixels with 0 in {ds_name} to 0 in waterMask')
ds = readfile.read(geom_file, datasetName=ds_name)[0]
water_mask[ds == 0] = 0
atr['FILE_TYPE'] = 'waterMask'
writefile.write(water_mask, out_file=water_mask_file, metadata=atr)
# 2) modify network
iargs = [stack_file, '-t', self.templateFile]
print('\nmodify_network.py', ' '.join(iargs))
import mintpy.cli.modify_network
mintpy.cli.modify_network.main(iargs)
# 3) plot network
iargs = [stack_file, '-t', self.templateFile, '--nodisplay']
dsNames = readfile.get_dataset_list(stack_file)
if any('phase' in i.lower() for i in dsNames):
iargs += ['-d', 'coherence', '-v', '0.2', '1.0']
elif any('offset' in i.lower() for i in dsNames):
iargs += ['-d', 'offsetSNR', '-v', '0', '20']
print('\nplot_network.py', ' '.join(iargs))
# run
if self.template['mintpy.plot']:
if ut.run_or_skip(out_file=net_fig,
in_file=[stack_file, coh_txt, self.templateFile],
readable=False) == 'run':
import mintpy.cli.plot_network
mintpy.cli.plot_network.main(iargs)
else:
print('mintpy.plot is turned OFF, skip plotting network.')
def generate_ifgram_aux_file(self):
"""Generate auxiliary files from ifgramStack file"""
stack_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[0]
dsNames = readfile.get_dataset_list(stack_file)
mask_file = os.path.join(self.workDir, 'maskConnComp.h5')
coh_file = os.path.join(self.workDir, 'avgSpatialCoh.h5')
snr_file = os.path.join(self.workDir, 'avgSpatialSNR.h5')
# 1) generate mask file from the common connected components
if any('phase' in i.lower() for i in dsNames):
iargs = [stack_file, '--nonzero', '-o', mask_file, '--update']
print('\ngenerate_mask.py', ' '.join(iargs))
import mintpy.cli.generate_mask
mintpy.cli.generate_mask.main(iargs)
# 2) generate average spatial coherence
if any('phase' in i.lower() for i in dsNames):
iargs = [stack_file, '--dataset', 'coherence', '-o', coh_file, '--update']
elif any('offset' in i.lower() for i in dsNames):
iargs = [stack_file, '--dataset', 'offsetSNR', '-o', snr_file, '--update']
print('\ntemporal_average.py', ' '.join(iargs))
import mintpy.cli.temporal_average
mintpy.cli.temporal_average.main(iargs)
def run_reference_point(self, step_name):
"""Select reference point.
It 1) generate mask file from common conn comp
2) generate average spatial coherence and its mask
3) add REF_X/Y and/or REF_LAT/LON attribute to stack file
"""
# 1-2) aux files: maskConnComp and avgSpatialCoh
self.generate_ifgram_aux_file()
# 3) add REF_X/Y(/LAT/LON) of the reference point
stack_file, _, lookup_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[:3]
coh_file = os.path.join(self.workDir, 'avgSpatialCoh.h5')
iargs = [stack_file, '-t', self.templateFile, '-c', coh_file]
if lookup_file is not None:
iargs += ['--lookup', lookup_file]
print('\nreference_point.py', ' '.join(iargs))
import mintpy.cli.reference_point
mintpy.cli.reference_point.main(iargs)
def run_quick_overview(self, step_name):
"""A quick overview on the interferogram stack for:
1) avgPhaseVelocity.h5: possible ground deformation through interferogram stacking
2) numTriNonzeroIntAmbiguity.h5: phase unwrapping errors
through the integer ambiguity of phase closure
"""
# check the existence of ifgramStack.h5
stack_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[0]
# 1) stack interferograms
pha_vel_file = os.path.join(self.workDir, 'avgPhaseVelocity.h5')
iargs = [stack_file, '--dataset', 'unwrapPhase', '-o', pha_vel_file, '--update']
print('\ntemporal_average.py', ' '.join(iargs))
import mintpy.cli.temporal_average
mintpy.cli.temporal_average.main(iargs)
# 2) calculate the number of interferogram triplets with non-zero integer ambiguity
water_mask_file = os.path.join(self.workDir, 'waterMask.h5')
iargs = [stack_file, '--water-mask', water_mask_file, '--action', 'calculate', '--update']
print('\nunwrap_error_phase_closure.py', ' '.join(iargs))
import mintpy.cli.unwrap_error_phase_closure
mintpy.cli.unwrap_error_phase_closure.main(iargs)
def run_unwrap_error_correction(self, step_name):
"""Correct phase-unwrapping errors"""
method = self.template['mintpy.unwrapError.method']
if not method:
print('phase-unwrapping error correction is OFF.')
return
# check required input files
stack_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[0]
mask_file = os.path.join(self.workDir, 'maskConnComp.h5')
iargs_bridge = [stack_file, '--template', self.templateFile, '--update']
iargs_closure = iargs_bridge + ['--cc-mask', mask_file]
import mintpy.cli.unwrap_error_bridging
import mintpy.cli.unwrap_error_phase_closure
if method == 'bridging':
print('\nunwrap_error_bridging.py', ' '.join(iargs_bridge))
mintpy.cli.unwrap_error_bridging.main(iargs_bridge)
elif method == 'phase_closure':
print('\nunwrap_error_phase_closure.py', ' '.join(iargs_closure))
mintpy.cli.unwrap_error_phase_closure.main(iargs_closure)
elif method == 'bridging+phase_closure':
iargs_bridge += ['-i', 'unwrapPhase',
'-o', 'unwrapPhase_bridging']
print('\nunwrap_error_bridging.py', ' '.join(iargs_bridge))
mintpy.cli.unwrap_error_bridging.main(iargs_bridge)
iargs_closure += ['-i', 'unwrapPhase_bridging',
'-o', 'unwrapPhase_bridging_phaseClosure']
print('\nunwrap_error_phase_closure.py', ' '.join(iargs_closure))
mintpy.cli.unwrap_error_phase_closure.main(iargs_closure)
else:
raise ValueError(f'un-recognized method: {method}')
def run_network_inversion(self, step_name):
"""Invert network of interferograms for raw phase time-series.
1) network inversion --> timeseries.h5, temporalCoherence.h5, numInvIfgram.h5
2) temporalCoherence.h5 --> maskTempCoh.h5
"""
# check the existence of ifgramStack.h5
stack_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[0]
# 1) invert ifgramStack for time-series
iargs = [stack_file, '-t', self.templateFile, '--update']
print('\nifgram_inversion.py', ' '.join(iargs))
import mintpy.cli.ifgram_inversion
mintpy.cli.ifgram_inversion.main(iargs)
# 2) get reliable pixel mask: maskTempCoh.h5
self.generate_temporal_coherence_mask()
def generate_temporal_coherence_mask(self):
"""Generate reliable pixel mask from temporal coherence"""
geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1]
tcoh_file = os.path.join(self.workDir, 'temporalCoherence.h5')
mask_file = os.path.join(self.workDir, 'maskTempCoh.h5')
tcoh_min = self.template['mintpy.networkInversion.minTempCoh']
# compose list of arguments
iargs = [tcoh_file, '-m', tcoh_min, '-o', mask_file]
# exclude pixels in shadow if shadowMask dataset is available
if (self.template['mintpy.networkInversion.shadowMask'] is True
and 'shadowMask' in readfile.get_dataset_list(geom_file)):
iargs += ['--base', geom_file, '--base-dataset', 'shadowMask', '--base-value', '1']
print('\ngenerate_mask.py', ' '.join(iargs))
# update mode: run only if:
# 1) output file exists and newer than input file, AND
# 2) all config keys are the same
config_keys = [f'mintpy.networkInversion.{i}' for i in ['minTempCoh','shadowMask']]
print('update mode: ON')
flag = 'skip'
if ut.run_or_skip(out_file=mask_file, in_file=tcoh_file, print_msg=False) == 'run':
flag = 'run'
else:
print(f'1) output file: {mask_file} already exists and newer than input file: {tcoh_file}')
atr = readfile.read_attribute(mask_file)
if any(str(self.template[i]) != atr.get(i, 'False') for i in config_keys):
flag = 'run'
print(f'2) NOT all key configuration parameters are the same: {config_keys}')
else:
print(f'2) all key configuration parameters are the same: {config_keys}')
print(f'run or skip: {flag}')
if flag == 'run':
import mintpy.cli.generate_mask
mintpy.cli.generate_mask.main(iargs)
# update config_keys
atr = {}
for key in config_keys:
atr[key] = self.template[key]
ut.add_attribute(mask_file, atr)
# check number of pixels selected in mask file for following analysis
num_pixel = np.sum(readfile.read(mask_file)[0] != 0.)
print(f'number of reliable pixels: {num_pixel}')
min_num_pixel = float(self.template['mintpy.networkInversion.minNumPixel'])
if num_pixel < min_num_pixel:
msg = f"Not enough reliable pixels (minimum of {int(min_num_pixel)}). "
msg += "Try the following:\n"
msg += "1) Check the reference pixel and make sure it's not in areas with unwrapping errors\n"
msg += "2) Check the network and make sure it's fully connected without subsets"
print(f'ERROR: {msg}')
# populate the pic folder to facilate the trouble shooting
self.plot_result(print_aux=False)
# terminate the program
raise RuntimeError(msg)
@staticmethod
def get_timeseries_filename(template, work_dir='./'):
"""Get input/output time-series filename for each step
Parameters: template : dict, content of smallbaselineApp.cfg
Returns: steps : dict of dicts, input/output filenames for each step
"""
work_dir = os.path.abspath(work_dir)
fname0 = os.path.join(work_dir, 'timeseries.h5')
fname1 = os.path.join(work_dir, 'timeseries.h5')
atr = readfile.read_attribute(fname0)
phase_correction_steps = ['correct_LOD',
'correct_SET',
'correct_troposphere',
'deramp',
'correct_topography']
# loop for all steps
steps = dict()
for sname in phase_correction_steps:
# fname0 == fname1 if no valid correction method is set.
fname0 = fname1
if sname == 'correct_LOD':
if atr['PLATFORM'].lower().startswith('env'):
fname1 = f'{os.path.splitext(fname0)[0]}_LOD.h5'
elif sname == 'correct_SET':
method = template['mintpy.solidEarthTides']
if method:
fname1 = f'{os.path.splitext(fname0)[0]}_SET.h5'
elif sname == 'correct_troposphere':
method = template['mintpy.troposphericDelay.method']
model = template['mintpy.troposphericDelay.weatherModel']
if method:
if method == 'height_correlation':
fname1 = f'{os.path.splitext(fname0)[0]}_tropHgt.h5'
elif method == 'gacos':
fname1 = f'{os.path.splitext(fname0)[0]}_GACOS.h5'
elif method == 'pyaps':
fname1 = f'{os.path.splitext(fname0)[0]}_{model}.h5'
else:
msg = f'un-recognized tropospheric correction method: {method}'
raise ValueError(msg)
elif sname == 'deramp':
method = template['mintpy.deramp']
if method:
if method in RAMP_LIST:
fname1 = f'{os.path.splitext(fname0)[0]}_ramp.h5'
else:
msg = f'un-recognized phase ramp type: {method}'
msg += f'\navailable ramp types:\n{RAMP_LIST}'
raise ValueError(msg)
elif sname == 'correct_topography':
method = template['mintpy.topographicResidual']
if method:
fname1 = f'{os.path.splitext(fname0)[0]}_demErr.h5'
step = dict()
step['input'] = fname0
step['output'] = fname1
steps[sname] = step
# step - reference_date
fnames = [steps[sname]['output'] for sname in phase_correction_steps]
fnames += [steps[sname]['input'] for sname in phase_correction_steps]
fnames = sorted(list(set(fnames)))
step = dict()
step['input'] = fnames
steps['reference_date'] = step
# step - velocity / geocode
step = dict()
step['input'] = steps['reference_date']['input'][-1]
steps['velocity'] = step
steps['geocode'] = step
# step - hdfeos5
if 'Y_FIRST' not in atr.keys():
step = dict()
fdir = os.path.dirname(steps['reference_date']['input'][-1])
fbase = os.path.basename(steps['reference_date']['input'][-1])
step['input'] = os.path.join(fdir, f'geo/geo_{fbase}')
steps['hdfeos5'] = step
return steps
def run_local_oscillator_drift_correction(self, step_name):
"""Correct local oscillator drift (LOD).
Automatically applied for Envisat data.
Automatically skipped for all the other data.
"""
geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1]
fnames = self.get_timeseries_filename(self.template, self.workDir)[step_name]
in_file = fnames['input']
out_file = fnames['output']
if in_file != out_file:
iargs = [in_file, geom_file, '-o', out_file]
print('\nlocal_oscilator_drift.py', ' '.join(iargs))
if ut.run_or_skip(out_file=out_file, in_file=in_file) == 'run':
import mintpy.cli.local_oscilator_drift
mintpy.cli.local_oscilator_drift.main(iargs)
else:
atr = readfile.read_attribute(in_file)
sat = atr.get('PLATFORM', None)
print(f'No local oscillator drift correction is needed for {sat}.')
def run_solid_earth_tides_correction(self, step_name):
"""Correct solid Earth tides (SET)."""
geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1]
fnames = self.get_timeseries_filename(self.template, self.workDir)[step_name]
in_file = fnames['input']
out_file = fnames['output']
if in_file != out_file:
iargs = [in_file, '-g', geom_file, '-o', out_file, '--update']
print('\nsolid_earth_tides.py', ' '.join(iargs))
if ut.run_or_skip(out_file=out_file, in_file=in_file) == 'run':
import mintpy.cli.solid_earth_tides
mintpy.cli.solid_earth_tides.main(iargs)
else:
print('No solid Earth tides correction.')
def run_tropospheric_delay_correction(self, step_name):
"""Correct tropospheric delays."""
geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1]
mask_file = os.path.join(self.workDir, 'maskTempCoh.h5')
fnames = self.get_timeseries_filename(self.template, self.workDir)[step_name]
in_file = fnames['input']
out_file = fnames['output']
if in_file != out_file:
poly_order = self.template['mintpy.troposphericDelay.polyOrder']
tropo_model = self.template['mintpy.troposphericDelay.weatherModel'].upper()
weather_dir = self.template['mintpy.troposphericDelay.weatherDir']
method = self.template['mintpy.troposphericDelay.method']
def get_dataset_size(fname):
atr = readfile.read_attribute(fname)
return (atr['LENGTH'], atr['WIDTH'])
# Phase/Elevation Ratio (Doin et al., 2009)
if method == 'height_correlation':
tropo_look = self.template['mintpy.troposphericDelay.looks']
tropo_min_cor = self.template['mintpy.troposphericDelay.minCorrelation']
iargs = [in_file,
'-g', geom_file,
'-p', poly_order,
'-m', mask_file,
'-o', out_file,
'-l', tropo_look,
'-t', tropo_min_cor]
print('tropospheric delay correction with height-correlation approach')
print('\ntropo_phase_elevation.py', ' '.join(iargs))
if ut.run_or_skip(out_file=out_file, in_file=in_file) == 'run':
import mintpy.cli.tropo_phase_elevation
mintpy.cli.tropo_phase_elevation.main(iargs)
# Weather re-analysis data with iterative tropospheric decomposition (GACOS)
# Yu et al. (2018, JGR)
elif method == 'gacos':
GACOS_dir = self.template['mintpy.troposphericDelay.gacosDir']
iargs = ['-f', in_file, '-g', geom_file, '-o', out_file, '--dir', GACOS_dir]
print('tropospheric delay correction with gacos approach')
print('\ntropo_gacos.py', ' '.join(iargs))
if ut.run_or_skip(out_file=out_file, in_file=in_file) == 'run':
import mintpy.cli.tropo_gacos
mintpy.cli.tropo_gacos.main(iargs)
# Weather Re-analysis Data (Jolivet et al., 2011;2014)
elif method == 'pyaps':
iargs = ['-f', in_file, '--model', tropo_model, '-g', geom_file, '-w', weather_dir]
print('Atmospheric correction using Weather Re-analysis dataset (PyAPS, Jolivet et al., 2011)')
print('Weather Re-analysis dataset:', tropo_model)
tropo_file = f'./inputs/{tropo_model}.h5'
if ut.run_or_skip(out_file=out_file, in_file=[in_file, tropo_file]) == 'run':
if os.path.isfile(tropo_file) and get_dataset_size(tropo_file) == get_dataset_size(in_file):
iargs = [in_file, tropo_file, '-o', out_file, '--force']
print('--------------------------------------------')
print(f'Use existed tropospheric delay file: {tropo_file}')
print('\ndiff.py', ' '.join(iargs))
import mintpy.cli.diff
mintpy.cli.diff.main(iargs)
else:
if tropo_model in ['ERA5']:
print('\ntropo_pyaps3.py', ' '.join(iargs))
import mintpy.cli.tropo_pyaps3
mintpy.cli.tropo_pyaps3.main(iargs)
elif tropo_model in ['MERRA', 'NARR']:
print('\ntropo_pyaps.py', ' '.join(iargs))
import mintpy.legacy.tropo_pyaps
mintpy.legacy.tropo_pyaps.main(iargs)
else:
raise ValueError(f'un-recognized dataset name: {tropo_model}.')
else:
print('No tropospheric delay correction.')
def run_phase_deramping(self, step_name):
"""Estimate and remove phase ramp from each acquisition."""
mask_file = self.template['mintpy.deramp.maskFile']
method = self.template['mintpy.deramp']
fnames = self.get_timeseries_filename(self.template, self.workDir)[step_name]
in_file = fnames['input']
out_file = fnames['output']
if in_file != out_file:
print(f'Remove for each acquisition a phase ramp: {method}')
iargs = [in_file, '-s', method, '-m', mask_file, '-o', out_file, '--update']
print('\nremove_ramp.py', ' '.join(iargs))
import mintpy.cli.remove_ramp
mintpy.cli.remove_ramp.main(iargs)
else:
print('No phase ramp removal.')
def run_topographic_residual_correction(self, step_name):
"""step - correct_topography
Topographic residual (DEM error) correction (optional).
"""
geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1]
fnames = self.get_timeseries_filename(self.template, self.workDir)[step_name]
in_file = fnames['input']
out_file = fnames['output']
if in_file != out_file:
iargs = [in_file, '-t', self.templateFile, '-o', out_file, '--update']
if self.template['mintpy.topographicResidual.pixelwiseGeometry']:
iargs += ['-g', geom_file]
print('\ndem_error.py', ' '.join(iargs))
import mintpy.cli.dem_error
mintpy.cli.dem_error.main(iargs)
else:
print('No topographic residual correction.')
def run_residual_phase_rms(self, step_name):
"""Noise evaluation based on the phase residual."""
res_file = 'timeseriesResidual.h5'
if os.path.isfile(res_file):
iargs = [res_file, '-t', self.templateFile]
print('\ntimeseries_rms.py', ' '.join(iargs))
import mintpy.cli.timeseries_rms
mintpy.cli.timeseries_rms.main(iargs)
else:
print('No residual phase file found! Skip residual RMS analysis.')
def run_reference_date(self, step_name):
"""Change reference date for all time-series files (optional)."""
if self.template['mintpy.reference.date']:
iargs = ['-t', self.templateFile]
in_files = self.get_timeseries_filename(self.template, self.workDir)[step_name]['input']
for in_file in in_files:
iargs += [in_file]
print('\nreference_date.py', ' '.join(iargs))
import mintpy.cli.reference_date
mintpy.cli.reference_date.main(iargs)
else:
print('No reference date change.')
def run_timeseries2velocity(self, step_name):
"""Estimate average velocity from displacement time-series"""
ts_file = self.get_timeseries_filename(self.template, self.workDir)[step_name]['input']
vel_file = os.path.join(self.workDir, 'velocity.h5')
iargs = [ts_file, '-t', self.templateFile, '-o', vel_file, '--update']
print('\ntimeseries2velocity.py', ' '.join(iargs))
import mintpy.cli.timeseries2velocity
mintpy.cli.timeseries2velocity.main(iargs)
# Velocity from estimated tropospheric delays
tropo_model = self.template['mintpy.troposphericDelay.weatherModel'].upper()
tropo_file = os.path.join(self.workDir, f'inputs/{tropo_model}.h5')
if os.path.isfile(tropo_file):
suffix = os.path.splitext(os.path.basename(tropo_file))[0]
tropo_vel_file = f'{os.path.splitext(vel_file)[0]}{suffix}.h5'
tropo_vel_file = os.path.join(self.workDir, tropo_vel_file)
iargs = [tropo_file, '-t', self.templateFile, '-o', tropo_vel_file, '--update']
# add reference info for a meaningful velocity to assess the impact of tropo delay on velocity
atr = readfile.read_attribute(vel_file)
iargs += ['--ref-date', atr['REF_DATE'], '--ref-yx', atr['REF_Y'], atr['REF_X']]
print('\ntimeseries2velocity.py', ' '.join(iargs))
mintpy.cli.timeseries2velocity.main(iargs)
def run_geocode(self, step_name):
"""geocode data files in radar coordinates into ./geo folder."""
if self.template['mintpy.geocode']:
ts_file = self.get_timeseries_filename(self.template, self.workDir)[step_name]['input']
atr = readfile.read_attribute(ts_file)
if 'Y_FIRST' not in atr.keys():
# 1. geocode
out_dir = os.path.join(self.workDir, 'geo')
os.makedirs(out_dir, exist_ok=True)
geom_file, lookup_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1:3]
in_files = [geom_file, 'temporalCoherence.h5', 'avgSpatialCoh.h5', ts_file, 'velocity.h5']
iargs = in_files + ['-l', lookup_file, '-t', self.templateFile, '--outdir', out_dir, '--update']
print('\ngeocode.py', ' '.join(iargs))
import mintpy.cli.geocode
mintpy.cli.geocode.main(iargs)
# 2. generate reliable pixel mask in geo coordinate
geom_file = os.path.join(out_dir, f'geo_{os.path.basename(geom_file)}')
tcoh_file = os.path.join(out_dir, 'geo_temporalCoherence.h5')
mask_file = os.path.join(out_dir, 'geo_maskTempCoh.h5')
tcoh_min = self.template['mintpy.networkInversion.minTempCoh']
iargs = [tcoh_file, '-m', tcoh_min, '-o', mask_file]
# exclude pixels in shadow if shadowMask dataset is available
if (self.template['mintpy.networkInversion.shadowMask'] is True
and 'shadowMask' in readfile.get_dataset_list(geom_file)):
iargs += ['--base', geom_file, '--base-dataset', 'shadowMask', '--base-value', '1']
print('\ngenerate_mask.py', ' '.join(iargs))
if ut.run_or_skip(out_file=mask_file, in_file=tcoh_file) == 'run':
import mintpy.cli.generate_mask
mintpy.cli.generate_mask.main(iargs)
else:
print('dataset is geocoded, skip geocoding and continue.')
else:
print('geocoding is OFF')
def run_save2google_earth(self, step_name):
"""Save velocity file in geo coordinates into Google Earth raster image."""
if self.template['mintpy.save.kmz'] is True:
print('creating Google Earth KMZ file for geocoded velocity file: ...')
# input
vel_file = os.path.join(self.workDir, 'velocity.h5')
atr = readfile.read_attribute(vel_file)
if 'Y_FIRST' not in atr.keys():
vel_file = os.path.join(self.workDir, 'geo/geo_velocity.h5')
# output
kmz_file = f'{os.path.splitext(vel_file)[0]}.kmz'
iargs = [vel_file, '-o', kmz_file]
print('\nsave_kmz.py', ' '.join(iargs))
# update mode
fbase = os.path.basename(kmz_file)
kmz_files = [i for i in [fbase, f'./geo/{fbase}', f'./pic/{fbase}']
if os.path.isfile(i)]
kmz_file = kmz_files[0] if len(kmz_files) > 0 else None
if ut.run_or_skip(out_file=kmz_file, in_file=vel_file, readable=False) == 'run':
import mintpy.cli.save_kmz
mintpy.cli.save_kmz.main(iargs)
else:
print('save velocity to Google Earth format is OFF.')
def run_save2hdfeos5(self, step_name):
"""Save displacement time-series and its aux data in geo coordinate into HDF-EOS5 format"""
if self.template['mintpy.save.hdfEos5'] is True:
# input
ts_file = self.get_timeseries_filename(self.template, self.workDir)[step_name]['input']
# Add attributes from custom template to timeseries file
if self.customTemplate is not None:
ut.add_attribute(ts_file, self.customTemplate)
tcoh_file = os.path.join(self.workDir, 'temporalCoherence.h5')
scoh_file = os.path.join(self.workDir, 'avgSpatialCoh.h5')
mask_file = os.path.join(self.workDir, 'maskTempCoh.h5')
geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1]
if 'geo' in ts_file:
tcoh_file = os.path.join(self.workDir, 'geo/geo_temporalCoherence.h5')
scoh_file = os.path.join(self.workDir, 'geo/geo_avgSpatialCoh.h5')
mask_file = os.path.join(self.workDir, 'geo/geo_maskTempCoh.h5')
geom_file = os.path.join(self.workDir, f'geo/geo_{os.path.basename(geom_file)}')
# cmd
print('--------------------------------------------')
iargs = [ts_file,
'--tc', tcoh_file,
'--asc', scoh_file,
'-m', mask_file,
'-g', geom_file,
'-t', self.templateFile]
print('\nsave_hdfeos5.py', ' '.join(iargs))
# output (check existing file)
atr = readfile.read_attribute(ts_file)
SAT = sensor.get_unavco_mission_name(atr)
hdfeos5_files = ut.get_file_list(f'{SAT}_*.he5')
hdfeos5_file = hdfeos5_files[0] if len(hdfeos5_files) > 0 else None
if ut.run_or_skip(out_file=hdfeos5_file,
in_file=[ts_file, tcoh_file,
scoh_file, mask_file,
geom_file]) == 'run':
import mintpy.cli.save_hdfeos5
mintpy.cli.save_hdfeos5.main(iargs)
else:
print('save time-series to HDF-EOS5 format is OFF.')
def run(self, steps):
"""run the chosen steps."""
for sname in steps:
print(f'\n\n******************** step - {sname} ********************')
if sname == 'load_data':
self.run_load_data(sname)
elif sname == 'modify_network':
self.run_network_modification(sname)
elif sname == 'reference_point':
self.run_reference_point(sname)
elif sname == 'quick_overview':
self.run_quick_overview(sname)
elif sname == 'correct_unwrap_error':
self.run_unwrap_error_correction(sname)
elif sname == 'invert_network':
self.run_network_inversion(sname)
elif sname == 'correct_LOD':
self.run_local_oscillator_drift_correction(sname)
elif sname == 'correct_SET':
self.run_solid_earth_tides_correction(sname)
elif sname == 'correct_troposphere':
self.run_tropospheric_delay_correction(sname)
elif sname == 'deramp':
self.run_phase_deramping(sname)
elif sname == 'correct_topography':
self.run_topographic_residual_correction(sname)
elif sname == 'residual_RMS':
self.run_residual_phase_rms(sname)
elif sname == 'reference_date':
self.run_reference_date(sname)
elif sname == 'velocity':
self.run_timeseries2velocity(sname)
elif sname == 'geocode':
self.run_geocode(sname)
elif sname == 'google_earth':
self.run_save2google_earth(sname)
elif sname == 'hdfeos5':
self.run_save2hdfeos5(sname)
def plot_result(self, print_aux=True):
"""Plot data files and save to figures in pic folder"""
print('\n******************** plot & save to pic ********************')
tropo_model = self.template['mintpy.troposphericDelay.weatherModel'].upper()
max_plot_memory = abs(float(self.template['mintpy.plot.maxMemory']))
max_memory = abs(float(self.template['mintpy.compute.maxMemory']))
fig_dpi = int(self.template['mintpy.plot.dpi'])
stack_file, geom_file, lookup_file, ion_file = ut.check_loaded_dataset(
self.workDir,
print_msg=False)[:4]
mask_file = os.path.join(self.workDir, 'maskTempCoh.h5')
geo_dir = os.path.join(self.workDir, 'geo')
pic_dir = os.path.join(self.workDir, 'pic')
# use relative path for shorter and cleaner printout view command
stack_file = os.path.relpath(stack_file) if stack_file else stack_file
ion_file = os.path.relpath(ion_file) if ion_file else ion_file
geom_file = os.path.relpath(geom_file) if geom_file else geom_file
lookup_file = os.path.relpath(lookup_file) if lookup_file else lookup_file
mask_file = os.path.relpath(mask_file) if mask_file else mask_file
geo_dir = os.path.relpath(geo_dir) if geo_dir else geo_dir
# view options
# for each element list:
# the 1st item is the data file
# the 2nd item is the dataset if applicable
opt4ts = ['--noaxis', '-u', 'cm', '--wrap', '--wrap-range', '-5', '5']
iargs_list0 = [
# key files
['velocity.h5', '--dem', geom_file, '--mask', mask_file],
['temporalCoherence.h5', '-c', 'gray', '-v', '0', '1'],
['maskTempCoh.h5', '-c', 'gray', '-v', '0', '1'],
# geometry
[geom_file],
[lookup_file],
# ifgramStack
[stack_file, 'unwrapPhase-', '--zero-mask', '--wrap', '-c', 'cmy'],
[stack_file, 'unwrapPhase-', '--zero-mask'],
[stack_file, 'coherence-', '--mask', 'no', '-v', '0', '1'],
[stack_file, 'connectComponent-', '--mask', 'no'],
# ifgramStack - unwrapping error correction
[stack_file, 'unwrapPhase_bridging-', '--zero-mask'],
[stack_file, 'unwrapPhase_phaseClosure-', '--zero-mask'],
[stack_file, 'unwrapPhase_bridging_phaseClosure-', '--zero-mask'],
# ifgramStack - auxliary files
['avgPhaseVelocity.h5'] ,
['avgSpatialCoh.h5', '-c', 'gray', '-v', '0', '1'],
['maskConnComp.h5', '-c', 'gray', '-v', '0', '1'],
# time-series
['timeseries.h5'] + opt4ts,
['timeseries_*.h5'] + opt4ts,
# files from geocoding
[os.path.join(geo_dir, 'geo_maskTempCoh.h5'), '-c', 'gray'],
[os.path.join(geo_dir, 'geo_temporalCoherence.h5'), '-c', 'gray'],
[os.path.join(geo_dir, 'geo_avgSpatialCoh.h5'), '-c', 'gray'],
[os.path.join(geo_dir, 'geo_velocity.h5'), 'velocity'],
[os.path.join(geo_dir, 'geo_timeseries*.h5')] + opt4ts,
# all the other files
[f'velocity{tropo_model}.h5', '--mask', 'no'],
['numInvIfgram.h5', '--mask', 'no'],
]
if ion_file:
iargs_list0 += [
[ion_file, 'unwrapPhase-', '--zero-mask'],
[ion_file, 'coherence-', '--mask', 'no', '-v', '0', '1'],
]
# translate element list whose file path has *
iargs_list = []
for iargs in iargs_list0:
fname, args = iargs[0], iargs[1:]
if not fname:
continue
if '*' in fname:
fnames = sorted(glob.glob(fname))
if len(fnames) > 0:
for fname in fnames:
iargs_list.append([fname] + args)
elif iargs not in iargs_list:
iargs_list.append(iargs)
# remove element list - file does not exists
iargs_list = [iargs for iargs in iargs_list if os.path.isfile(iargs[0])]
# remote element list - file is ifgramStack and the dataset of interest does not exist
stack_dset_list = readfile.get_dataset_list(stack_file)
stack_dset_list += [f'{x}-' for x in stack_dset_list]
iargs_list = [iargs for iargs in iargs_list
if (iargs[0] != stack_file
or (iargs[0] == stack_file
and iargs[1] in stack_dset_list))]
# add the following common options to all element lists
opt_common = ['--dpi', str(fig_dpi), '--noverbose', '--nodisplay',
'--update', '--memory', str(max_plot_memory)]
iargs_list = [opt_common + iargs for iargs in iargs_list]
# run view
start_time = time.time()
run_parallel = False
cluster_type = self.template['mintpy.compute.cluster']
if cluster_type:
num_workers = self.template['mintpy.compute.numWorker']
num_workers = cluster.DaskCluster.format_num_worker(cluster_type, num_workers)
# limit number of parallel processes based on available CPU
num_cores, run_parallel, Parallel, delayed = ut.check_parallel(
len(iargs_list),
print_msg=False,
maxParallelNum=num_workers)
# limit number of parallel processes based on memory limit
# default view.py call could use up to 1.5 GB reserved memory (~700 MB actual memory)
# scale based on utils.plot.auto_multilook_num()
plot_memory = 1.5 if 2.0 < max_plot_memory <= 4.0 else 1.5 * (max_plot_memory / 4.0)
num_cores = min(num_cores, max(int(max_memory / plot_memory), 1))
import mintpy.cli.view
if run_parallel and num_cores > 1:
print(f"parallel processing using {num_cores} cores ...")
Parallel(n_jobs=num_cores)(delayed(mintpy.cli.view.main)(iargs) for iargs in iargs_list)
else:
for iargs in iargs_list:
mintpy.cli.view.main(iargs)
# copy text files to pic
print('copy *.txt files into ./pic directory.')
tfiles = glob.glob('*.txt')
for tfile in tfiles:
shutil.copy2(tfile, pic_dir)
# move picture files to pic
print('move *.png/pdf/kmz files to ./pic directory.')
pfiles = glob.glob('*.png')
pfiles += glob.glob('*.pdf')
pfiles += glob.glob('*.kmz')
pfiles += glob.glob(os.path.join(geo_dir, '*.kmz'))
for pfile in pfiles:
shutil.move(pfile, os.path.join(pic_dir, os.path.basename(pfile)))
# time info
m, s = divmod(time.time()-start_time, 60)
print(f'time used: {m:02.0f} mins {s:02.1f} secs.')
# message for more visualization scripts
msg = """Explore more info & visualization options with the following scripts:
info.py #check HDF5 file structure and metadata
view.py #2D map view
tsview.py #1D point time-series (interactive)
transect.py #1D profile (interactive)
plot_coherence_matrix.py #plot coherence matrix for one pixel (interactive)
plot_network.py #plot network configuration of the dataset
plot_transection.py #plot 1D profile along a line of a 2D matrix (interactive)
save_kmz.py #generate Google Earth KMZ file in raster image
save_kmz_timeseries.py #generate Google Earth KMZ file in points for time-series (interactive)
"""
if print_aux:
print(msg)
def close(self, normal_end=True):
# go back to original directory
print('Go back to directory:', self.cwd)
os.chdir(self.cwd)
# message
if normal_end:
msg = '\n################################################'
msg += '\n Normal end of smallbaselineApp processing!'
msg += '\n################################################'
print(msg)
def run_smallbaselineApp(inps):
"""Run the small baseline time series analsysis workflow."""
start_time = time.time()
# open and run
app = TimeSeriesAnalysis(inps.customTemplateFile, inps.workDir)
app.open()
app.run(steps=inps.runSteps)
# plot if:
# a) --plot in command line OR
# b) template['mintpy.plot'] = yes AND runSteps > 1
if inps.plot or (app.template['mintpy.plot'] and len(inps.runSteps) > 1):
app.plot_result()
# close
app.close()
# used time
m, s = divmod(time.time()-start_time, 60)
print(f'Time used: {m:02.0f} mins {s:02.1f} secs\n')
return
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...