[Feature] Silent all-zero output from invert_network with custom LiCSAR GeoTIFF reader

[Liuqi2002519]

Description of the desired feature
I am attempting to process LiCSAR interferometric data with MintPy as an alternative to the standard LiCSBAS workflow. Since MintPy does not natively support LiCSAR inputs, I wrote a custom prep_licsar script modeled after prep_gmtsar, placed it in bin/ and mintpy/cli/, and verified it executes globally and correctly generates .rsc metadata files.
Pre-processing steps applied to the LiCSAR products:

Unwrapped interferograms: read directly from GeoTIFF (.tif).
Coherence (CC): rescaled from the native 0–255 integer range to 0–1 float before ingestion.
Baseline: retrieved from the ASF baseline tool and reformatted to be MintPy-compatible.

Is your feature request related to a problem? Please describe
The pipeline runs without errors through load_data → modify_network → reference_point → quick_overview and produces a visually plausible avgPhaseVelocity.png that closely matches the velocity map independently generated with LiCSBAS.
However, invert_network (default method: LSQR) completes without any error or warning, yet both timeseries.h5 and temporalCoherence.h5 contain only zeros, and all corresponding .png outputs are blank.
Describe the solution you'd like
I would like to understand why invert_network silently produces an all-zero time series. Two hypotheses in order of suspicion:
Hypothesis 1 (more likely) — all pixels are silently masked before inversion.
MintPy uses 0 as its internal no-data/masked value. It is possible that the coherence-based mask (or the combined mask from maskPS.h5 / maskTempCoh.h5) is discarding every pixel before the LSQR solver runs, causing it to write zeros without raising an exception. Potential causes:

The coherence threshold in smallbaselineApp.cfg is higher than the effective values in my rescaled CC data.
A unit or dtype mismatch in the .tif coherence files causes MintPy to misread the 0–1 float values (e.g., reading them as integers, making everything effectively 0).

Hypothesis 2 (less likely) — a subtle format issue specific to the inversion stage.
avgPhaseVelocity.png is computed directly from the raw interferogram stack before masking or inversion, so its looking correct does not rule out a format problem that only surfaces when the LSQR network system is assembled (e.g., wrong phase unit, byte order issue in the GeoTIFF reader, or an inconsistency between the interferogram and coherence spatial extents).
Ideally this could also serve as a starting point for an official prep_licsar reader in MintPy, given that LiCSAR is a major publicly available Sentinel-1 archive.
Describe alternatives you have considered

Using LiCSBAS directly — works and produces consistent velocity maps, confirming the raw data is valid. I specifically want MintPy's inversion and time series analysis capabilities.
Adapting prep_aria instead of prep_gmtsar as the template — not yet attempted.

Additional context
Key diagnostic observations:
StageResultload_dataCompletes normally; HDF5 stack looks correctavgPhaseVelocity.pngVisually consistent with LiCSBAS outputinvert_network (LSQR, default)Exits cleanly, no errors or warningstimeseries.h5 / temporalCoherence.h5All zeros
Suggested debugging steps I have not yet been able to isolate:

Is there a way to print the number of valid (unmasked) pixels entering the LSQR solver, to confirm whether masking is the cause?
Could a float32 vs float64 mismatch in the rescaled coherence TIFFs silently produce an all-zero mask?

Happy to share the prep_licsar script, example input files, smallbaselineApp.cfg, and the full mintpy.log if that helps reproduce the issue.
Are you willing to help implement and maintain this feature?

Yes

prep_licsar.py
#!/usr/bin/env python
############################################################

Derived from MintPy prep_gmtsar.py and prep_SARscape.py

Baseline: GMTSAR baseline_table.dat format

Metadata: SARscape .sml / .prm format

Data: GeoTIFF (.tif) files

############################################################

import argparse
import datetime
import glob
import os
import sys

import numpy as np

try:
from osgeo import gdal
gdal.UseExceptions()
except ImportError:
raise ImportError('Cannot import gdal!')

from mintpy.utils import ptime, readfile, utils as ut, writefile

##############################################################

-------------------- Baseline (GMTSAR format) -----------

##############################################################

def read_baseline_table(fname):
print(f'read baseline time-series from file: {fname}')
fc = np.loadtxt(fname, dtype=str, usecols=(1, 4))

if fc.ndim == 1:
    fc = fc[np.newaxis, :]

# 用datetime转换，避免ptime.yyyyddd2yyyymmdd的off-by-one问题
date_list = []
for x in fc[:, 0]:
    yyyyddd = x[:7]
    year = int(yyyyddd[:4])
    doy  = int(yyyyddd[4:])
    dt   = datetime.datetime(year, 1, 1) + datetime.timedelta(days=doy - 1)
    date_list.append(dt.strftime('%Y%m%d'))

bperp_list = fc[:, 1].astype(np.float32)

bdict = {}
for date_str, bperp in zip(date_list, bperp_list):
    bdict[date_str] = bperp

print(f'  Found {len(bdict)} acquisitions: '
      f'{sorted(bdict.keys())[0]} ... {sorted(bdict.keys())[-1]}')
return bdict

##############################################################

-------------------- Metadata (SARscape PRM format) -----

##############################################################

SARscape SML key -> MintPy standard key

_SARSCAPE_KEY_MAP = {
# Orbit / geometry
'OrbitDirection' : 'ORBIT_DIRECTION',
'PassDirection' : 'ORBIT_DIRECTION', # fallback alias
# Sensor
'SensorName' : 'PLATFORM',
'NadirHeight' : 'HEIGHT',
'OrbitHeight' : 'HEIGHT', # fallback alias
'CenterFrequency' : None, # handled separately -> WAVELENGTH
'Wavelength' : 'WAVELENGTH',
'RangePixelSpacing' : 'RANGE_PIXEL_SIZE',
'AzimuthPixelSpacing' : 'AZIMUTH_PIXEL_SIZE',
'NearSlantRange' : 'STARTING_RANGE',
'EarthRadius' : 'EARTH_RADIUS',
# Looks
'RangeLooks' : 'RLOOKS',
'AzimuthLooks' : 'ALOOKS',
# Dates
'MasterAcquisitionDate' : 'DATE12', # will be post-processed
'SlaveAcquisitionDate' : None,
}

Speed of light (m/s) for freq -> wavelength conversion

_C = 299_792_458.0

def read_sarscape_prm(prm_file):
"""Parse a SARscape SML/PRM metadata file (key = value format)."""
meta = {}
print(f' Reading SARscape PRM: {prm_file}')
with open(prm_file, 'r', encoding='utf-8', errors='replace') as fh:
for line in fh:
line = line.strip()
if not line or line.startswith('#') or line.startswith(';'):
continue
if '=' in line:
k, _, v = line.partition('=')
meta[k.strip()] = v.strip()
print(f' Read {len(meta)} raw keys from PRM file')
return meta

def _sarscape_to_mintpy(raw):
"""Convert raw SARscape PRM dict to MintPy-standard metadata dict."""
meta = {}

for src_key, dst_key in _SARSCAPE_KEY_MAP.items():
    if src_key not in raw:
        continue
    val = raw[src_key]

    if src_key == 'CenterFrequency' and val:
        try:
            freq = float(val)
            meta['WAVELENGTH'] = str(_C / freq)
        except ValueError:
            pass
        continue

    if dst_key and val:
        if dst_key == 'ORBIT_DIRECTION':
            v_low = val.lower()
            if v_low.startswith('asc'):
                meta[dst_key] = 'ascending'
            elif v_low.startswith('desc'):
                meta[dst_key] = 'descending'
            else:
                meta[dst_key] = val
        else:
            meta[dst_key] = val

if 'HEIGHT' not in meta or not meta.get('HEIGHT'):
    for alt_key in ('NadirHeight', 'OrbitHeight'):
        if raw.get(alt_key):
            meta['HEIGHT'] = raw[alt_key]
            print(f'  HEIGHT set from {alt_key}: {meta["HEIGHT"]}')
            break

if 'HEIGHT' not in meta or not meta.get('HEIGHT'):
    _platform_height = {
        'SENTINEL': '693000.0',
        'S1'      : '693000.0',
        'ALOS2'   : '628370.0',
        'ALOS-2'  : '628370.0',
        'TSX'     : '514000.0',
        'TDX'     : '514000.0',
        'COSMO'   : '619600.0',
        'CSK'     : '619600.0',
        'ICEYE'   : '561000.0',
        'NISAR'   : '747000.0',
    }
    platform = meta.get('PLATFORM', raw.get('SensorName', '')).upper()
    for kw, ht in _platform_height.items():
        if kw in platform:
            meta['HEIGHT'] = ht
            print(f'  HEIGHT inferred from PLATFORM ({platform}): {ht} m')
            break

return meta

def extract_sarscape_metadata(prm_file):
"""High-level wrapper: read SARscape PRM and return MintPy metadata dict."""
raw = read_sarscape_prm(prm_file)
meta = _sarscape_to_mintpy(raw)
meta['PROCESSOR'] = 'licsar'
return meta

##############################################################

-------------------- TIF helpers -----------------------

##############################################################

def fix_nodata_tif(tif_file):
"""原地将tif中的nodata值替换为NaN（float32）。

读取GDAL记录的NoDataValue，将所有等于该值的像素替换为NaN，
并更新波段的NoDataValue为NaN，避免MintPy读入时把无效值当真值。

Parameters: tif_file - str, path to the GeoTIFF (modified in-place)
"""
ds = gdal.Open(tif_file, gdal.GA_Update)
if ds is None:
    print(f'  WARNING: Cannot open for update: {tif_file}')
    return

for b in range(1, ds.RasterCount + 1):
    band   = ds.GetRasterBand(b)
    nodata = band.GetNoDataValue()
    if nodata is None:
        continue

    data = band.ReadAsArray().astype(np.float32)
    n    = np.sum(data == nodata)
    if n > 0:
        data[data == nodata] = np.nan
        band.WriteArray(data)
        band.SetNoDataValue(float('nan'))
        print(f'  {os.path.basename(tif_file)} band{b}: '
              f'{n} nodata({nodata}) pixels -> NaN')

ds.FlushCache()
ds = None

def read_tif_spatial_meta(tif_file):
"""Extract spatial metadata from a GeoTIFF via GDAL."""
ds = gdal.Open(tif_file, gdal.GA_ReadOnly)
if ds is None:
raise FileNotFoundError(f'GDAL cannot open: {tif_file}')

gt  = ds.GetGeoTransform()
nx  = ds.RasterXSize
ny  = ds.RasterYSize
ds  = None

meta = {
    'WIDTH'  : str(nx),
    'LENGTH' : str(ny),
    'X_FIRST': str(gt[0]),
    'Y_FIRST': str(gt[3]),
    'X_STEP' : str(gt[1]),
    'Y_STEP' : str(gt[5]),
    'X_UNIT' : 'degrees',
    'Y_UNIT' : 'degrees',
}
return meta

def get_image_corners(tif_file, orbit_direction):
"""Compute LAT/LON_REF1/2/3/4 corner metadata from a GeoTIFF."""
ds = gdal.Open(tif_file, gdal.GA_ReadOnly)
gt = ds.GetGeoTransform()
W = gt[0]
N = gt[3]
dx = abs(gt[1])
dy = abs(gt[5])
E = W + dx * ds.RasterXSize
S = N - dy * ds.RasterYSize
ds = None

meta = {}
if orbit_direction.lower().startswith('asc'):
    meta.update(LAT_REF1=str(S), LAT_REF2=str(S),
                LAT_REF3=str(N), LAT_REF4=str(N),
                LON_REF1=str(W), LON_REF2=str(E),
                LON_REF3=str(W), LON_REF4=str(E))
else:
    meta.update(LAT_REF1=str(N), LAT_REF2=str(N),
                LAT_REF3=str(S), LAT_REF4=str(S),
                LON_REF1=str(E), LON_REF2=str(W),
                LON_REF3=str(E), LON_REF4=str(W))
return meta

def compute_incidence_angle(meta):
"""Estimate scene-centre incidence angle from orbit geometry."""
try:
Re = float(meta['EARTH_RADIUS'])
H = float(meta['HEIGHT'])
Rn = float(meta['STARTING_RANGE'])
dr = float(meta['RANGE_PIXEL_SIZE'])
w = float(meta['WIDTH'])
except (KeyError, ValueError) as e:
print(f'WARNING: Cannot compute incidence angle ({e}), skipping.')
return meta

Rg  = Rn + dr * w / 2.0
inc = (np.pi - np.arccos((Re**2 + Rg**2 - (Re + H)**2) / (2 * Re * Rg))) \
      * 180.0 / np.pi
meta['SLANT_RANGE_DISTANCE'] = str(Rg)
meta['INCIDENCE_ANGLE']      = str(inc)
return meta

##############################################################

-------------------- Geometry preparation ---------------

##############################################################

GEOM_KEYS = [
'demFile', 'lookupYFile', 'lookupXFile',
'incAngleFile', 'azAngleFile', 'shadowMaskFile', 'waterMaskFile',
]

def _resolve_file(pattern):
"""Resolve a file path or glob pattern to a single existing file path."""
if not pattern or str(pattern).lower() in ('auto', 'none', '', 'false', 'true'):
return None
pattern = str(pattern).strip("'"")
if os.path.isfile(pattern):
return pattern
matches = sorted(glob.glob(pattern))
return matches[0] if matches else None

def prepare_geometry(template, common_meta, update_mode=True):
"""Write .rsc sidecar files for all geometry TIF files in the template."""
prefix = 'mintpy.load.'
geom_files = []
for key in GEOM_KEYS:
fpath = _resolve_file(template.get(prefix + key, ''))
if fpath:
geom_files.append(fpath)
print(f' Found geometry file ({key}): {fpath}')
else:
print(f' WARNING: No file for {key}')

if not geom_files:
    print('WARNING: No geometry files found – skipping geometry preparation.')
    return

for geom_file in geom_files:
    # 原地将nodata值替换为NaN，防止MintPy把-9999当真值读入
    fix_nodata_tif(geom_file)

    geom_meta = dict(common_meta)
    geom_meta.update(read_tif_spatial_meta(geom_file))

    rsc_file = geom_file + '.rsc'
    writefile.write_roipac_rsc(
        geom_meta, rsc_file,
        update_mode=update_mode, print_msg=True,
    )

##############################################################

-------------------- Interferogram stack preparation ----

##############################################################

def dates_from_dir(ifg_dir):
"""Extract (date1, date2) from interferogram directory name."""
name = os.path.basename(ifg_dir)
tokens = name.replace('-', '').split('_')
dates = [t for t in tokens if len(t) == 8 and t.isdigit()]
if len(dates) >= 2:
return dates[0], dates[1]
return None

def prepare_stack(unw_files, common_meta, bperp_dict, update_mode=True):
"""Write .rsc sidecar files for all unwrapped TIF interferograms."""
n_total = len(unw_files)
if n_total == 0:
raise FileNotFoundError('No unwrapped interferogram files found!')

bperp_dates = sorted(bperp_dict.keys())
print(f'  Baseline data spans: {bperp_dates[0]} ... {bperp_dates[-1]}')
print(f'  Processing {n_total} interferograms ...')

prog_bar  = ptime.progressBar(maxValue=n_total)
n_ok      = 0
n_skip    = 0

for i, unw_file in enumerate(unw_files):
    ifg_dir = os.path.dirname(unw_file)
    pair    = _dates_from_dir(ifg_dir)

    if pair is None:
        print(f'\n  WARNING: Cannot extract dates from: {os.path.basename(ifg_dir)}')
        n_skip += 1
        prog_bar.update(i + 1)
        continue

    date1, date2 = pair

    missing = [d for d in (date1, date2) if d not in bperp_dict]
    if missing:
        print(f'\n  WARNING: Missing baseline for {missing} in pair {date1}_{date2}')
        n_skip += 1
        prog_bar.update(i + 1)
        continue

    ifg_meta = dict(common_meta)
    ifg_meta.update(read_tif_spatial_meta(unw_file))

    ifg_meta['DATE12'] = f'{ptime.yymmdd(date1)}-{ptime.yymmdd(date2)}'

    bperp = float(bperp_dict[date2]) - float(bperp_dict[date1])
    ifg_meta['P_BASELINE_TOP_HDR']    = str(bperp)
    ifg_meta['P_BASELINE_BOTTOM_HDR'] = str(bperp)

    rsc_file = unw_file + '.rsc'
    writefile.write_roipac_rsc(
        ifg_meta, rsc_file,
        update_mode=update_mode, print_msg=False,
    )

    n_ok += 1
    prog_bar.update(i + 1, suffix=f'{date1}_{date2}')

prog_bar.close()
print(f'\n  Summary: {n_ok} written, {n_skip} skipped (of {n_total} total)')

if n_ok == 0:
    raise RuntimeError('No interferograms were successfully processed!')

##############################################################

-------------------- Main entry point ------------------

##############################################################

EXAMPLE = """Example usage:
prep_licsar.py smallbaselineApp.cfg
prep_licsar.py smallbaselineApp.cfg --no-update
prep_licsar.py -t smallbaselineApp.cfg

Required template keys (add to your .cfg):
mintpy.load.processor = licsar

GMTSAR-format baseline table (single file)

mintpy.load.baselineDir = /data/baseline_table.dat

Interferograms (TIF, one per directory named YYYYMMDD_YYYYMMDD)

mintpy.load.unwFile = /data/ifg//filt_fine.unw.tif
mintpy.load.corFile = /data/ifg//filt_fine.cor.tif

Geometry files (GeoTIFF)

mintpy.load.demFile = /data/geometry/dem.tif
mintpy.load.incAngleFile = /data/geometry/incLocal.tif
mintpy.load.azAngleFile = /data/geometry/az.tif
"""

def cmd_line_parse(iargs=None):
parser = argparse.ArgumentParser(
description='Prepare .rsc sidecar files for a LiCSAR/SARscape '
'GeoTIFF interferogram stack.',
formatter_class=argparse.RawDescriptionHelpFormatter,
epilog=EXAMPLE,
)
parser.add_argument(
'template_file',
nargs='?',
default=None,
help='Path to the MintPy template (.cfg) file',
)
parser.add_argument(
'-t', '--template',
dest='template_file_opt',
default=None,
help='Path to the MintPy template (.cfg) file (alternative to positional)',
)
parser.add_argument(
'--no-update',
dest='update_mode',
action='store_false',
default=True,
help='Force re-write of existing .rsc files',
)
inps = parser.parse_args(args=iargs)

# 合并两种传入方式：-t 优先，其次位置参数
if inps.template_file_opt:
    inps.template_file = inps.template_file_opt

if inps.template_file is None:
    parser.error('Please provide the template .cfg file.')

return inps

def _get_template_value(template, key):
"""Get a string value from template, stripping quotes and ignoring auto/False/None."""
val = template.get(key, '')
if not val:
return ''
val = str(val).strip("'"")
if val.lower() in ('auto', 'none', 'false', 'true', ''):
return ''
return val

def prep_licsar(inps):
"""Main processing pipeline."""
print('=' * 60)
print('prep_licsar: preparing MintPy-compatible sidecar files')
print('=' * 60)

# ------------------------------------------------------------------ #
# 0. Read template                                                     #
# ------------------------------------------------------------------ #
template_raw = readfile.read_template(inps.template_file)
template     = ut.check_template_auto_value(template_raw)

# check_template_auto_value 可能把路径键清空，从原始模板恢复
for key in ['mintpy.load.baselineDir', 'mintpy.load.metaFile',
            'mintpy.load.unwFile',     'mintpy.load.corFile',
            'mintpy.load.demFile',     'mintpy.load.incAngleFile',
            'mintpy.load.azAngleFile']:
    if not template.get(key) or str(template.get(key)).lower() in ('auto', 'none', 'false'):
        raw_val = template_raw.get(key, '')
        if raw_val and str(raw_val).lower() not in ('auto', 'none', ''):
            template[key] = raw_val

# ------------------------------------------------------------------ #
# 1. Resolve file paths from template                                  #
# ------------------------------------------------------------------ #
prm_file    = _get_template_value(template, 'mintpy.load.metaFile')
bperp_file  = _get_template_value(template, 'mintpy.load.baselineDir')
unw_pattern = _get_template_value(template, 'mintpy.load.unwFile')

print(f'\nPRM file         : {prm_file!r}')
print(f'Baseline file    : {bperp_file!r}')
print(f'Unwrapped pattern: {unw_pattern!r}')

# metaFile 是可选的
if prm_file and not os.path.isfile(prm_file):
    print(f'  WARNING: mintpy.load.metaFile not found: {prm_file!r}, '
          f'will read metadata from GeoTIFF')
    prm_file = ''

# baselineDir 是必须的
if not bperp_file:
    raise FileNotFoundError(
        'mintpy.load.baselineDir is not set in the template.\n'
        'Please add: mintpy.load.baselineDir = /path/to/baseline_table.dat')
if not os.path.isfile(bperp_file):
    raise FileNotFoundError(
        f'Baseline file not found: {bperp_file!r}\n'
        'Please check mintpy.load.baselineDir in the template.')

# unwFile 是必须的
if not unw_pattern:
    raise FileNotFoundError(
        'mintpy.load.unwFile is not set in the template.')
unw_files = sorted(glob.glob(unw_pattern))
if not unw_files:
    raise FileNotFoundError(f'No unwrapped files matched: {unw_pattern!r}')

print(f'Unwrapped files  : {len(unw_files)} found')

# ------------------------------------------------------------------ #
# 2. Build common (stack-level) metadata                              #
# ------------------------------------------------------------------ #
print('\n--- Building common metadata ---')

common_meta = read_tif_spatial_meta(unw_files[0])
if prm_file:
    prm_meta = extract_sarscape_metadata(prm_file)
    common_meta.update(prm_meta)

orbit_dir = common_meta.get('ORBIT_DIRECTION', 'ascending')
corners   = get_image_corners(unw_files[0], orbit_dir)
common_meta.update(corners)

common_meta = compute_incidence_angle(common_meta)

if 'X_FIRST' in common_meta:
    x_first = float(common_meta['X_FIRST'])
    if x_first > 180.0:
        common_meta['X_FIRST'] = str(x_first - 360.0)

for k, v in list(common_meta.items()):
    common_meta[k] = str(v)

common_meta = readfile.standardize_metadata(common_meta)

defaults = {
    'WAVELENGTH':      '0.05546576',
    'EARTH_RADIUS':    '6371000.0',
    'HEIGHT':          '693000.0',
    'PLATFORM':        'Sen',
    'PROCESSOR':       'licsar',
    'ORBIT_DIRECTION': 'ascending',
}
for k, v in defaults.items():
    if k not in common_meta:
        common_meta[k] = v

print('\nKey metadata summary:')
for k in ['PLATFORM', 'WAVELENGTH', 'HEIGHT', 'ORBIT_DIRECTION',
          'PROCESSOR', 'EARTH_RADIUS', 'INCIDENCE_ANGLE']:
    if k in common_meta:
        print(f'  {k:25s} = {common_meta[k]}')

stack_rsc = os.path.join(os.path.dirname(unw_files[0]), '_stack.rsc')
print(f'\nWriting stack-level RSC: {stack_rsc}')
os.makedirs(os.path.dirname(os.path.abspath(stack_rsc)), exist_ok=True)
writefile.write_roipac_rsc(common_meta, stack_rsc)

# ------------------------------------------------------------------ #
# 3. Baseline table                                                    #
# ------------------------------------------------------------------ #
print('\n--- Reading GMTSAR baseline table ---')
bperp_dict = read_baseline_table(bperp_file)

# ------------------------------------------------------------------ #
# 4. Geometry files                                                    #
# ------------------------------------------------------------------ #
print('\n--- Preparing geometry sidecar files ---')
prepare_geometry(template, common_meta, update_mode=inps.update_mode)

# ------------------------------------------------------------------ #
# 5. Interferogram stack                                               #
# ------------------------------------------------------------------ #
print('\n--- Preparing interferogram sidecar files ---')
prepare_stack(unw_files, common_meta, bperp_dict,
              update_mode=inps.update_mode)

print('\n' + '=' * 60)
print('prep_licsar: Done.')
print('=' * 60)

##############################################################

-------------------- CLI entry -------------------------

##############################################################

def main(iargs=None):
inps = cmd_line_parse(iargs)
prep_licsar(inps)

if name == 'main':
main(sys.argv[1:])

avgPhaseVelocity.png

avgSpatialCoh.png

these pngs sem to be normal

temporalCoherence.png

this png seems to be strange.

Sorry to bother you, but I’d really appreciate any help anyone can offer.

[welcome]

👋 Thanks for opening your first issue here! Please filled out the template with as much details as possible. We appreciate that you took the time to contribute!
Make sure you read our contributing guidelines.