Prepare ISCE topsApp interferograms for processing with MintPy

• • Prepare the structured folders
  • Link files
  • generate baseline
  • prep_isce.py
  • write options file for MintPy

We often use topsApp.pyto process sentinel-1 data. However, MintPy can only read the results of topsStackdirectly at present. We need some post-processing to make MintPy can produce the result of topsApp.py.

Prepare the structured folders

The folder structure required by Minpy to handle tops is as follows ISCE / topsStack:

$DATA_DIR/GalapagosSenDT128├── baselines│   ├── 20141213_20141225│   │   └── 20141213_20141225.txt│   └── ...├── reference│   ├── data.rsc    #generated by prep_isce.py│   ├── IW1│   ├── IW1.xml│   ├── IW2│   └── IW2.xml├── merged│   ├── geom_reference│   │   ├── hgt.rdr│   │   ├── hgt.rdr.full.vrt│   │   ├── hgt.rdr.full.xml│   │   ├── hgt.rdr.vrt│   │   ├── hgt.rdr.xml│   │   ├── lat.rdr│   │   ├── lat.rdr.full.vrt│   │   ├── lat.rdr.full.xml│   │   ├── lat.rdr.vrt│   │   ├── lat.rdr.xml│   │   ├── lon.rdr│   │   ├── lon.rdr.full.vrt│   │   ├── lon.rdr.full.xml│   │   ├── lon.rdr.vrt│   │   ├── lon.rdr.xml│   │   ├── los.rdr│   │   ├── los.rdr.full.vrt│   │   ├── los.rdr.full.xml│   │   ├── los.rdr.vrt│   │   ├── los.rdr.xml│   │   ├── shadowMask.rdr│   │   ├── shadowMask.rdr.full.vrt│   │   ├── shadowMask.rdr.full.xml│   │   ├── shadowMask.rdr.vrt│   │   └── shadowMask.rdr.xml│   └── interferograms│       ├── 20141213_20141225│       │   ├── filt_fine.cor│       │   ├── filt_fine.cor.vrt│       │   ├── filt_fine.cor.xml│       │   ├── filt_fine.unw│       │   ├── filt_fine.unw.conncomp│       │   ├── filt_fine.unw.conncomp.vrt│       │   ├── filt_fine.unw.conncomp.xml│       │   ├── filt_fine.unw.vrt│       │   ├── filt_fine.unw.xml│       │   ├── ...│       ├── 20141213_20150307│       └── ...├── secondarys│   ├── 20141225│   │   ├── IW1│   │   ├── IW1.xml│   │   ├── IW2│   │   └── IW2.xml│   ├── 20150307│   └── ...└── mintpy    └── GalapagosSenDT128.txt

Suppose you process all interferograms in one folder, which has the folder structure as follows:

.|-- 20210917_20210929|-- 20210917_20211011|-- 20210917_20211023|-- 20210917_20211210|-- 20210929_20211011|-- 20210929_20211023|-- 20210929_20211104|-- ...

Firstly, import packages and specify the home directory.

import rasterioimport osfrom pathlib import Pathimport warningsfrom tqdm import tqdmimport numpy as np# the folder that contains the interferograms processed by topsApp.pyifg_home = Path('/data/sentinel1/result/interferogams')  # the folder that used to run MintPystack_home = Path('/data/sentinel1/result/stacks')

Now, let’s create the folder structure that MintPy required.

############  make dirs #############def ensure_folder(folder):    if not folder.is_dir():        folder.mkdir(parents=True)baseline_dir = stack_home / 'baselines'reference_dir = stack_home / 'reference'secondarys_dir = stack_home / 'secondarys'mintpy_dir = stack_home / 'mintpy'merged_dir = stack_home / 'merged'geom_reference_dir = merged_dir / 'geom_reference'interferograms_dir = merged_dir / 'interferograms'stack_fld_list = [baseline_dir, reference_dir,                  secondarys_dir, mintpy_dir,                  merged_dir, geom_reference_dir,                  interferograms_dir]for fld in stack_fld_list:    ensure_folder(fld)

Link files

We have created the folders that MintPy required. But there are no files in those folders. We need to link files from ISCE processed into those folders.

def multilook(array, nlook_r, nlook_c, nodata, n_valid_thre=0.5):    """    Nodata in input array must be filled with nan beforehand.    if the number of valid data is less than n_valid_thre*nlook_r*nlook_c, return nan.    """    length, width = array.shape    length_ml = int(np.floor(length / nlook_r))    width_ml = int(np.floor(width / nlook_c))    array_reshape = array[: length_ml * nlook_r, : width_ml * nlook_c].reshape(        length_ml, nlook_r, width_ml, nlook_c    )    with warnings.catch_warnings():  # To silence RuntimeWarning: Mean of empty slice        warnings.simplefilter("ignore", RuntimeWarning)        array_ml = np.nanmean(array_reshape, axis=(1, 3))    n_valid = np.sum(~np.isnan(array_reshape), axis=(1, 3))    bool_invalid = n_valid < n_valid_thre * nlook_r * nlook_c    array_ml[bool_invalid] = nodata    return array_mldef vrt2rdr(vrt_file, rdr_file, ref_file, nlook_r, nlook_c):    with rasterio.open(vrt_file) as src:        arr = src.read(1)        arr_ml = multilook(arr, nlook_r, nlook_c, src.nodata)        meta = src.meta        with rasterio.open(ref_file) as ref_ds:            meta.update({                'crs': ref_ds.crs,                'transform': ref_ds.transform,                'width': ref_ds.width,                'height': ref_ds.height,                'driver': 'ISCE'})        with rasterio.open(rdr_file, 'w', **meta) as dst:            dst.write(arr_ml,1)ifgs_files = sorted(ifg_home.iterdir())for file in tqdm(ifgs_files, desc='link files'):    name = file.name    ref_name, sec_name = name.split('_')    merged = file / 'merged'    i = 0    if merged.exists():        os.chdir(file)        if i == 0:            ####################################################################            # For the folder 'reference' and 'merged/geom_reference', you only             # need to get it from the first interferogram                        #####################################################################            i += 1                        ######### link to 'reference' folder ######            ref_dir = file / ref_name            ref_cmd = f'ln -sf {reference_dir}/'            os.system(ref_cmd)            ####################################################################            #     merge and convert geo_info to 'geom_reference' folder            #####################################################################            # For the files that contains the geo info of ifgs, only 'los' was             # processed well (multi-looked and with '.rdr' format). We need             # to convert the full resolution geo info file to desired format            # (multi-look and convert to '.rdr' format). This function was             # named vrt2rdr. in the script, azimuth/range looks are: 4*20.            # the 'los.rdr' file was used to get the geo properties.                     #####################################################################            vrt2rdr(merged.joinpath('z.rdr.full.vrt'),                    geom_reference_dir.joinpath('hgt.rdr'),                    merged.joinpath('los.rdr'),                    4,20)            vrt2rdr(merged.joinpath('lat.rdr.full.vrt'),                    geom_reference_dir.joinpath('lat.rdr'),                    merged.joinpath('los.rdr'),                    4,20)            vrt2rdr(merged.joinpath('lon.rdr.full.vrt'),                    geom_reference_dir.joinpath('lon.rdr'),                    merged.joinpath('los.rdr'),                    4,20)            vrt2rdr(merged.joinpath('los.rdr.full.vrt'),                    geom_reference_dir.joinpath('los.rdr'),                    merged.joinpath('los.rdr'),                    4,20)        # ########## link ifgs  ################        ifg_out = interferograms_dir / name        ensure_folder(ifg_out)        # ifgs        ln_cmd1 = f'ln -sf {ifg_out}/'        os.system(ln_cmd1)        # cor        ln_cmd2 = f'ln -sf {ifg_out}/'        os.system(ln_cmd2)        # ########## link secondarys  ################        sec_dir = file / sec_name        sec_cmd = f'ln -sf {secondarys_dir}/'        os.system(sec_cmd)

generate baseline

We can use isce2/contrib/stack/topsStack/computeBaseline.pyto generate baselines. This script was not installed by default, for example, installed by conda. You need to find its location from the source code of ISCE.

########## generate_baseline  #############def get_dates_from_ifgs(ifgs):    dates_all = []    for ifg in ifgs:        dates_all.extend(ifg.split('_'))    dates = sorted(set(dates_all))    return datesifgs = [i.name for i in ifgs_files]dates = get_dates_from_ifgs(ifgs)ref_name = dates[0]ref_dir = ifgs_files[0] / ref_namefor sec_name in tqdm(dates[1:], desc='generate baselines'):    sec_dir = secondarys_dir / sec_name    bl_name = f'{sec_name}'    bl_dir = baseline_dir / bl_name    ensure_folder(bl_dir)    bl_file = bl_dir / f'{bl_name}.txt'    bl_cmd = f'/data/isce2/contrib/stack/topsStack/computeBaseline.py -m {bl_file}'    os.system(bl_cmd)

prep_isce.py

For the folder reference, we need to use prep_isce.pyto generate the data.rscfile.

os.chdir(reference_dir)prep_cmd = ('/root/miniconda3/bin/prep_isce.py'            f' -f "{interferograms_dir}/*/filt_*.unw"'            f' -m {reference_dir}/IW2.xml'            f' -b {baseline_dir}/'            f' -g {geom_reference_dir}/')print(prep_cmd)os.system(prep_cmd)

If everything is ok, no errors will be reported.

write options file for MintPy

cfg_file = mintpy_dir / 'mintpy_data_options.cfg'cfg_info = f'''mintpy.load.processor        = isce##---------for ISCE only:mintpy.load.metaFile         = {reference_dir}/IW*.xmlmintpy.load.baselineDir      = {baseline_dir}##---------interferogram datasets:mintpy.load.unwFile          = {interferograms_dir}/*/filt_*.unwmintpy.load.corFile          = {interferograms_dir}/*/phsig.cormintpy.load.connCompFile     = {interferograms_dir}/*/filt_*.unw.conncomp##---------geometry datasets:mintpy.load.demFile          = {geom_reference_dir}/hgt.rdrmintpy.load.lookupYFile      = {geom_reference_dir}/lat.rdrmintpy.load.lookupXFile      = {geom_reference_dir}/lon.rdrmintpy.load.incAngleFile     = {geom_reference_dir}/los.rdrmintpy.load.azAngleFile      = {geom_reference_dir}/los.rdr# mintpy.load.shadowMaskFile   = {geom_reference_dir}/shadowMask.rdr'''with open(cfg_file, 'w') as f:    f.write(cfg_info)