From eec888271c90cdc971558722e17e409bbc966538 Mon Sep 17 00:00:00 2001 From: Bekaert Date: Mon, 2 Mar 2026 09:43:26 +0000 Subject: [PATCH 01/24] add basic entries to enable OPERA tropo support for parsing in MintPy --- src/mintpy/cli/tropo_opera.py | 109 ++++++++++++++ src/mintpy/tropo_opera.py | 265 ++++++++++++++++++++++++++++++++++ 2 files changed, 374 insertions(+) create mode 100755 src/mintpy/cli/tropo_opera.py create mode 100644 src/mintpy/tropo_opera.py diff --git a/src/mintpy/cli/tropo_opera.py b/src/mintpy/cli/tropo_opera.py new file mode 100755 index 000000000..be0d91468 --- /dev/null +++ b/src/mintpy/cli/tropo_opera.py @@ -0,0 +1,109 @@ +#!/usr/bin/env python3 +############################################################ +# Program is part of MintPy # +# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi # +# Author: David Bekaert, March 2026 # +############################################################ + + +import os +import sys + +from mintpy.utils.arg_utils import create_argument_parser + +############################################################################ +REFERENCE = """references: + in prerp to be added. +""" + +DIR_DEMO = """--dir ./OPERA + To be added + ... +""" + +EXAMPLE = """example: + tropo_opera.py -f timeseries.h5 -g inputs/geometryRadar.h5 --dir ./OPERA + tropo_opera.py -f geo/geo_timeseries.h5 -g geo/geo_geometryRadar.h5 --dir ./OPERA +""" + + +def create_parser(subparsers=None): + synopsis = 'Tropospheric correction using OPERA computed Zenith Tropopsheric Delays from ECMWF HRES data (https://www.earthdata.nasa.gov/data/catalog/asf-opera-l4-tropo-zenith-v1-1)' + epilog = REFERENCE + '\n' + DIR_DEMO + '\n' + EXAMPLE + name = __name__.split('.')[-1] + parser = create_argument_parser( + name, synopsis=synopsis, description=synopsis, epilog=epilog, subparsers=subparsers) + parser.add_argument('-f', '--file', dest='dis_file', required=True, + help='timeseries HDF5 file, i.e. timeseries.h5') + parser.add_argument('-g', '--geom', dest='geom_file', required=True, + help='geometry file.') + parser.add_argument('--dir','--opera-dir', dest='opera_dir', default='./OPERA', + help='directory to store downloaded OPERA delays data (default: %(default)s).') + parser.add_argument('-o', dest='cor_dis_file', + help='Output file name for trospheric corrected timeseries.') + + return parser + + +def cmd_line_parse(iargs=None): + """Command line parser.""" + # parse + parser = create_parser() + inps = parser.parse_args(args=iargs) + + # import + from mintpy.utils import readfile + + # check: --opera-dir (use absolute path) + inps.opera_dir = os.path.abspath(inps.opera_dir) + print('Use OPERA products at directory:', inps.opera_dir) + + # check: exsitence of input files + for fname in [inps.dis_file, inps.geom_file]: + if fname and not os.path.isfile(fname): + raise FileNotFoundError(f'input file not exist: {fname}') + + # check: processors & coordinates of input files + atr1 = readfile.read_attribute(inps.dis_file) + atr2 = readfile.read_attribute(inps.geom_file) + coord1 = 'geo' if 'Y_FIRST' in atr1.keys() else 'radar' + coord2 = 'geo' if 'Y_FIRST' in atr2.keys() else 'radar' + proc = atr1.get('PROCESSOR', 'isce') + + # check: radar-coord product from gamma and roipac is not supported + if coord1 == 'radar' and proc in ['gamma', 'roipac']: + msg = f'Radar-coded file from {proc} is NOT supported!' + msg += '\n Try to geocode the time-series and geometry files and re-run with them instead.' + raise ValueError(msg) + + # check: coordinate system must be consistent btw. displacement and geometry files + if coord1 != coord2: + n = max(len(os.path.basename(i)) for i in [inps.dis_file, inps.geom_file]) + msg = 'Input time-series and geometry file are NOT in the same coordinate!' + msg += f'\n file {os.path.basename(inps.dis_file):<{n}} coordinate: {coord1}' + msg += f'\n file {os.path.basename(inps.geom_file):<{n}} coordinate: {coord2}' + raise ValueError(msg) + + # default: output tropo and corrected displacement file names + inps.tropo_file = os.path.join(os.path.dirname(inps.geom_file), 'OPERA.h5') + if not inps.cor_dis_file: + inps.cor_dis_file = inps.dis_file.split('.')[0] + '_OPERA.h5' + + return inps + + +############################################################################ +def main(iargs=None): + # parse + inps = cmd_line_parse(iargs) + + # import + from mintpy.tropo_opera import run_tropo_opera + + # run + run_tropo_opera(inps) + + +############################################################################ +if __name__ == '__main__': + main(sys.argv[1:]) diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py new file mode 100644 index 000000000..1b66f1ee6 --- /dev/null +++ b/src/mintpy/tropo_opera.py @@ -0,0 +1,265 @@ +############################################################ +# Program is part of MintPy # +# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi # +# Author: Sara Mirzaee, Zhang Yunjun, Bhuvan Varugu, 2018 # +############################################################ + + +import os +import re + +import h5py +import numpy as np +from scipy.interpolate import RegularGridInterpolator +from skimage.transform import resize + +import mintpy.cli.diff +from mintpy.objects import timeseries +from mintpy.utils import ptime, readfile, utils as ut, writefile + + +############################################################################ +def get_delay_geo(ztd_file, atr, cos_inc_angle): + """calc single path tropo delay in line-of-sight direction + + Parameters: ztd_file - str, path of zenith delay file + atr - dict, dictionary of attribute for output file + cos_inc_angle - 2D np.ndarray in float32, cos(inc_angle) + Returns: delay - 2D np.ndarray in float32, LOS delay + """ + + # get geo_box from ts_file attributes + length, width = int(atr['LENGTH']), int(atr['WIDTH']) + geo_box = ut.coordinate(atr).box_pixel2geo(pixel_box=(0, 0, width, length)) + + # geo_box --> pix_box in ztd file + atr_ztd = readfile.read_attribute(ztd_file) + pix_box = ut.coordinate(atr_ztd).box_geo2pixel(geo_box) + + # read ztd file + delay = readfile.read(ztd_file, box=pix_box)[0] + + # interpolate/resample into the same resolution as ts_file + delay = resize( + delay, (length, width), + order=1, + mode='constant', + anti_aliasing=True, + preserve_range=True, + ) + + # project from zenith to line-of-sight + delay /= cos_inc_angle + + # reverse the sign for consistency between different phase correction steps/methods + delay *= -1 + + return delay + + +def get_delay_radar(ztd_file, cos_inc_angle, pts_new): + """calc single path tropo delay in line-of-sight direction + + Parameters: ztd_file - str, path of zenith delay file + cos_inc_angle - 2D np.ndarray in (len, wid) in float32, cos(inc_angle) + pts_new - 2D np.ndarray in (len*wid, 2) in float32 + Returns: delay - 2D np.ndarray in float32, LOS delay + """ + # read ztd file + delay_ztd, atr_ztd = readfile.read(ztd_file) + # flip to be consistent with the reversed lats + delay_ztd = np.flipud(delay_ztd) + + # pixel coordinates in ztd file + lats, lons = ut.get_lat_lon(atr_ztd, dimension=1) + # set lats in ascending order as required by RegularGridInterpolator + lats = np.flipud(lats) + pts_ztd = ((lats.flatten(), + lons.flatten())) + + # resample in pts_new coordinates + interp_func = RegularGridInterpolator( + pts_ztd, + delay_ztd, + method='nearest', + bounds_error=False, + fill_value=0, + ) + delay = interp_func(pts_new) + delay = delay.reshape(cos_inc_angle.shape) + + # project from zenith to line-of-sight + delay /= cos_inc_angle + + # reverse the sign for consistency between different phase correction steps/methods + delay *= -1 + + return delay + + +############################################################################ +def calculate_delay_timeseries(tropo_file, dis_file, geom_file, opera_dir): + """calculate delay time-series and write to HDF5 file""" + + ## get list of dates + atr = readfile.read_attribute(dis_file) + ftype = atr['FILE_TYPE'] + if ftype == 'timeseries': + date_list = timeseries(dis_file).get_date_list() + + elif ftype == '.unw': + date12 = readfile.read_attribute(dis_file)['DATE12'] + date_list = ptime.yyyymmdd(date12.split('-')) + + else: + raise ValueError(f'un-supported displacement file type: {ftype}') + + # list of dates --> list of ztd files + ztd_files = [] + flag = np.ones(len(date_list), dtype=np.bool_) + for i, date_str in enumerate(date_list): + fnames = [os.path.join(opera_dir, f'{date_str}{fext}') for fext in ['.ztd', '.ztd.tif']] + fnames = [f for f in fnames if os.path.exists(f)] + if len(fnames) > 0: + ztd_files.append(fnames[0]) + else: + print(f'WARNING: NO ztd file found for {date_str}! Continue without it.') + flag[i] = False + + # update date_list to be consistent with ztd_files + if np.any(flag == 0): + date_list = np.array(date_list)[flag].tolist() + + ## update_mode + def get_dataset_size(fname): + atr = readfile.read_attribute(fname) + return (atr['LENGTH'], atr['WIDTH']) + + def run_or_skip(ztd_files, tropo_file, geom_file): + print('update mode: ON') + print(f'output file: {tropo_file}') + flag = 'skip' + + # check existence and modification time + if ut.run_or_skip(out_file=tropo_file, in_file=ztd_files, print_msg=False) == 'run': + flag = 'run' + print('1) output file either do NOT exist or is NOT newer than all ZTD files.') + + else: + print('1) output file exists and is newer than all ZTD files.') + + # check dataset size in space / time + date_list = [str(re.findall(r'\d{8}', i)[0]) for i in ztd_files] + if (get_dataset_size(tropo_file) != get_dataset_size(geom_file) + or any(i not in timeseries(tropo_file).get_date_list() for i in date_list)): + flag = 'run' + print(('2) output file does NOT have the same len/wid as the geometry file {}' + ' or does NOT contain all dates').format(geom_file)) + else: + print('2) output file has the same len/wid as the geometry file and contains all dates') + + # check if output file is fully written + with h5py.File(tropo_file, 'r') as f: + if np.all(f['timeseries'][-1,:,:] == 0): + flag = 'run' + print('3) output file is NOT fully written.') + else: + print('3) output file is fully written.') + + # result + print(f'run or skip: {flag}') + return flag + + if run_or_skip(ztd_files, tropo_file, geom_file) == 'skip': + return + + + ## prepare output file + + # metadata + atr['FILE_TYPE'] = 'timeseries' + atr['UNIT'] = 'm' + + # remove metadata related with double reference + # because absolute delay is calculated and saved + for key in ['REF_DATE', 'REF_X', 'REF_Y', 'REF_LAT', 'REF_LON']: + if key in atr.keys(): + atr.pop(key) + + # instantiate time-series + length, width = int(atr['LENGTH']), int(atr['WIDTH']) + num_date = len(date_list) + dates = np.array(date_list, dtype=np.bytes_) + ds_name_dict = { + "date" : [dates.dtype, (num_date,), dates], + "timeseries" : [np.float32, (num_date, length, width), None], + } + writefile.layout_hdf5(tropo_file, ds_name_dict, metadata=atr) + + + ## calculate phase delay + + # read geometry + print(f'read incidenceAngle from file: {geom_file}') + inc_angle = readfile.read(geom_file, datasetName='incidenceAngle')[0] + cos_inc_angle = np.cos(inc_angle * np.pi / 180.0) + + if 'Y_FIRST' in atr.keys(): + # No need for data in geo-coordinates + pts_new = None + + else: + # Get pixel lat/lon for data in radar-coordinates + print('get pixel coordinates in geometry file') + lats, lons = ut.get_lat_lon(atr, geom_file) + pts_new = np.hstack((lats.reshape(-1, 1), + lons.reshape(-1, 1))) + + # loop for date-by-date IO + prog_bar = ptime.progressBar(maxValue=num_date) + for i in range(num_date): + date_str = date_list[i] + ztd_file = ztd_files[i] + + # calc delay + if 'Y_FIRST' in atr.keys(): + delay = get_delay_geo(ztd_file, atr, cos_inc_angle) + + else: + delay = get_delay_radar(ztd_file, cos_inc_angle, pts_new) + + # write delay to file + block = [i, i+1, 0, length, 0, width] + writefile.write_hdf5_block( + tropo_file, + data=delay, + datasetName='timeseries', + block=block, + print_msg=False, + ) + + prog_bar.update(i + 1, suffix=os.path.basename(ztd_file)) + prog_bar.close() + + return tropo_file + + +############################################################################ +def run_tropo_opera(inps): + + # calculate tropo delay and savee to h5 file + calculate_delay_timeseries( + tropo_file=inps.tropo_file, + dis_file=inps.dis_file, + geom_file=inps.geom_file, + opera_dir=inps.opera_dir) + + # correct tropo delay (using diff.py) + # diff.py can handle different reference in space and time + # e.g. the absolute delay and the double referenced time-series + print('correcting delay for using diff.py') + iargs = [inps.dis_file, inps.tropo_file, '-o', inps.cor_dis_file, '--force'] + print('diff.py', ' '.join(iargs)) + mintpy.cli.diff.main(iargs) + + return From dc0bcdb5d973eb7861fc3151ffbc476a98c127ca Mon Sep 17 00:00:00 2001 From: David Date: Wed, 4 Mar 2026 16:59:48 +0100 Subject: [PATCH 02/24] first attempt on tropo_opera.py --- src/mintpy/tropo_opera.py | 787 ++++++++++++++++++++++++++++++-------- 1 file changed, 619 insertions(+), 168 deletions(-) diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index 1b66f1ee6..b065b6f7e 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -1,244 +1,632 @@ ############################################################ # Program is part of MintPy # # Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi # -# Author: Sara Mirzaee, Zhang Yunjun, Bhuvan Varugu, 2018 # +# Author: David Bekaert, March 2026 # ############################################################ +import glob +import netrc import os import re +from datetime import datetime, timedelta import h5py import numpy as np from scipy.interpolate import RegularGridInterpolator -from skimage.transform import resize import mintpy.cli.diff from mintpy.objects import timeseries from mintpy.utils import ptime, readfile, utils as ut, writefile -############################################################################ -def get_delay_geo(ztd_file, atr, cos_inc_angle): - """calc single path tropo delay in line-of-sight direction +OPERA_MODEL_HOURS = (0, 6, 12, 18) +OPERA_MIN_ACQ_DATE = '20160101' +OPERA_TROPO_COLLECTIONS = [ + 'C3717139408-ASF', + 'C1273910987-ASF', +] - Parameters: ztd_file - str, path of zenith delay file - atr - dict, dictionary of attribute for output file - cos_inc_angle - 2D np.ndarray in float32, cos(inc_angle) - Returns: delay - 2D np.ndarray in float32, LOS delay - """ - # get geo_box from ts_file attributes - length, width = int(atr['LENGTH']), int(atr['WIDTH']) - geo_box = ut.coordinate(atr).box_pixel2geo(pixel_box=(0, 0, width, length)) +def _model_time_token(date_str, hour_str): + """Build OPERA model timestamp token as YYYYMMDDTHHMMSSZ.""" + return f'{date_str}T{int(hour_str):02d}0000Z' - # geo_box --> pix_box in ztd file - atr_ztd = readfile.read_attribute(ztd_file) - pix_box = ut.coordinate(atr_ztd).box_geo2pixel(geo_box) - # read ztd file - delay = readfile.read(ztd_file, box=pix_box)[0] +def _is_valid_opera_file(opera_file): + """Return True if OPERA file is readable and has required datasets.""" + try: + with h5py.File(opera_file, 'r') as fobj: + req_dsets = ['latitude', 'longitude', 'height', 'wet_delay', 'hydrostatic_delay'] + return all(k in fobj for k in req_dsets) + except Exception: + return False - # interpolate/resample into the same resolution as ts_file - delay = resize( - delay, (length, width), - order=1, - mode='constant', - anti_aliasing=True, - preserve_range=True, - ) - # project from zenith to line-of-sight - delay /= cos_inc_angle +def get_geom_lat_lon_bounds(geom_file): + """Get geometry file lat/lon bounds in degrees.""" + from mintpy.tropo_pyaps3 import get_bounding_box + + meta = readfile.read_attribute(geom_file) + lat0, lat1, lon0, lon1 = get_bounding_box(meta, geom_file=geom_file) + return min(lat0, lat1), max(lat0, lat1), min(lon0, lon1), max(lon0, lon1) + + +def get_opera_crop_indices(lat, lon, geom_file, pad_cells=3): + """Get OPERA cube crop indices from geometry bounds with fixed ±0.2 deg padding. + + The pad_cells argument is kept for backward compatibility and ignored. + """ + if lat.size < 2 or lon.size < 2: + raise ValueError('Invalid OPERA latitude/longitude arrays: size < 2') + + lat_min, lat_max, lon_min, lon_max = get_geom_lat_lon_bounds(geom_file) + + buffer_deg = 0.2 + lat_min_pad = lat_min - buffer_deg + lat_max_pad = lat_max + buffer_deg + lon_min_pad = lon_min - buffer_deg + lon_max_pad = lon_max + buffer_deg + + lat_idx = np.where((lat >= lat_min_pad) & (lat <= lat_max_pad))[0] + lon_idx = np.where((lon >= lon_min_pad) & (lon <= lon_max_pad))[0] + if lat_idx.size == 0 or lon_idx.size == 0: + raise ValueError('No overlap between OPERA grid and MintPy geometry bounds') + + i0, i1 = int(lat_idx[0]), int(lat_idx[-1]) + 1 + j0, j1 = int(lon_idx[0]), int(lon_idx[-1]) + 1 + + crop_info = { + 'lat_min': lat_min, + 'lat_max': lat_max, + 'lon_min': lon_min, + 'lon_max': lon_max, + 'lat_min_pad': lat_min_pad, + 'lat_max_pad': lat_max_pad, + 'lon_min_pad': lon_min_pad, + 'lon_max_pad': lon_max_pad, + 'buffer_deg': buffer_deg, + } + return (i0, i1, j0, j1), crop_info - # reverse the sign for consistency between different phase correction steps/methods - delay *= -1 - return delay +def read_opera_total_delay_cube(opera_file, geom_file, pad_cells=3): + """Read and crop OPERA delay cube, returning total zenith delay. + Returns a dict with keys: + latitude, longitude, height, wet_delay, hydro_delay, total_delay, crop_info + """ + with h5py.File(opera_file, 'r') as fobj: + req_dsets = ['latitude', 'longitude', 'height', 'wet_delay', 'hydrostatic_delay'] + missing = [k for k in req_dsets if k not in fobj] + if missing: + raise ValueError(f'OPERA file missing expected datasets: {missing}') + + lat = fobj['latitude'][:] + lon = fobj['longitude'][:] + height = fobj['height'][:] + + (i0, i1, j0, j1), crop_info = get_opera_crop_indices(lat, lon, geom_file, pad_cells=pad_cells) + + lat_crop = lat[i0:i1] + lon_crop = lon[j0:j1] + + wet = fobj['wet_delay'][0, :, i0:i1, j0:j1] + hydro = fobj['hydrostatic_delay'][0, :, i0:i1, j0:j1] + total = wet + hydro + + return { + 'latitude': lat_crop, + 'longitude': lon_crop, + 'height': height, + 'wet_delay': wet, + 'hydro_delay': hydro, + 'total_delay': total, + 'crop_info': crop_info, + } -def get_delay_radar(ztd_file, cos_inc_angle, pts_new): - """calc single path tropo delay in line-of-sight direction - Parameters: ztd_file - str, path of zenith delay file - cos_inc_angle - 2D np.ndarray in (len, wid) in float32, cos(inc_angle) - pts_new - 2D np.ndarray in (len*wid, 2) in float32 - Returns: delay - 2D np.ndarray in float32, LOS delay +def get_geom_lat_lon_dem(geom_file): + """Read DEM and latitude/longitude grids from geometry file. + + Returns: lat2d/lon2d/dem in 2D np.ndarray with same shape """ - # read ztd file - delay_ztd, atr_ztd = readfile.read(ztd_file) - # flip to be consistent with the reversed lats - delay_ztd = np.flipud(delay_ztd) - - # pixel coordinates in ztd file - lats, lons = ut.get_lat_lon(atr_ztd, dimension=1) - # set lats in ascending order as required by RegularGridInterpolator - lats = np.flipud(lats) - pts_ztd = ((lats.flatten(), - lons.flatten())) - - # resample in pts_new coordinates - interp_func = RegularGridInterpolator( - pts_ztd, - delay_ztd, - method='nearest', + dem = readfile.read(geom_file, datasetName='height')[0].astype(np.float32) + atr = readfile.read_attribute(geom_file) + + dset_list = readfile.get_dataset_list(geom_file) + if 'latitude' in dset_list and 'longitude' in dset_list: + lat2d = readfile.read(geom_file, datasetName='latitude')[0].astype(np.float64) + lon2d = readfile.read(geom_file, datasetName='longitude')[0].astype(np.float64) + elif 'Y_FIRST' in atr and 'X_FIRST' in atr: + length, width = int(atr['LENGTH']), int(atr['WIDTH']) + lat0 = float(atr['Y_FIRST']) + lon0 = float(atr['X_FIRST']) + dlat = float(atr['Y_STEP']) + dlon = float(atr['X_STEP']) + lat1d = lat0 + dlat * np.arange(length, dtype=np.float64) + lon1d = lon0 + dlon * np.arange(width, dtype=np.float64) + lon2d, lat2d = np.meshgrid(lon1d, lat1d) + else: + raise ValueError('Unable to get latitude/longitude grids from geometry file!') + + return lat2d, lon2d, dem + + +def calc_zenith_delay_from_opera_file(opera_file, geom_file, pad_cells=3): + """Calculate 2D zenith tropospheric delay map intersected with DEM.""" + cube = read_opera_total_delay_cube(opera_file, geom_file, pad_cells=pad_cells) + lat2d, lon2d, dem = get_geom_lat_lon_dem(geom_file) + + z_axis = np.asarray(cube['height'], dtype=np.float64) + y_axis = np.asarray(cube['latitude'], dtype=np.float64) + x_axis = np.asarray(cube['longitude'], dtype=np.float64) + data = np.asarray(cube['total_delay'], dtype=np.float32) + + # RegularGridInterpolator requires strictly ascending axes. + if z_axis[1] < z_axis[0]: + z_axis = z_axis[::-1] + data = data[::-1, :, :] + if y_axis[1] < y_axis[0]: + y_axis = y_axis[::-1] + data = data[:, ::-1, :] + if x_axis[1] < x_axis[0]: + x_axis = x_axis[::-1] + data = data[:, :, ::-1] + + interp = RegularGridInterpolator( + (z_axis, y_axis, x_axis), + data, + method='cubic', bounds_error=False, - fill_value=0, + fill_value=np.nan, ) - delay = interp_func(pts_new) - delay = delay.reshape(cos_inc_angle.shape) - # project from zenith to line-of-sight - delay /= cos_inc_angle + points = np.column_stack([ + dem.reshape(-1).astype(np.float64), + lat2d.reshape(-1).astype(np.float64), + lon2d.reshape(-1).astype(np.float64), + ]) + ztd = interp(points).reshape(dem.shape).astype(np.float32) - # reverse the sign for consistency between different phase correction steps/methods - delay *= -1 + # mask invalid pixels from geometry + ztd[~np.isfinite(dem)] = np.nan + ztd[dem == 0] = np.nan - return delay + return ztd, cube ############################################################################ -def calculate_delay_timeseries(tropo_file, dis_file, geom_file, opera_dir): - """calculate delay time-series and write to HDF5 file""" +def read_inps2date_time(inps): + """Read acquisition date/time info from input arguments. - ## get list of dates - atr = readfile.read_attribute(dis_file) + Parameters: inps - Namespace for input arguments + Returns: date_list - list(str), acquisition dates in YYYYMMDD + utc_sec - float, acquisition UTC time in seconds + """ + if not inps.dis_file: + raise ValueError('input displacement file is required to derive OPERA date/time info.') + + print(f'read date/time info from file: {inps.dis_file}') + atr = readfile.read_attribute(inps.dis_file) ftype = atr['FILE_TYPE'] - if ftype == 'timeseries': - date_list = timeseries(dis_file).get_date_list() + if ftype == 'timeseries': + date_list = timeseries(inps.dis_file).get_date_list() elif ftype == '.unw': - date12 = readfile.read_attribute(dis_file)['DATE12'] - date_list = ptime.yyyymmdd(date12.split('-')) - + date_list = ptime.yyyymmdd(atr['DATE12'].split('-')) else: raise ValueError(f'un-supported displacement file type: {ftype}') - # list of dates --> list of ztd files - ztd_files = [] - flag = np.ones(len(date_list), dtype=np.bool_) - for i, date_str in enumerate(date_list): - fnames = [os.path.join(opera_dir, f'{date_str}{fext}') for fext in ['.ztd', '.ztd.tif']] - fnames = [f for f in fnames if os.path.exists(f)] - if len(fnames) > 0: - ztd_files.append(fnames[0]) - else: - print(f'WARNING: NO ztd file found for {date_str}! Continue without it.') - flag[i] = False + if 'CENTER_LINE_UTC' not in atr: + raise ValueError(f'CENTER_LINE_UTC is missing in metadata of file: {inps.dis_file}') + + utc_sec = float(atr['CENTER_LINE_UTC']) + + inps.date_list = date_list + inps.utc_sec = utc_sec + + return date_list, utc_sec + + +def nearest_opera_product_time(date_str, utc_sec): + """Round one acquisition date/time to the nearest OPERA model datetime. + + Parameters: date_str - str, acquisition date in YYYYMMDD + utc_sec - float, acquisition UTC time in seconds + Returns: out_date - str, OPERA model date in YYYYMMDD + out_hour - str, OPERA model hour in HH format + """ + acq_dt = datetime.strptime(date_str, '%Y%m%d') + timedelta(seconds=float(utc_sec)) + + candidates = [] + for day_offset in (-1, 0, 1): + base_day = (acq_dt + timedelta(days=day_offset)).replace(hour=0, minute=0, second=0, microsecond=0) + for hour in OPERA_MODEL_HOURS: + cand_dt = base_day + timedelta(hours=hour) + candidates.append(cand_dt) + + opera_dt = min(candidates, key=lambda x: (abs(x - acq_dt), x)) + out_date = opera_dt.strftime('%Y%m%d') + out_hour = opera_dt.strftime('%H') + return out_date, out_hour + + +def get_opera_date_time_list(date_list, utc_sec): + """Build required OPERA date/hour list for all acquisitions. + + Parameters: date_list - list(str), acquisition dates in YYYYMMDD + utc_sec - float, acquisition UTC time in seconds + Returns: opera_date_list - list(str), OPERA model date in YYYYMMDD + opera_hour_list - list(str), OPERA model hour in HH format + """ + opera_date_list = [] + opera_hour_list = [] + for date_str in date_list: + opera_date, opera_hour = nearest_opera_product_time(date_str, utc_sec) + opera_date_list.append(opera_date) + opera_hour_list.append(opera_hour) + + return opera_date_list, opera_hour_list + + +def get_expected_opera_file_patterns(opera_date_list, opera_hour_list): + """Create expected OPERA ZTD filename patterns for model times. + + The production timestamp is left as wildcard. + + Returns: expected_patterns - list(str), e.g. + OPERA_L4_TROPO-ZENITH_YYYYMMDDTHHMMSSZ_*_HRES_v* + """ + expected_patterns = sorted({ + f'OPERA_L4_TROPO-ZENITH_{_model_time_token(date_str, hour_str)}_*_HRES_v*' + for date_str, hour_str in zip(opera_date_list, opera_hour_list) + }) + return expected_patterns + + +def get_opera_file_status(opera_date_list, opera_hour_list, opera_dir): + """Get matched OPERA files and missing model date/hour list. + + Returns: expected_patterns - list(str) + matched_files - dict(str, list(str)) + missing_date_hour_list - list(tuple(str, str)) in (YYYYMMDD, HH) + """ + expected_patterns = get_expected_opera_file_patterns(opera_date_list, opera_hour_list) + + matched_files = {} + for pattern in expected_patterns: + files = sorted(glob.glob(os.path.join(opera_dir, pattern))) + valid_files = [f for f in files if _is_valid_opera_file(f)] + bad_files = [f for f in files if f not in valid_files] + if len(bad_files) > 0: + print(f'WARNING: ignore corrupted/invalid OPERA files for pattern {pattern}:') + for f in bad_files: + print(f' {os.path.basename(f)}') + matched_files[pattern] = valid_files + + missing_list = [] + seen = set() + for date_str, hour_str in zip(opera_date_list, opera_hour_list): + pattern = f'OPERA_L4_TROPO-ZENITH_{_model_time_token(date_str, hour_str)}_*_HRES_v*' + if len(matched_files.get(pattern, [])) == 0: + key = (date_str, f'{int(hour_str):02d}') + if key not in seen: + missing_list.append(key) + seen.add(key) + return expected_patterns, matched_files, missing_list + + +def _asf_product_text(product): + """Get searchable text from ASF product object.""" + text_items = [str(product)] + + for key in ['sceneName', 'fileName', 'granuleName', 'url', 'downloadUrl', 'fileID']: + val = getattr(product, key, None) + if val: + text_items.append(str(val)) + + props = getattr(product, 'properties', None) + if isinstance(props, dict): + for key in ['sceneName', 'fileName', 'granuleName', 'url', 'downloadUrl', 'fileID']: + val = props.get(key, None) + if val: + text_items.append(str(val)) + + add_urls = props.get('additionalUrls', None) + if isinstance(add_urls, (list, tuple)): + for val in add_urls: + if val: + text_items.append(str(val)) + + return ' '.join(text_items) + + +def _parse_opera_time_tokens(text): + """Parse model/production time tokens from OPERA filename-like text.""" + mobj = re.search( + r'OPERA_L4_TROPO-ZENITH_(\d{8}T\d{6}Z)_(\d{8}T\d{6}Z)_HRES_v', + text, + ) + if not mobj: + return None, None + return mobj.group(1), mobj.group(2) + + +def _create_asf_session(asf): + """Create authenticated ASF session using ~/.netrc when possible.""" + if not hasattr(asf, 'ASFSession'): + return None + + session = asf.ASFSession() + + # Newer asf_search API + if hasattr(session, 'auth_with_netrc'): + session.auth_with_netrc() + return session + + # Backward-compatible fallback + if hasattr(session, 'auth_with_creds'): + host = 'urs.earthdata.nasa.gov' + parsed = netrc.netrc().authenticators(host) + if not parsed: + raise ValueError(f'No authentication credentials found in ~/.netrc for {host}') + + username = parsed[0] + password = parsed[2] + session.auth_with_creds(username, password) + return session + + return None + - # update date_list to be consistent with ztd_files - if np.any(flag == 0): - date_list = np.array(date_list)[flag].tolist() +def dload_opera_files(missing_date_hour_list, opera_dir): + """Download missing OPERA TROPO-ZENITH files via ASF Search. + + Uses per-day ASF queries constrained by OPERA collection IDs and + date windows, then down-filters by required model-time tokens. + """ + if len(missing_date_hour_list) == 0: + return + + try: + import asf_search as asf + except Exception as exc: + raise ImportError( + 'Missing dependency: asf_search. Install it to enable OPERA download.' + ) from exc + + os.makedirs(opera_dir, exist_ok=True) + + print('-' * 50) + print('search/download missing OPERA TROPO-ZENITH products from ASF ...') + + # Assume user has ~/.netrc configured. Create session once when available. + session = None + try: + session = _create_asf_session(asf) + if session is not None: + print(' ASF authentication: enabled (from ~/.netrc)') + except Exception as exc: + print(f' WARNING: ASF authentication setup failed: {exc}') + session = None + + missing_tokens = sorted({_model_time_token(d, h) for d, h in missing_date_hour_list}) + date_to_tokens = {} + for token in missing_tokens: + date_to_tokens.setdefault(token[:8], []).append(token) + + n_selected = 0 + prog_bar = ptime.progressBar(maxValue=len(missing_tokens), prefix=' download ') + + print(' token match summary (token: candidates -> selected):') + n_done = 0 + for date_str in sorted(date_to_tokens.keys()): + day_tokens = date_to_tokens[date_str] + start_str = f'{date_str[0:4]}-{date_str[4:6]}-{date_str[6:8]}T00:00:00Z' + end_str = f'{date_str[0:4]}-{date_str[4:6]}-{date_str[6:8]}T23:59:59Z' + + opts = asf.ASFSearchOptions(**{ + 'maxResults': 250, + 'collections': OPERA_TROPO_COLLECTIONS, + 'collectionAlias': False, + 'start': start_str, + 'end': end_str, + }) + if session is not None: + opts.session = session + + results = asf.search(opts=opts) + + token_best_product = {} + token_best_prod_time = {} + token_match_count = {token: 0 for token in day_tokens} + + for product in results: + text = _asf_product_text(product) + if 'OPERA_L4_TROPO-ZENITH_' not in text: + continue + + model_token, prod_token = _parse_opera_time_tokens(text) + if model_token not in token_match_count: + continue + + token_match_count[model_token] += 1 + new_prod_time = prod_token or '' + old_prod_time = token_best_prod_time.get(model_token, '') + if model_token not in token_best_product or new_prod_time > old_prod_time: + token_best_product[model_token] = product + token_best_prod_time[model_token] = new_prod_time + + for token in day_tokens: + if token not in token_best_product: + n_done += 1 + prog_bar.update(n_done, suffix=f'{token} (no match)') + print(f' {token}: {token_match_count[token]} -> no') + continue + + dload_results = asf.ASFSearchResults([token_best_product[token]]) + if session is not None: + dload_results.download(path=opera_dir, session=session) + else: + dload_results.download(path=opera_dir) + + n_selected += 1 + n_done += 1 + prog_bar.update(n_done, suffix=f'{token} (downloaded)') + print(f' {token}: {token_match_count[token]} -> yes') + + prog_bar.close() + + print(f' missing model-time tokens downloaded : {n_selected} / {len(missing_tokens)}') + print(f' downloaded files to: {opera_dir}') + print('-' * 50) + + +############################################################################ +def _get_best_local_opera_file(date_str, hour_str, opera_dir): + token = _model_time_token(date_str, hour_str) + pattern = os.path.join(opera_dir, f'OPERA_L4_TROPO-ZENITH_{token}_*_HRES_v*.nc') + files = sorted(f for f in glob.glob(pattern) if _is_valid_opera_file(f)) + if len(files) == 0: + return None + return files[-1] + + +def _project_zenith_to_los(zenith_delay, inc_angle_deg): + cos_inc = np.cos(np.deg2rad(inc_angle_deg.astype(np.float64))) + out = zenith_delay.astype(np.float32).copy() + invalid = (~np.isfinite(cos_inc)) | (np.abs(cos_inc) < 1.0e-8) + out[invalid] = np.nan + out[~invalid] = out[~invalid] / cos_inc[~invalid] + out *= -1 + return out + + +def _split_supported_dates(date_list, min_date=OPERA_MIN_ACQ_DATE): + """Split acquisition dates into OPERA-supported and unsupported sets.""" + supported = [d for d in date_list if d >= min_date] + unsupported = [d for d in date_list if d < min_date] + return supported, unsupported + + +def _run_or_skip_opera(opera_files, used_dates, tropo_file, geom_file): + """Run/skip decision for OPERA delay file generation (MintPy-style update mode).""" - ## update_mode def get_dataset_size(fname): atr = readfile.read_attribute(fname) return (atr['LENGTH'], atr['WIDTH']) - def run_or_skip(ztd_files, tropo_file, geom_file): - print('update mode: ON') - print(f'output file: {tropo_file}') - flag = 'skip' + print('update mode: ON') + print(f'output file: {tropo_file}') + flag = 'skip' - # check existence and modification time - if ut.run_or_skip(out_file=tropo_file, in_file=ztd_files, print_msg=False) == 'run': - flag = 'run' - print('1) output file either do NOT exist or is NOT newer than all ZTD files.') + if ut.run_or_skip(out_file=tropo_file, in_file=opera_files, print_msg=False) == 'run': + flag = 'run' + print('1) output file either does NOT exist or is NOT newer than all OPERA files.') + else: + print('1) output file exists and is newer than all OPERA files.') + if (get_dataset_size(tropo_file) != get_dataset_size(geom_file) + or any(i not in timeseries(tropo_file).get_date_list() for i in used_dates)): + flag = 'run' + print(('2) output file does NOT have the same len/wid as the geometry file {}' + ' or does NOT contain all required dates').format(geom_file)) else: - print('1) output file exists and is newer than all ZTD files.') - - # check dataset size in space / time - date_list = [str(re.findall(r'\d{8}', i)[0]) for i in ztd_files] - if (get_dataset_size(tropo_file) != get_dataset_size(geom_file) - or any(i not in timeseries(tropo_file).get_date_list() for i in date_list)): - flag = 'run' - print(('2) output file does NOT have the same len/wid as the geometry file {}' - ' or does NOT contain all dates').format(geom_file)) - else: - print('2) output file has the same len/wid as the geometry file and contains all dates') + print('2) output file has the same len/wid as the geometry file and contains all required dates.') + with h5py.File(tropo_file, 'r') as fobj: + if np.all(fobj['timeseries'][-1, :, :] == 0): + flag = 'run' + print('3) output file is NOT fully written.') + else: + print('3) output file is fully written.') - # check if output file is fully written - with h5py.File(tropo_file, 'r') as f: - if np.all(f['timeseries'][-1,:,:] == 0): - flag = 'run' - print('3) output file is NOT fully written.') - else: - print('3) output file is fully written.') + print(f'run or skip: {flag}') + return flag - # result - print(f'run or skip: {flag}') - return flag - if run_or_skip(ztd_files, tropo_file, geom_file) == 'skip': - return +############################################################################ +def calculate_delay_timeseries(tropo_file, dis_file, geom_file, opera_dir): + """Calculate OPERA delay time-series and write to HDF5 file.""" + atr = readfile.read_attribute(dis_file) + ftype = atr['FILE_TYPE'] + if ftype == 'timeseries': + date_list = timeseries(dis_file).get_date_list() + elif ftype == '.unw': + date12 = atr['DATE12'] + date_list = ptime.yyyymmdd(date12.split('-')) + else: + raise ValueError(f'un-supported displacement file type: {ftype}') - ## prepare output file + if 'CENTER_LINE_UTC' not in atr: + raise ValueError(f'CENTER_LINE_UTC is missing in metadata of file: {dis_file}') + utc_sec = float(atr['CENTER_LINE_UTC']) - # metadata - atr['FILE_TYPE'] = 'timeseries' - atr['UNIT'] = 'm' + supported_dates, unsupported_dates = _split_supported_dates(date_list, min_date=OPERA_MIN_ACQ_DATE) + if len(unsupported_dates) > 0: + print((f'WARNING: {len(unsupported_dates)} acquisition date(s) are earlier than ' + f'OPERA availability ({OPERA_MIN_ACQ_DATE}) and will be skipped.')) + for d in unsupported_dates: + print(f' {d}') - # remove metadata related with double reference - # because absolute delay is calculated and saved - for key in ['REF_DATE', 'REF_X', 'REF_Y', 'REF_LAT', 'REF_LON']: - if key in atr.keys(): - atr.pop(key) + if len(supported_dates) == 0: + raise RuntimeError( + f'No acquisition dates are supported by OPERA (requires >= {OPERA_MIN_ACQ_DATE}).' + ) - # instantiate time-series - length, width = int(atr['LENGTH']), int(atr['WIDTH']) - num_date = len(date_list) - dates = np.array(date_list, dtype=np.bytes_) - ds_name_dict = { - "date" : [dates.dtype, (num_date,), dates], - "timeseries" : [np.float32, (num_date, length, width), None], - } - writefile.layout_hdf5(tropo_file, ds_name_dict, metadata=atr) + model_dates, model_hours = get_opera_date_time_list(supported_dates, utc_sec) + opera_files = [] + used_dates = [] + for acq_date, model_date, model_hour in zip(supported_dates, model_dates, model_hours): + fname = _get_best_local_opera_file(model_date, model_hour, opera_dir) + if fname is None: + print(f'WARNING: NO OPERA file found for acquisition {acq_date} -> {model_date} {model_hour}:00') + continue + used_dates.append(acq_date) + opera_files.append(fname) - ## calculate phase delay + if len(opera_files) == 0: + raise RuntimeError('No OPERA files found for requested acquisitions. Cannot create delay time-series.') - # read geometry - print(f'read incidenceAngle from file: {geom_file}') - inc_angle = readfile.read(geom_file, datasetName='incidenceAngle')[0] - cos_inc_angle = np.cos(inc_angle * np.pi / 180.0) + if _run_or_skip_opera(opera_files, used_dates, tropo_file, geom_file) == 'skip': + return tropo_file - if 'Y_FIRST' in atr.keys(): - # No need for data in geo-coordinates - pts_new = None + atr_out = atr.copy() + atr_out['FILE_TYPE'] = 'timeseries' + atr_out['UNIT'] = 'm' + for key in ['REF_DATE', 'REF_X', 'REF_Y', 'REF_LAT', 'REF_LON']: + if key in atr_out: + atr_out.pop(key) - else: - # Get pixel lat/lon for data in radar-coordinates - print('get pixel coordinates in geometry file') - lats, lons = ut.get_lat_lon(atr, geom_file) - pts_new = np.hstack((lats.reshape(-1, 1), - lons.reshape(-1, 1))) + length, width = int(atr_out['LENGTH']), int(atr_out['WIDTH']) + num_date = len(opera_files) + dates = np.array(used_dates, dtype=np.bytes_) + ds_name_dict = { + 'date': [dates.dtype, (num_date,), dates], + 'timeseries': [np.float32, (num_date, length, width), None], + } + writefile.layout_hdf5(tropo_file, ds_name_dict, metadata=atr_out) + + print(f'read incidenceAngle from file: {geom_file}') + inc_angle = readfile.read(geom_file, datasetName='incidenceAngle')[0].astype(np.float32) - # loop for date-by-date IO prog_bar = ptime.progressBar(maxValue=num_date) for i in range(num_date): - date_str = date_list[i] - ztd_file = ztd_files[i] + ztd, _ = calc_zenith_delay_from_opera_file(opera_files[i], geom_file, pad_cells=3) + delay_los = _project_zenith_to_los(ztd, inc_angle) - # calc delay - if 'Y_FIRST' in atr.keys(): - delay = get_delay_geo(ztd_file, atr, cos_inc_angle) - - else: - delay = get_delay_radar(ztd_file, cos_inc_angle, pts_new) - - # write delay to file - block = [i, i+1, 0, length, 0, width] + block = [i, i + 1, 0, length, 0, width] writefile.write_hdf5_block( tropo_file, - data=delay, + data=delay_los, datasetName='timeseries', block=block, print_msg=False, ) - prog_bar.update(i + 1, suffix=os.path.basename(ztd_file)) + prog_bar.update(i + 1, suffix=os.path.basename(opera_files[i])) prog_bar.close() return tropo_file @@ -246,17 +634,80 @@ def run_or_skip(ztd_files, tropo_file, geom_file): ############################################################################ def run_tropo_opera(inps): + """Run OPERA-based tropospheric correction workflow.""" + print('tropospheric delay correction with OPERA approach') + print(f'input displacement file : {inps.dis_file}') + print(f'input geometry file : {inps.geom_file}') + print(f'input OPERA directory : {inps.opera_dir}') + print(f'output tropo file : {inps.tropo_file}') + print(f'output corrected file : {inps.cor_dis_file}') + + date_list, utc_sec = read_inps2date_time(inps) + supported_dates, unsupported_dates = _split_supported_dates(date_list, min_date=OPERA_MIN_ACQ_DATE) + if len(unsupported_dates) > 0: + print((f'WARNING: {len(unsupported_dates)} acquisition date(s) are earlier than ' + f'OPERA availability ({OPERA_MIN_ACQ_DATE}) and will be skipped.')) + for d in unsupported_dates: + print(f' {d}') + + if len(supported_dates) == 0: + print(f'No acquisition dates >= {OPERA_MIN_ACQ_DATE}. Skip OPERA correction.') + return + + opera_date_list, opera_hour_list = get_opera_date_time_list(supported_dates, utc_sec) + expected_patterns, matched_files, missing_date_hour_list = get_opera_file_status( + opera_date_list=opera_date_list, + opera_hour_list=opera_hour_list, + opera_dir=inps.opera_dir, + ) + + inps.opera_date_list = opera_date_list + inps.opera_hour_list = opera_hour_list + inps.expected_opera_patterns = expected_patterns + inps.matched_opera_files = matched_files + inps.missing_opera_date_hour_list = missing_date_hour_list + + if len(missing_date_hour_list) > 0: + dload_opera_files(missing_date_hour_list, inps.opera_dir) + expected_patterns, matched_files, missing_date_hour_list = get_opera_file_status( + opera_date_list=opera_date_list, + opera_hour_list=opera_hour_list, + opera_dir=inps.opera_dir, + ) + inps.expected_opera_patterns = expected_patterns + inps.matched_opera_files = matched_files + inps.missing_opera_date_hour_list = missing_date_hour_list + + print(f'number of acquisitions (input) : {len(date_list)}') + print(f'number of acquisitions (supported): {len(supported_dates)}') + print(f'acquisition UTC seconds: {utc_sec}') + if len(set(opera_hour_list)) == 1: + print(f'nearest OPERA model time: {opera_hour_list[0]}:00 UTC') + else: + print(f'nearest OPERA model hours: {sorted(set(opera_hour_list))}') + + n_total = len(expected_patterns) + n_found = sum(len(i) > 0 for i in matched_files.values()) + n_missing = len(missing_date_hour_list) + print(f'expected OPERA model files (unique): {n_total}') + print(f'expected OPERA model files found : {n_found}') + print(f'OPERA model files still missing : {n_missing}') + if n_missing > 0: + print('missing model date/hour list (YYYYMMDD, HH):') + for date_str, hour_str in missing_date_hour_list: + print(f' {date_str}, {hour_str}') + + if n_found == 0: + print('No OPERA files available to calculate delay. Stop.') + return - # calculate tropo delay and savee to h5 file calculate_delay_timeseries( tropo_file=inps.tropo_file, dis_file=inps.dis_file, geom_file=inps.geom_file, - opera_dir=inps.opera_dir) + opera_dir=inps.opera_dir, + ) - # correct tropo delay (using diff.py) - # diff.py can handle different reference in space and time - # e.g. the absolute delay and the double referenced time-series print('correcting delay for using diff.py') iargs = [inps.dis_file, inps.tropo_file, '-o', inps.cor_dis_file, '--force'] print('diff.py', ' '.join(iargs)) From 0283efae2288b496d02f8f6562b228a18fc82ac4 Mon Sep 17 00:00:00 2001 From: David Date: Wed, 4 Mar 2026 19:33:41 +0100 Subject: [PATCH 03/24] add opera functions to enable its method for tropo correction --- pyproject.toml | 1 + src/mintpy/__main__.py | 8 ++++++++ src/mintpy/defaults/smallbaselineApp.cfg | 7 ++++++- src/mintpy/defaults/smallbaselineApp_auto.cfg | 3 +++ src/mintpy/smallbaselineApp.py | 13 +++++++++++++ src/mintpy/workflow/__init__.py | 1 + 6 files changed, 32 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index e0cacae14..22880a5cf 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -96,6 +96,7 @@ Issues = "https://github.com/insarlab/MintPy/issues" "timeseries2velocity.py" = "mintpy.cli.timeseries2velocity:main" "timeseries_rms.py" = "mintpy.cli.timeseries_rms:main" "tropo_gacos.py" = "mintpy.cli.tropo_gacos:main" +"tropo_opera.py" = "mintpy.cli.tropo_opera:main" "tropo_phase_elevation.py" = "mintpy.cli.tropo_phase_elevation:main" "tropo_pyaps3.py" = "mintpy.cli.tropo_pyaps3:main" "tsview.py" = "mintpy.cli.tsview:main" diff --git a/src/mintpy/__main__.py b/src/mintpy/__main__.py index f6b3e6fc8..ded28f330 100644 --- a/src/mintpy/__main__.py +++ b/src/mintpy/__main__.py @@ -497,6 +497,13 @@ def get_tropo_gacos_parser(subparsers=None): return parser +def get_tropo_opera_parser(subparsers=None): + from mintpy.cli import tropo_opera + parser = tropo_opera.create_parser(subparsers) + parser.set_defaults(func=tropo_opera.main) + return parser + + def get_tropo_phase_elevation_parser(subparsers=None): from mintpy.cli import tropo_phase_elevation parser = tropo_phase_elevation.create_parser(subparsers) @@ -629,6 +636,7 @@ def get_parser(): get_s1ab_range_bias_parser(sp) get_solid_earth_tides_parser(sp) get_tropo_gacos_parser(sp) + get_tropo_opera_parser(sp) get_tropo_phase_elevation_parser(sp) get_tropo_pyaps3_parser(sp) get_unwrap_error_bridging_parser(sp) diff --git a/src/mintpy/defaults/smallbaselineApp.cfg b/src/mintpy/defaults/smallbaselineApp.cfg index 4d8269768..76e883959 100644 --- a/src/mintpy/defaults/smallbaselineApp.cfg +++ b/src/mintpy/defaults/smallbaselineApp.cfg @@ -224,7 +224,8 @@ mintpy.ionosphericDelay.excludeDate12 = auto #[20080520_20090817 / no], auto fo ## NARR - NARR from NOAA [need to install PyAPS from Caltech/EarthDef; recommended for N America] ## c. gacos - use GACOS with the iterative tropospheric decomposition model (Yu et al., 2018, JGR) ## need to manually download GACOS products at http://www.gacos.net for all acquisitions before running this step -mintpy.troposphericDelay.method = auto #[pyaps / height_correlation / gacos / no], auto for pyaps +## d. opera - use OPERA L4 tropospheric zenith delay products +mintpy.troposphericDelay.method = auto #[pyaps / height_correlation / gacos / opera / no], auto for pyaps ## Notes for pyaps: ## a. GAM data latency: with the most recent SAR data, there will be GAM data missing, the correction @@ -249,6 +250,10 @@ mintpy.troposphericDelay.minCorrelation = auto #[0.0-1.0], auto for 0 ## Set the path below to directory that contains the downloaded *.ztd* files mintpy.troposphericDelay.gacosDir = auto # [path2directory], auto for "./GACOS" +## Notes for opera: +## Set the path below to directory that contains downloaded OPERA *.nc files +mintpy.troposphericDelay.operaDir = auto # [path2directory], auto for "./OPERA" + ########## 9. deramp (optional) ## Estimate and remove a phase ramp for each acquisition based on the reliable pixels. diff --git a/src/mintpy/defaults/smallbaselineApp_auto.cfg b/src/mintpy/defaults/smallbaselineApp_auto.cfg index 6763d77bc..255d2397a 100644 --- a/src/mintpy/defaults/smallbaselineApp_auto.cfg +++ b/src/mintpy/defaults/smallbaselineApp_auto.cfg @@ -107,6 +107,9 @@ mintpy.troposphericDelay.minCorrelation = 0 ## gacos mintpy.troposphericDelay.gacosDir = ./GACOS +## opera +mintpy.troposphericDelay.operaDir = ./OPERA + ########## deramp mintpy.deramp = no diff --git a/src/mintpy/smallbaselineApp.py b/src/mintpy/smallbaselineApp.py index 9e607093a..4d17e7805 100644 --- a/src/mintpy/smallbaselineApp.py +++ b/src/mintpy/smallbaselineApp.py @@ -502,6 +502,9 @@ def get_timeseries_filename(template, work_dir='./'): elif method == 'gacos': fname1 = f'{os.path.splitext(fname0)[0]}_GACOS.h5' + elif method == 'opera': + fname1 = f'{os.path.splitext(fname0)[0]}_OPERA.h5' + elif method == 'pyaps': fname1 = f'{os.path.splitext(fname0)[0]}_{model}.h5' @@ -648,6 +651,16 @@ def get_dataset_size(fname): import mintpy.cli.tropo_phase_elevation mintpy.cli.tropo_phase_elevation.main(iargs) + # OPERA ZTD products derived from ECMWF HRES data and available through the ASF DAAC + elif method == 'opera': + opera_dir = self.template['mintpy.troposphericDelay.operaDir'] + iargs = ['-f', in_file, '-g', geom_file, '-o', out_file, '--dir', opera_dir] + print('tropospheric delay correction with OPERA ZTD products') + print('\ntropo_opera.py', ' '.join(iargs)) + if ut.run_or_skip(out_file=out_file, in_file=in_file) == 'run': + import mintpy.cli.tropo_opera + mintpy.cli.tropo_opera.main(iargs) + # Weather re-analysis data with iterative tropospheric decomposition (GACOS) # Yu et al. (2018, JGR) elif method == 'gacos': diff --git a/src/mintpy/workflow/__init__.py b/src/mintpy/workflow/__init__.py index 6957f3f0e..50d00ed49 100644 --- a/src/mintpy/workflow/__init__.py +++ b/src/mintpy/workflow/__init__.py @@ -31,6 +31,7 @@ 'timeseries2velocity', 'timeseries_rms', 'tropo_gacos', + 'tropo_opera', 'tropo_phase_elevation', 'tropo_pyaps3', 'unwrap_error_bridging', From 16411d268112ac543de80f1dd9cfa5a601c1f659 Mon Sep 17 00:00:00 2001 From: dbekaert Date: Thu, 5 Mar 2026 00:13:22 +0100 Subject: [PATCH 04/24] Update src/mintpy/tropo_opera.py Co-authored-by: sourcery-ai[bot] <58596630+sourcery-ai[bot]@users.noreply.github.com> --- src/mintpy/tropo_opera.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index b065b6f7e..f19282f55 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -189,9 +189,8 @@ def calc_zenith_delay_from_opera_file(opera_file, geom_file, pad_cells=3): ]) ztd = interp(points).reshape(dem.shape).astype(np.float32) - # mask invalid pixels from geometry + # mask invalid pixels from DEM (NaN/Inf); preserve valid zero-elevation pixels ztd[~np.isfinite(dem)] = np.nan - ztd[dem == 0] = np.nan return ztd, cube From c059ae8d8c0ccfdf01feb975cbd0cf0fdc63fae3 Mon Sep 17 00:00:00 2001 From: David Bekaert Date: Tue, 31 Mar 2026 10:05:04 +0200 Subject: [PATCH 05/24] Added subsetting instead of full file download. This speeds up drastically --- src/mintpy/tropo_opera.py | 394 ++++++++++++++++++++++++++++++++++---- 1 file changed, 360 insertions(+), 34 deletions(-) diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index f19282f55..17417c80c 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -90,9 +90,57 @@ def get_opera_crop_indices(lat, lon, geom_file, pad_cells=3): return (i0, i1, j0, j1), crop_info -def read_opera_total_delay_cube(opera_file, geom_file, pad_cells=3): +def get_opera_height_crop_indices(height_levels, dem_min, dem_max): + """Get height level indices spanning the DEM range. + + The vertical interpolation is linear, so only the levels that bracket + [dem_min, dem_max] are needed — no extra buffer beyond that. + + For example, if dem_min=-5, dem_max=10 and height levels are + [-20, -10, 0, 10, 20], the returned slice covers [-10, 0, 10] + (the first level at-or-below dem_min through the first level + at-or-above dem_max). + + Parameters: height_levels - 1D np.ndarray, OPERA height levels (metres) + dem_min - float, minimum DEM elevation (metres) + dem_max - float, maximum DEM elevation (metres) + Returns: k0, k1 - int, start/end slice indices into *height_levels* + """ + h = np.asarray(height_levels, dtype=np.float64) + n = len(h) + if n < 2: + return 0, n + + ascending = h[1] > h[0] + if not ascending: + h = h[::-1] + + # first level at or below dem_min + k_low = int(np.searchsorted(h, dem_min, side='right')) - 1 + k_low = max(0, k_low) + + # first level at or above dem_max + k_high = int(np.searchsorted(h, dem_max, side='left')) + k_high = min(n - 1, k_high) + + k0, k1 = k_low, k_high + 1 # +1 for Python slice end + + if not ascending: + k0, k1 = n - k1, n - k0 + + return k0, k1 + + +def read_opera_total_delay_cube(opera_file, geom_file, dem_range=None, pad_cells=3): """Read and crop OPERA delay cube, returning total zenith delay. + Parameters: opera_file - str, path to OPERA netCDF file + geom_file - str, path to MintPy geometry HDF5 file + dem_range - tuple(float, float) or None, (dem_min, dem_max) in + metres. When provided the height dimension is also + cropped to the levels spanning [dem_min, dem_max] + plus a one-level buffer on each side. + pad_cells - int, kept for backward compatibility (ignored) Returns a dict with keys: latitude, longitude, height, wet_delay, hydro_delay, total_delay, crop_info """ @@ -108,17 +156,24 @@ def read_opera_total_delay_cube(opera_file, geom_file, pad_cells=3): (i0, i1, j0, j1), crop_info = get_opera_crop_indices(lat, lon, geom_file, pad_cells=pad_cells) + # vertical subsetting + if dem_range is not None: + k0, k1 = get_opera_height_crop_indices(height, dem_range[0], dem_range[1]) + else: + k0, k1 = 0, len(height) + lat_crop = lat[i0:i1] lon_crop = lon[j0:j1] + height_crop = height[k0:k1] - wet = fobj['wet_delay'][0, :, i0:i1, j0:j1] - hydro = fobj['hydrostatic_delay'][0, :, i0:i1, j0:j1] + wet = fobj['wet_delay'][0, k0:k1, i0:i1, j0:j1] + hydro = fobj['hydrostatic_delay'][0, k0:k1, i0:i1, j0:j1] total = wet + hydro return { 'latitude': lat_crop, 'longitude': lon_crop, - 'height': height, + 'height': height_crop, 'wet_delay': wet, 'hydro_delay': hydro, 'total_delay': total, @@ -154,40 +209,155 @@ def get_geom_lat_lon_dem(geom_file): def calc_zenith_delay_from_opera_file(opera_file, geom_file, pad_cells=3): - """Calculate 2D zenith tropospheric delay map intersected with DEM.""" - cube = read_opera_total_delay_cube(opera_file, geom_file, pad_cells=pad_cells) + """Calculate 2D zenith tropospheric delay map intersected with DEM. + + Two-step interpolation: + 1. Linear interpolation in the vertical (height -> DEM elevation) + at each OPERA lat/lon grid node. + 2. Cubic interpolation in the lateral (lat/lon -> pixel locations). + """ lat2d, lon2d, dem = get_geom_lat_lon_dem(geom_file) + # derive DEM range for vertical subsetting + valid_dem = dem[np.isfinite(dem)] + if valid_dem.size > 0: + dem_range = (float(np.nanmin(valid_dem)), float(np.nanmax(valid_dem))) + else: + dem_range = None + + cube = read_opera_total_delay_cube(opera_file, geom_file, dem_range=dem_range, pad_cells=pad_cells) + z_axis = np.asarray(cube['height'], dtype=np.float64) y_axis = np.asarray(cube['latitude'], dtype=np.float64) x_axis = np.asarray(cube['longitude'], dtype=np.float64) - data = np.asarray(cube['total_delay'], dtype=np.float32) + data = np.asarray(cube['total_delay'], dtype=np.float64) - # RegularGridInterpolator requires strictly ascending axes. - if z_axis[1] < z_axis[0]: + # Ensure strictly ascending axes for interpolation. + if z_axis.size > 1 and z_axis[1] < z_axis[0]: z_axis = z_axis[::-1] data = data[::-1, :, :] - if y_axis[1] < y_axis[0]: + if y_axis.size > 1 and y_axis[1] < y_axis[0]: y_axis = y_axis[::-1] data = data[:, ::-1, :] - if x_axis[1] < x_axis[0]: + if x_axis.size > 1 and x_axis[1] < x_axis[0]: x_axis = x_axis[::-1] data = data[:, :, ::-1] - interp = RegularGridInterpolator( - (z_axis, y_axis, x_axis), - data, - method='cubic', - bounds_error=False, - fill_value=np.nan, - ) - - points = np.column_stack([ - dem.reshape(-1).astype(np.float64), - lat2d.reshape(-1).astype(np.float64), - lon2d.reshape(-1).astype(np.float64), + # --- Step 1: Linear interpolation in the vertical ----------------------- + # At each OPERA (lat, lon) grid node, interpolate along height to the + # DEM elevation of every output pixel. This produces a 2-D slab per + # (lat_node, lon_node) → we collect them into a 2-D grid (ny, nx) for + # each output pixel, then do step 2. + # + # Efficient approach: for every output pixel, the DEM elevation is known. + # Build a 1-D linear interpolator per (lat, lon) column, evaluate at the + # pixel's DEM elevation → gives delay_2d[iy, ix] on the OPERA lat/lon grid. + # Then do cubic 2-D interpolation of that 2-D field to the pixel lat/lon. + + nrows, ncols = dem.shape + dem_flat = dem.ravel().astype(np.float64) + lat_flat = lat2d.ravel().astype(np.float64) + lon_flat = lon2d.ravel().astype(np.float64) + + ny, nx = len(y_axis), len(x_axis) + + # For each OPERA (lat, lon) column, build a 1-D linear interpolator + # along height, then evaluate at *all* output pixel DEM values at once. + # Result: delay_at_dem[j, i, pixel] but that's memory-heavy. + # Instead, use RegularGridInterpolator in 1-D (height) per column → too slow. + # + # Better: use scipy interp1d broadcast or manual linear interp. + # Use np.interp per column vectorised over the output pixels. + # + # Most efficient: do the vertical interp as a single + # RegularGridInterpolator(method='linear') in 3-D, then re-interpolate + # the result in 2-D with cubic. But that couples axes again. + # + # Cleanest two-step: collapse vertical first on the OPERA grid, then + # interpolate laterally. + # For each output pixel p with DEM elevation h_p: + # delay_on_opera_grid[iy, ix] = interp1d(z_axis, data[:, iy, ix])(h_p) + # Then: ztd[p] = interp2d(y_axis, x_axis, delay_on_opera_grid)(lat_p, lon_p) + # + # This is O(npixels * ny * nx) if done naively. + # Vectorise: for all pixels at once, linear interp along z for all columns. + + # Pre-compute vertical interpolation weights for each pixel + # np.interp works on 1-D; we need to interp data[:, iy, ix] at dem values. + # Reshape data to (nz, ny*nx), interp each column at all dem values → large. + # Instead: interp all columns at each unique dem value. + + # Practical approach: loop over OPERA grid columns (ny*nx is small, ~hundreds) + # and use np.interp for all pixels at once per column. + # Result: ztd_on_grid shape (ny, nx, npixels) → then for each pixel, do 2-D interp. + # That's still large. Better: for each pixel, do 1-D vertical interp at + # 1 height value across all columns, giving (ny, nx), then 2-D interp. + # But looping over pixels is slow. + + # Most practical vectorised approach: + # Step 1: vertical linear interp → 2D field on OPERA grid for each pixel + # Use broadcasting: find bracketing indices and weights once per pixel. + idx = np.searchsorted(z_axis, dem_flat, side='right') - 1 + idx = np.clip(idx, 0, len(z_axis) - 2) + # fractional weight + dz = z_axis[idx + 1] - z_axis[idx] + dz[dz == 0] = 1.0 # avoid division by zero + w = (dem_flat - z_axis[idx]) / dz + w = np.clip(w, 0.0, 1.0) + + # data shape: (nz, ny, nx) + # For each pixel p: delay_2d[:, :, p] = data[idx[p], :, :] * (1-w[p]) + data[idx[p]+1, :, :] * w[p] + # Vectorise: gather the two bracketing slabs for all pixels + # data_lo[p, iy, ix] = data[idx[p], iy, ix] → shape (npixels, ny, nx) + npixels = len(dem_flat) + data_lo = data[idx] # shape (npixels, ny, nx) + data_hi = data[idx + 1] # shape (npixels, ny, nx) + # Weighted blend: shape (npixels, ny, nx) + w3 = w[:, np.newaxis, np.newaxis] + delay_2d = data_lo * (1.0 - w3) + data_hi * w3 # (npixels, ny, nx) + + # --- Step 2: Cubic interpolation in lat/lon per pixel ------------------- + # For each pixel, we have a 2-D delay field on the OPERA (y_axis, x_axis) + # grid. Interpolate to the pixel's (lat, lon) using cubic. + # + # Batch approach: all pixels share the same OPERA grid, and each pixel has + # its own 2-D field. Use RegularGridInterpolator per pixel → too slow. + # + # Faster: the lateral interpolation weights depend only on (lat, lon) and + # are the same for all pixels at the same location. Use + # RegularGridInterpolator once with a dummy z-axis to vectorise: + # Build a 3-D field (npixels, ny, nx) with "pixel index" as the first + # axis → but that's not a regular grid in the first axis. + # + # Best practical approach: pre-compute cubic interpolation weights for + # each pixel's (lat, lon) on the OPERA grid, then apply them. + # scipy.interpolate doesn't expose weights directly, so use map_coordinates + # or manual cubic. + # + # Simple and efficient: use scipy.ndimage.map_coordinates with order=3 + # (cubic) on each pixel's 2-D field. But looping over pixels is slow. + # + # Actually, since all pixels share the OPERA grid, we can convert + # (lat, lon) to fractional grid indices once, then use them for all pixels. + + from scipy.ndimage import map_coordinates + + # Convert pixel lat/lon to fractional indices in the OPERA grid + # y_axis is ascending, x_axis is ascending + iy_frac = np.interp(lat_flat, y_axis, np.arange(ny)) + ix_frac = np.interp(lon_flat, x_axis, np.arange(nx)) + + # map_coordinates on delay_2d: shape (npixels, ny, nx) + # We want to evaluate delay_2d[p, iy_frac[p], ix_frac[p]] for each p. + # Coordinates: axis0=p (exact integer), axis1=iy_frac, axis2=ix_frac + coords = np.array([ + np.arange(npixels, dtype=np.float64), + iy_frac, + ix_frac, ]) - ztd = interp(points).reshape(dem.shape).astype(np.float32) + ztd_flat = map_coordinates(delay_2d, coords, order=3, mode='nearest') + + ztd = ztd_flat.reshape(dem.shape).astype(np.float32) # mask invalid pixels from DEM (NaN/Inf); preserve valid zero-elevation pixels ztd[~np.isfinite(dem)] = np.nan @@ -379,11 +549,128 @@ def _create_asf_session(asf): return None -def dload_opera_files(missing_date_hour_list, opera_dir): +def _get_product_url(product): + """Extract the data download URL from an ASF search product.""" + # Try the most common attribute paths + url = getattr(product, 'url', None) + if url: + return str(url) + + props = getattr(product, 'properties', None) + if isinstance(props, dict): + for key in ('url', 'downloadUrl'): + val = props.get(key) + if val: + return str(val) + + return None + + +def _get_product_filename(product): + """Extract the filename from an ASF search product.""" + props = getattr(product, 'properties', None) + if isinstance(props, dict): + fname = props.get('fileName') + if fname: + return str(fname) + + # Fallback: parse from text representation + text = _asf_product_text(product) + mobj = re.search(r'(OPERA_L4_TROPO-ZENITH_\S+\.nc)', text) + if mobj: + return mobj.group(1) + + return None + + +def dload_opera_file_subset(product, opera_dir, geom_file, dem_range=None, session=None): + """Download a spatially and vertically subsetted OPERA file. + + Uses ``fsspec`` + ``h5py`` to open the remote file via byte-range HTTP + requests, reads only the required spatial/vertical slice, and writes the + subset to a local HDF5 file that is compatible with + :func:`read_opera_total_delay_cube`. + + Parameters: product - ASF search product object + opera_dir - str, local directory for output files + geom_file - str, path to MintPy geometry file (for bounding box) + dem_range - tuple(float, float) or None, (min, max) DEM in metres + session - ASF session object (used for authentication cookies) + Returns: out_file - str, path to the written local subset file, or None + """ + import fsspec + + url = _get_product_url(product) + if url is None: + print(' WARNING: could not extract download URL from product') + return None + + filename = _get_product_filename(product) + if filename is None: + print(' WARNING: could not extract filename from product') + return None + + out_file = os.path.join(opera_dir, filename) + + # Build authenticated fsspec HTTP options from ASF session cookies + fs_kwargs = {} + if session is not None: + cookies = getattr(session, 'cookies', None) + if cookies is not None: + cookie_str = '; '.join(f'{k}={v}' for k, v in dict(cookies).items()) + fs_kwargs['headers'] = {'Cookie': cookie_str} + + storage_opts = {'https': fs_kwargs, 'http': fs_kwargs} + scheme = 'https' if url.startswith('https') else 'http' + + with fsspec.open(url, mode='rb', **storage_opts.get(scheme, {})) as f: + with h5py.File(f, 'r') as src: + lat = src['latitude'][:] + lon = src['longitude'][:] + height = src['height'][:] + + (i0, i1, j0, j1), _crop_info = get_opera_crop_indices(lat, lon, geom_file) + + if dem_range is not None: + k0, k1 = get_opera_height_crop_indices(height, dem_range[0], dem_range[1]) + else: + k0, k1 = 0, len(height) + + lat_crop = lat[i0:i1] + lon_crop = lon[j0:j1] + height_crop = height[k0:k1] + + wet_crop = src['wet_delay'][0, k0:k1, i0:i1, j0:j1] + hydro_crop = src['hydrostatic_delay'][0, k0:k1, i0:i1, j0:j1] + + # Write the subset locally + os.makedirs(opera_dir, exist_ok=True) + with h5py.File(out_file, 'w') as dst: + dst.create_dataset('latitude', data=lat_crop) + dst.create_dataset('longitude', data=lon_crop) + dst.create_dataset('height', data=height_crop) + dst.create_dataset('wet_delay', data=wet_crop[np.newaxis, ...]) + dst.create_dataset('hydrostatic_delay', data=hydro_crop[np.newaxis, ...]) + + return out_file + + +def dload_opera_files(missing_date_hour_list, opera_dir, geom_file=None, dem_range=None): """Download missing OPERA TROPO-ZENITH files via ASF Search. Uses per-day ASF queries constrained by OPERA collection IDs and date windows, then down-filters by required model-time tokens. + + When *geom_file* is provided (and ``fsspec`` is installed), files are + downloaded as spatial/vertical subsets using byte-range HTTP requests, + which significantly reduces bandwidth and disk usage. Falls back to + full-file download when subsetting is not possible. + + Parameters: missing_date_hour_list - list of (date_str, hour_str) tuples + opera_dir - str, local directory for OPERA files + geom_file - str or None, MintPy geometry file + dem_range - tuple(float, float) or None, + (min, max) DEM elevation in metres """ if len(missing_date_hour_list) == 0: return @@ -410,6 +697,15 @@ def dload_opera_files(missing_date_hour_list, opera_dir): print(f' WARNING: ASF authentication setup failed: {exc}') session = None + # Determine whether remote subsetting is possible + _can_subset = geom_file is not None + if _can_subset: + try: + import fsspec # noqa: F401 + except ImportError: + _can_subset = False + print(' NOTE: fsspec not installed; downloading full OPERA files.') + missing_tokens = sorted({_model_time_token(d, h) for d, h in missing_date_hour_list}) date_to_tokens = {} for token in missing_tokens: @@ -464,16 +760,33 @@ def dload_opera_files(missing_date_hour_list, opera_dir): print(f' {token}: {token_match_count[token]} -> no') continue - dload_results = asf.ASFSearchResults([token_best_product[token]]) - if session is not None: - dload_results.download(path=opera_dir, session=session) - else: - dload_results.download(path=opera_dir) - + best_product = token_best_product[token] + downloaded = False + + # Try remote subset download first + if _can_subset: + try: + out = dload_opera_file_subset( + best_product, opera_dir, geom_file, + dem_range=dem_range, session=session, + ) + if out is not None and _is_valid_opera_file(out): + downloaded = True + print(f' {token}: {token_match_count[token]} -> subset') + except Exception as exc: + print(f' {token}: subset download failed ({exc}), falling back to full download') + + # Fallback: download the full file via asf_search + if not downloaded: + dload_results = asf.ASFSearchResults([best_product]) + if session is not None: + dload_results.download(path=opera_dir, session=session) + else: + dload_results.download(path=opera_dir) + print(f' {token}: {token_match_count[token]} -> full') n_selected += 1 n_done += 1 - prog_bar.update(n_done, suffix=f'{token} (downloaded)') - print(f' {token}: {token_match_count[token]} -> yes') + prog_bar.update(n_done, suffix=f'{token} (done)') prog_bar.close() @@ -667,7 +980,20 @@ def run_tropo_opera(inps): inps.missing_opera_date_hour_list = missing_date_hour_list if len(missing_date_hour_list) > 0: - dload_opera_files(missing_date_hour_list, inps.opera_dir) + # Compute DEM range for remote subsetting + dem_range = None + try: + dem = readfile.read(inps.geom_file, datasetName='height')[0] + valid = dem[np.isfinite(dem)] + if valid.size > 0: + dem_range = (float(np.nanmin(valid)), float(np.nanmax(valid))) + except Exception: + pass + + dload_opera_files( + missing_date_hour_list, inps.opera_dir, + geom_file=inps.geom_file, dem_range=dem_range, + ) expected_patterns, matched_files, missing_date_hour_list = get_opera_file_status( opera_date_list=opera_date_list, opera_hour_list=opera_hour_list, From 99acda4bbad74420816faa9278227a0be6da3677 Mon Sep 17 00:00:00 2001 From: David Bekaert Date: Tue, 31 Mar 2026 12:18:50 +0200 Subject: [PATCH 06/24] fix: address review feedback (isort, memory scaling, docstring, CLI typos) - Remove extra blank line after imports to satisfy isort - Rewrite interpolation: cubic-lateral per height level then linear-vertical to reduce memory from O(npixels*ny*nx) to O(nz*npixels) - Fix docstring in read_opera_total_delay_cube (one-level buffer -> bracketing) - Fix CLI typos: Tropopsheric, trospheric, prerp, placeholder descriptions - Add _ensure_ascending helper for axis/data orientation --- src/mintpy/cli/tropo_opera.py | 11 +- src/mintpy/tropo_opera.py | 191 ++++++++++------------------------ 2 files changed, 63 insertions(+), 139 deletions(-) diff --git a/src/mintpy/cli/tropo_opera.py b/src/mintpy/cli/tropo_opera.py index be0d91468..1250109ff 100755 --- a/src/mintpy/cli/tropo_opera.py +++ b/src/mintpy/cli/tropo_opera.py @@ -13,12 +13,13 @@ ############################################################################ REFERENCE = """references: - in prerp to be added. + Bekaert, D. et al., OPERA L4 Tropospheric Zenith Delay products + derived from ECMWF HRES model (https://www.earthdata.nasa.gov/data/catalog/asf-opera-l4-tropo-zenith-v1-1). """ DIR_DEMO = """--dir ./OPERA - To be added - ... + Path to the directory for downloading and storing OPERA tropospheric delay + products. Files are automatically downloaded from ASF on demand. """ EXAMPLE = """example: @@ -28,7 +29,7 @@ def create_parser(subparsers=None): - synopsis = 'Tropospheric correction using OPERA computed Zenith Tropopsheric Delays from ECMWF HRES data (https://www.earthdata.nasa.gov/data/catalog/asf-opera-l4-tropo-zenith-v1-1)' + synopsis = 'Tropospheric correction using OPERA Zenith Tropospheric Delays from ECMWF HRES data (https://www.earthdata.nasa.gov/data/catalog/asf-opera-l4-tropo-zenith-v1-1)' epilog = REFERENCE + '\n' + DIR_DEMO + '\n' + EXAMPLE name = __name__.split('.')[-1] parser = create_argument_parser( @@ -40,7 +41,7 @@ def create_parser(subparsers=None): parser.add_argument('--dir','--opera-dir', dest='opera_dir', default='./OPERA', help='directory to store downloaded OPERA delays data (default: %(default)s).') parser.add_argument('-o', dest='cor_dis_file', - help='Output file name for trospheric corrected timeseries.') + help='Output file name for tropospheric corrected timeseries.') return parser diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index 17417c80c..68c16ef71 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -13,13 +13,11 @@ import h5py import numpy as np -from scipy.interpolate import RegularGridInterpolator import mintpy.cli.diff from mintpy.objects import timeseries from mintpy.utils import ptime, readfile, utils as ut, writefile - OPERA_MODEL_HOURS = (0, 6, 12, 18) OPERA_MIN_ACQ_DATE = '20160101' OPERA_TROPO_COLLECTIONS = [ @@ -138,8 +136,7 @@ def read_opera_total_delay_cube(opera_file, geom_file, dem_range=None, pad_cells geom_file - str, path to MintPy geometry HDF5 file dem_range - tuple(float, float) or None, (dem_min, dem_max) in metres. When provided the height dimension is also - cropped to the levels spanning [dem_min, dem_max] - plus a one-level buffer on each side. + cropped to the levels bracketing [dem_min, dem_max]. pad_cells - int, kept for backward compatibility (ignored) Returns a dict with keys: latitude, longitude, height, wet_delay, hydro_delay, total_delay, crop_info @@ -208,14 +205,30 @@ def get_geom_lat_lon_dem(geom_file): return lat2d, lon2d, dem +def _ensure_ascending(axis, data, dim): + """Reverse axis and corresponding data dimension if descending.""" + if axis.size > 1 and axis[1] < axis[0]: + axis = axis[::-1] + slices = [slice(None)] * data.ndim + slices[dim] = slice(None, None, -1) + data = data[tuple(slices)] + return axis, data + + def calc_zenith_delay_from_opera_file(opera_file, geom_file, pad_cells=3): """Calculate 2D zenith tropospheric delay map intersected with DEM. - Two-step interpolation: - 1. Linear interpolation in the vertical (height -> DEM elevation) - at each OPERA lat/lon grid node. - 2. Cubic interpolation in the lateral (lat/lon -> pixel locations). + Two-step interpolation (cubic lateral, linear vertical): + 1. At each OPERA height level, interpolate the 2D delay field from the + OPERA (lat, lon) grid to the pixel (lat, lon) locations using cubic + interpolation. This yields delay values at the pixel locations for + every height level — shape (nz, npixels). + 2. At each pixel, linearly interpolate along height to the DEM elevation. + + Memory usage is O(nz * npixels), where nz is typically 20–80 levels. """ + from scipy.interpolate import RegularGridInterpolator + lat2d, lon2d, dem = get_geom_lat_lon_dem(geom_file) # derive DEM range for vertical subsetting @@ -230,136 +243,44 @@ def calc_zenith_delay_from_opera_file(opera_file, geom_file, pad_cells=3): z_axis = np.asarray(cube['height'], dtype=np.float64) y_axis = np.asarray(cube['latitude'], dtype=np.float64) x_axis = np.asarray(cube['longitude'], dtype=np.float64) - data = np.asarray(cube['total_delay'], dtype=np.float64) - - # Ensure strictly ascending axes for interpolation. - if z_axis.size > 1 and z_axis[1] < z_axis[0]: - z_axis = z_axis[::-1] - data = data[::-1, :, :] - if y_axis.size > 1 and y_axis[1] < y_axis[0]: - y_axis = y_axis[::-1] - data = data[:, ::-1, :] - if x_axis.size > 1 and x_axis[1] < x_axis[0]: - x_axis = x_axis[::-1] - data = data[:, :, ::-1] - - # --- Step 1: Linear interpolation in the vertical ----------------------- - # At each OPERA (lat, lon) grid node, interpolate along height to the - # DEM elevation of every output pixel. This produces a 2-D slab per - # (lat_node, lon_node) → we collect them into a 2-D grid (ny, nx) for - # each output pixel, then do step 2. - # - # Efficient approach: for every output pixel, the DEM elevation is known. - # Build a 1-D linear interpolator per (lat, lon) column, evaluate at the - # pixel's DEM elevation → gives delay_2d[iy, ix] on the OPERA lat/lon grid. - # Then do cubic 2-D interpolation of that 2-D field to the pixel lat/lon. - - nrows, ncols = dem.shape + data = np.asarray(cube['total_delay'], dtype=np.float64) # (nz, ny, nx) + + # ensure strictly ascending axes + z_axis, data = _ensure_ascending(z_axis, data, 0) + y_axis, data = _ensure_ascending(y_axis, data, 1) + x_axis, data = _ensure_ascending(x_axis, data, 2) + + nz = len(z_axis) dem_flat = dem.ravel().astype(np.float64) lat_flat = lat2d.ravel().astype(np.float64) lon_flat = lon2d.ravel().astype(np.float64) + npixels = len(dem_flat) - ny, nx = len(y_axis), len(x_axis) - - # For each OPERA (lat, lon) column, build a 1-D linear interpolator - # along height, then evaluate at *all* output pixel DEM values at once. - # Result: delay_at_dem[j, i, pixel] but that's memory-heavy. - # Instead, use RegularGridInterpolator in 1-D (height) per column → too slow. - # - # Better: use scipy interp1d broadcast or manual linear interp. - # Use np.interp per column vectorised over the output pixels. - # - # Most efficient: do the vertical interp as a single - # RegularGridInterpolator(method='linear') in 3-D, then re-interpolate - # the result in 2-D with cubic. But that couples axes again. - # - # Cleanest two-step: collapse vertical first on the OPERA grid, then - # interpolate laterally. - # For each output pixel p with DEM elevation h_p: - # delay_on_opera_grid[iy, ix] = interp1d(z_axis, data[:, iy, ix])(h_p) - # Then: ztd[p] = interp2d(y_axis, x_axis, delay_on_opera_grid)(lat_p, lon_p) - # - # This is O(npixels * ny * nx) if done naively. - # Vectorise: for all pixels at once, linear interp along z for all columns. - - # Pre-compute vertical interpolation weights for each pixel - # np.interp works on 1-D; we need to interp data[:, iy, ix] at dem values. - # Reshape data to (nz, ny*nx), interp each column at all dem values → large. - # Instead: interp all columns at each unique dem value. - - # Practical approach: loop over OPERA grid columns (ny*nx is small, ~hundreds) - # and use np.interp for all pixels at once per column. - # Result: ztd_on_grid shape (ny, nx, npixels) → then for each pixel, do 2-D interp. - # That's still large. Better: for each pixel, do 1-D vertical interp at - # 1 height value across all columns, giving (ny, nx), then 2-D interp. - # But looping over pixels is slow. - - # Most practical vectorised approach: - # Step 1: vertical linear interp → 2D field on OPERA grid for each pixel - # Use broadcasting: find bracketing indices and weights once per pixel. + # Step 1: cubic lateral interpolation at each height level + # Evaluate the 2D (lat, lon) field at pixel locations → (nz, npixels) + pts_2d = np.column_stack([lat_flat, lon_flat]) + delay_at_pixels = np.empty((nz, npixels), dtype=np.float64) + for k in range(nz): + interp_2d = RegularGridInterpolator( + (y_axis, x_axis), data[k, :, :], + method='cubic', bounds_error=False, fill_value=np.nan, + ) + delay_at_pixels[k, :] = interp_2d(pts_2d) + + # Step 2: linear vertical interpolation at each pixel + # Find bracketing height indices and fractional weights idx = np.searchsorted(z_axis, dem_flat, side='right') - 1 - idx = np.clip(idx, 0, len(z_axis) - 2) - # fractional weight + idx = np.clip(idx, 0, nz - 2) dz = z_axis[idx + 1] - z_axis[idx] - dz[dz == 0] = 1.0 # avoid division by zero - w = (dem_flat - z_axis[idx]) / dz - w = np.clip(w, 0.0, 1.0) - - # data shape: (nz, ny, nx) - # For each pixel p: delay_2d[:, :, p] = data[idx[p], :, :] * (1-w[p]) + data[idx[p]+1, :, :] * w[p] - # Vectorise: gather the two bracketing slabs for all pixels - # data_lo[p, iy, ix] = data[idx[p], iy, ix] → shape (npixels, ny, nx) - npixels = len(dem_flat) - data_lo = data[idx] # shape (npixels, ny, nx) - data_hi = data[idx + 1] # shape (npixels, ny, nx) - # Weighted blend: shape (npixels, ny, nx) - w3 = w[:, np.newaxis, np.newaxis] - delay_2d = data_lo * (1.0 - w3) + data_hi * w3 # (npixels, ny, nx) - - # --- Step 2: Cubic interpolation in lat/lon per pixel ------------------- - # For each pixel, we have a 2-D delay field on the OPERA (y_axis, x_axis) - # grid. Interpolate to the pixel's (lat, lon) using cubic. - # - # Batch approach: all pixels share the same OPERA grid, and each pixel has - # its own 2-D field. Use RegularGridInterpolator per pixel → too slow. - # - # Faster: the lateral interpolation weights depend only on (lat, lon) and - # are the same for all pixels at the same location. Use - # RegularGridInterpolator once with a dummy z-axis to vectorise: - # Build a 3-D field (npixels, ny, nx) with "pixel index" as the first - # axis → but that's not a regular grid in the first axis. - # - # Best practical approach: pre-compute cubic interpolation weights for - # each pixel's (lat, lon) on the OPERA grid, then apply them. - # scipy.interpolate doesn't expose weights directly, so use map_coordinates - # or manual cubic. - # - # Simple and efficient: use scipy.ndimage.map_coordinates with order=3 - # (cubic) on each pixel's 2-D field. But looping over pixels is slow. - # - # Actually, since all pixels share the OPERA grid, we can convert - # (lat, lon) to fractional grid indices once, then use them for all pixels. - - from scipy.ndimage import map_coordinates - - # Convert pixel lat/lon to fractional indices in the OPERA grid - # y_axis is ascending, x_axis is ascending - iy_frac = np.interp(lat_flat, y_axis, np.arange(ny)) - ix_frac = np.interp(lon_flat, x_axis, np.arange(nx)) - - # map_coordinates on delay_2d: shape (npixels, ny, nx) - # We want to evaluate delay_2d[p, iy_frac[p], ix_frac[p]] for each p. - # Coordinates: axis0=p (exact integer), axis1=iy_frac, axis2=ix_frac - coords = np.array([ - np.arange(npixels, dtype=np.float64), - iy_frac, - ix_frac, - ]) - ztd_flat = map_coordinates(delay_2d, coords, order=3, mode='nearest') + dz[dz == 0] = 1.0 + w = np.clip((dem_flat - z_axis[idx]) / dz, 0.0, 1.0) - ztd = ztd_flat.reshape(dem.shape).astype(np.float32) + # linear blend between bracketing levels + val_lo = delay_at_pixels[idx, np.arange(npixels)] + val_hi = delay_at_pixels[idx + 1, np.arange(npixels)] + ztd_flat = val_lo * (1.0 - w) + val_hi * w - # mask invalid pixels from DEM (NaN/Inf); preserve valid zero-elevation pixels + ztd = ztd_flat.reshape(dem.shape).astype(np.float32) ztd[~np.isfinite(dem)] = np.nan return ztd, cube @@ -701,9 +622,11 @@ def dload_opera_files(missing_date_hour_list, opera_dir, geom_file=None, dem_ran _can_subset = geom_file is not None if _can_subset: try: - import fsspec # noqa: F401 - except ImportError: + from importlib.util import find_spec + _can_subset = find_spec('fsspec') is not None + except Exception: _can_subset = False + if not _can_subset: print(' NOTE: fsspec not installed; downloading full OPERA files.') missing_tokens = sorted({_model_time_token(d, h) for d, h in missing_date_hour_list}) @@ -987,8 +910,8 @@ def run_tropo_opera(inps): valid = dem[np.isfinite(dem)] if valid.size > 0: dem_range = (float(np.nanmin(valid)), float(np.nanmax(valid))) - except Exception: - pass + except Exception as exc: + print(f' WARNING: could not read DEM for subsetting: {exc}') dload_opera_files( missing_date_hour_list, inps.opera_dir, From cfef0171044ee64442e72c4af6f39bb8e324fca1 Mon Sep 17 00:00:00 2001 From: David Bekaert Date: Tue, 31 Mar 2026 12:25:44 +0200 Subject: [PATCH 07/24] fix: pyupgrade f-string parens + add asf_search to requirements - Remove redundant parentheses around f-strings (pyupgrade --py36-plus) - Add asf_search to requirements.txt (needed for OPERA download) --- requirements.txt | 1 + src/mintpy/tropo_opera.py | 8 ++++---- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/requirements.txt b/requirements.txt index cce66eee5..ee082d8f3 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,6 @@ pip argcomplete +asf_search cartopy cvxopt dask>=1.0 diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index 68c16ef71..fe41b0dfd 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -800,8 +800,8 @@ def calculate_delay_timeseries(tropo_file, dis_file, geom_file, opera_dir): supported_dates, unsupported_dates = _split_supported_dates(date_list, min_date=OPERA_MIN_ACQ_DATE) if len(unsupported_dates) > 0: - print((f'WARNING: {len(unsupported_dates)} acquisition date(s) are earlier than ' - f'OPERA availability ({OPERA_MIN_ACQ_DATE}) and will be skipped.')) + print(f'WARNING: {len(unsupported_dates)} acquisition date(s) are earlier than ' + f'OPERA availability ({OPERA_MIN_ACQ_DATE}) and will be skipped.') for d in unsupported_dates: print(f' {d}') @@ -880,8 +880,8 @@ def run_tropo_opera(inps): date_list, utc_sec = read_inps2date_time(inps) supported_dates, unsupported_dates = _split_supported_dates(date_list, min_date=OPERA_MIN_ACQ_DATE) if len(unsupported_dates) > 0: - print((f'WARNING: {len(unsupported_dates)} acquisition date(s) are earlier than ' - f'OPERA availability ({OPERA_MIN_ACQ_DATE}) and will be skipped.')) + print(f'WARNING: {len(unsupported_dates)} acquisition date(s) are earlier than ' + f'OPERA availability ({OPERA_MIN_ACQ_DATE}) and will be skipped.') for d in unsupported_dates: print(f' {d}') From 9a3cbc4bee24c2fe74d2d2769e554afc2adf5e6e Mon Sep 17 00:00:00 2001 From: David Bekaert Date: Tue, 31 Mar 2026 12:28:36 +0200 Subject: [PATCH 08/24] add fsspec[http] to requirements for OPERA remote subsetting fsspec[http] includes aiohttp, enabling byte-range HTTP downloads for efficient spatial/vertical subsetting of OPERA ZTD files. --- requirements.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index ee082d8f3..c3243be52 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,10 +1,11 @@ pip argcomplete -asf_search +asf_search>=12.0.0 cartopy cvxopt dask>=1.0 dask-jobqueue>=0.3 +fsspec[http] h5py joblib lxml From 592494e423d347749796e5c2102c4c17ad3edc21 Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Wed, 1 Apr 2026 22:41:12 +0800 Subject: [PATCH 09/24] fsspec requires python=3.10+ --- Dockerfile | 2 +- docs/README.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Dockerfile b/Dockerfile index a8b9f983f..8e59b6d1d 100644 --- a/Dockerfile +++ b/Dockerfile @@ -33,7 +33,7 @@ ENV PATH=/opt/conda/bin:${PATH} ARG MINTPY_HOME=/home/mambauser/tools/MintPy COPY --chown=mambauser:mambauser . ${MINTPY_HOME}/ -ARG PYTHON_VERSION="3.9" +ARG PYTHON_VERSION="3.11" RUN micromamba install -y -n base -c conda-forge python=${PYTHON_VERSION} \ jupyterlab ipympl gdal">=3" isce2 -f ${MINTPY_HOME}/requirements.txt && \ python -m pip install --no-cache-dir ${MINTPY_HOME} && \ diff --git a/docs/README.md b/docs/README.md index d87251cce..ad7e8ce0b 100644 --- a/docs/README.md +++ b/docs/README.md @@ -1,4 +1,4 @@ -[![Language](https://img.shields.io/badge/python-3.8%2B-blue.svg?style=flat-square)](https://www.python.org/) +[![Language](https://img.shields.io/badge/python-3.10%2B-blue.svg?style=flat-square)](https://www.python.org/) [![Docs Status](https://readthedocs.org/projects/mintpy/badge/?version=latest&style=flat-square)](https://mintpy.readthedocs.io/?badge=latest) [![CircleCI](https://img.shields.io/circleci/build/github/insarlab/MintPy.svg?logo=circleci&label=tests&style=flat-square)](https://circleci.com/gh/insarlab/MintPy) [![Docker Status](https://img.shields.io/github/actions/workflow/status/insarlab/MintPy/build-docker.yml?label=docker&style=flat-square&logo=docker&logoColor=white)](https://github.com/insarlab/MintPy/pkgs/container/mintpy) From 49ea84272142779baa2ae1ccdf9382db2479675b Mon Sep 17 00:00:00 2001 From: Simran S Sangha Date: Wed, 1 Apr 2026 17:08:33 -0700 Subject: [PATCH 10/24] Implement geometry grid reprojection support for UTM TS --- src/mintpy/tropo_opera.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index fe41b0dfd..95bdaac7a 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -228,9 +228,27 @@ def calc_zenith_delay_from_opera_file(opera_file, geom_file, pad_cells=3): Memory usage is O(nz * npixels), where nz is typically 20–80 levels. """ from scipy.interpolate import RegularGridInterpolator + from pyproj import Transformer + from mintpy.utils import readfile lat2d, lon2d, dem = get_geom_lat_lon_dem(geom_file) + # Reprojection block + atr = readfile.read_attribute(geom_file) + epsg = atr.get('EPSG', '4326') + + if str(epsg) != '4326': + print(f"Reprojecting geometry grids from EPSG:{epsg} to EPSG:4326 (degrees)...") + # always_xy=True ensures the order is (X, Y) -> (Longitude, Latitude) + transformer = Transformer.from_crs(f"EPSG:{epsg}", "EPSG:4326", always_xy=True) + + # lon2d is X (Easting), lat2d is Y (Northing) + lon2d_deg, lat2d_deg = transformer.transform(lon2d, lat2d) + + # Replace the meter arrays with the degree arrays for the interpolator + lat2d = lat2d_deg + lon2d = lon2d_deg + # derive DEM range for vertical subsetting valid_dem = dem[np.isfinite(dem)] if valid_dem.size > 0: From e6feefb60caeb0b9f42d1dea88633333c37cbca5 Mon Sep 17 00:00:00 2001 From: David Bekaert Date: Thu, 2 Apr 2026 08:57:44 +0200 Subject: [PATCH 11/24] fix: remove duplicate readfile import and trailing whitespace --- src/mintpy/tropo_opera.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index 95bdaac7a..21dfd2fdf 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -229,22 +229,21 @@ def calc_zenith_delay_from_opera_file(opera_file, geom_file, pad_cells=3): """ from scipy.interpolate import RegularGridInterpolator from pyproj import Transformer - from mintpy.utils import readfile lat2d, lon2d, dem = get_geom_lat_lon_dem(geom_file) # Reprojection block atr = readfile.read_attribute(geom_file) epsg = atr.get('EPSG', '4326') - + if str(epsg) != '4326': print(f"Reprojecting geometry grids from EPSG:{epsg} to EPSG:4326 (degrees)...") # always_xy=True ensures the order is (X, Y) -> (Longitude, Latitude) transformer = Transformer.from_crs(f"EPSG:{epsg}", "EPSG:4326", always_xy=True) - + # lon2d is X (Easting), lat2d is Y (Northing) lon2d_deg, lat2d_deg = transformer.transform(lon2d, lat2d) - + # Replace the meter arrays with the degree arrays for the interpolator lat2d = lat2d_deg lon2d = lon2d_deg From d069e644aa8d7eddbd8cd80d653a60ad71881d58 Mon Sep 17 00:00:00 2001 From: David Bekaert Date: Thu, 2 Apr 2026 09:39:48 +0200 Subject: [PATCH 12/24] fix: sort imports (isort) in calc_zenith_delay_from_opera_file --- src/mintpy/tropo_opera.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index 21dfd2fdf..758417188 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -227,8 +227,8 @@ def calc_zenith_delay_from_opera_file(opera_file, geom_file, pad_cells=3): Memory usage is O(nz * npixels), where nz is typically 20–80 levels. """ - from scipy.interpolate import RegularGridInterpolator from pyproj import Transformer + from scipy.interpolate import RegularGridInterpolator lat2d, lon2d, dem = get_geom_lat_lon_dem(geom_file) From e3da8056ac7d0aed6826ce3cd24f0bdd9d65d133 Mon Sep 17 00:00:00 2001 From: Simran S Sangha Date: Sun, 5 Apr 2026 12:54:35 -0700 Subject: [PATCH 13/24] Refactor calc_zenith_delay_from_opera_file parameters --- src/mintpy/tropo_opera.py | 41 +++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 21 deletions(-) diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index 758417188..f3d25bbc1 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -215,7 +215,9 @@ def _ensure_ascending(axis, data, dim): return axis, data -def calc_zenith_delay_from_opera_file(opera_file, geom_file, pad_cells=3): +def calc_zenith_delay_from_opera_file( + opera_file, geom_file, lat2d, lon2d, dem, pad_cells=3 +): """Calculate 2D zenith tropospheric delay map intersected with DEM. Two-step interpolation (cubic lateral, linear vertical): @@ -227,27 +229,8 @@ def calc_zenith_delay_from_opera_file(opera_file, geom_file, pad_cells=3): Memory usage is O(nz * npixels), where nz is typically 20–80 levels. """ - from pyproj import Transformer from scipy.interpolate import RegularGridInterpolator - lat2d, lon2d, dem = get_geom_lat_lon_dem(geom_file) - - # Reprojection block - atr = readfile.read_attribute(geom_file) - epsg = atr.get('EPSG', '4326') - - if str(epsg) != '4326': - print(f"Reprojecting geometry grids from EPSG:{epsg} to EPSG:4326 (degrees)...") - # always_xy=True ensures the order is (X, Y) -> (Longitude, Latitude) - transformer = Transformer.from_crs(f"EPSG:{epsg}", "EPSG:4326", always_xy=True) - - # lon2d is X (Easting), lat2d is Y (Northing) - lon2d_deg, lat2d_deg = transformer.transform(lon2d, lat2d) - - # Replace the meter arrays with the degree arrays for the interpolator - lat2d = lat2d_deg - lon2d = lon2d_deg - # derive DEM range for vertical subsetting valid_dem = dem[np.isfinite(dem)] if valid_dem.size > 0: @@ -800,6 +783,8 @@ def get_dataset_size(fname): ############################################################################ def calculate_delay_timeseries(tropo_file, dis_file, geom_file, opera_dir): """Calculate OPERA delay time-series and write to HDF5 file.""" + from pyproj import Transformer + atr = readfile.read_attribute(dis_file) ftype = atr['FILE_TYPE'] @@ -864,9 +849,23 @@ def calculate_delay_timeseries(tropo_file, dis_file, geom_file, opera_dir): print(f'read incidenceAngle from file: {geom_file}') inc_angle = readfile.read(geom_file, datasetName='incidenceAngle')[0].astype(np.float32) + # Read geometry once + print(f'Reading and preparing geometry grids from {geom_file}...') + lat2d, lon2d, dem = get_geom_lat_lon_dem(geom_file) + + atr_geom = readfile.read_attribute(geom_file) + epsg = atr_geom.get('EPSG', '4326') + + if str(epsg) != '4326': + print(f"Reprojecting geometry grids from EPSG:{epsg} to EPSG:4326 (degrees)...") + transformer = Transformer.from_crs(f"EPSG:{epsg}", "EPSG:4326", always_xy=True) + lon2d, lat2d = transformer.transform(lon2d, lat2d) + prog_bar = ptime.progressBar(maxValue=num_date) for i in range(num_date): - ztd, _ = calc_zenith_delay_from_opera_file(opera_files[i], geom_file, pad_cells=3) + ztd, _ = calc_zenith_delay_from_opera_file( + opera_files[i], geom_file, lat2d, lon2d, dem, pad_cells=3 + ) delay_los = _project_zenith_to_los(ztd, inc_angle) block = [i, i + 1, 0, length, 0, width] From b12bc2ca3a0ca9040be3edf7dfd25cb6b36d7940 Mon Sep 17 00:00:00 2001 From: Simran S Sangha Date: Mon, 6 Apr 2026 01:47:57 -0700 Subject: [PATCH 14/24] Fix formatting in tropo_opera.py --- src/mintpy/tropo_opera.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index f3d25bbc1..374c6d794 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -852,7 +852,7 @@ def calculate_delay_timeseries(tropo_file, dis_file, geom_file, opera_dir): # Read geometry once print(f'Reading and preparing geometry grids from {geom_file}...') lat2d, lon2d, dem = get_geom_lat_lon_dem(geom_file) - + atr_geom = readfile.read_attribute(geom_file) epsg = atr_geom.get('EPSG', '4326') From bba8153f3727ae484e0a75ec6defe8d1b4242c79 Mon Sep 17 00:00:00 2001 From: Simran S Sangha Date: Mon, 6 Apr 2026 11:32:55 -0700 Subject: [PATCH 15/24] Refactor code for improved readability and structure --- src/mintpy/tropo_opera.py | 50 ++++++++++++++++++++++++++++----------- 1 file changed, 36 insertions(+), 14 deletions(-) diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index 374c6d794..1fabb5de9 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -151,7 +151,9 @@ def read_opera_total_delay_cube(opera_file, geom_file, dem_range=None, pad_cells lon = fobj['longitude'][:] height = fobj['height'][:] - (i0, i1, j0, j1), crop_info = get_opera_crop_indices(lat, lon, geom_file, pad_cells=pad_cells) + (i0, i1, j0, j1), crop_info = get_opera_crop_indices( + lat, lon, geom_file, pad_cells=pad_cells + ) # vertical subsetting if dem_range is not None: @@ -229,7 +231,7 @@ def calc_zenith_delay_from_opera_file( Memory usage is O(nz * npixels), where nz is typically 20–80 levels. """ - from scipy.interpolate import RegularGridInterpolator + from scipy.ndimage import map_coordinates # derive DEM range for vertical subsetting valid_dem = dem[np.isfinite(dem)] @@ -238,7 +240,9 @@ def calc_zenith_delay_from_opera_file( else: dem_range = None - cube = read_opera_total_delay_cube(opera_file, geom_file, dem_range=dem_range, pad_cells=pad_cells) + cube = read_opera_total_delay_cube(opera_file, geom_file, + dem_range=dem_range, + pad_cells=pad_cells) z_axis = np.asarray(cube['height'], dtype=np.float64) y_axis = np.asarray(cube['latitude'], dtype=np.float64) @@ -256,16 +260,25 @@ def calc_zenith_delay_from_opera_file( lon_flat = lon2d.ravel().astype(np.float64) npixels = len(dem_flat) - # Step 1: cubic lateral interpolation at each height level - # Evaluate the 2D (lat, lon) field at pixel locations → (nz, npixels) - pts_2d = np.column_stack([lat_flat, lon_flat]) + # Calculate the fixed grid spacing of the OPERA cube + dy = y_axis[1] - y_axis[0] + dx = x_axis[1] - x_axis[0] + + # Convert geographical lat/lon into fractional array indices for the C-engine + y_idx = (lat_flat - y_axis[0]) / dy + x_idx = (lon_flat - x_axis[0]) / dx + coords = np.vstack([y_idx, x_idx]) + + # Step 1: cubic lateral interpolation at each height level using map_coordinates delay_at_pixels = np.empty((nz, npixels), dtype=np.float64) for k in range(nz): - interp_2d = RegularGridInterpolator( - (y_axis, x_axis), data[k, :, :], - method='cubic', bounds_error=False, fill_value=np.nan, + delay_at_pixels[k, :] = map_coordinates( + data[k, :, :], + coords, + order=3, # order=3 means cubic interpolation + mode='constant', # If a pixel is outside the padded bounds... + cval=np.nan # ...fill it with NaN ) - delay_at_pixels[k, :] = interp_2d(pts_2d) # Step 2: linear vertical interpolation at each pixel # Find bracketing height indices and fractional weights @@ -800,7 +813,9 @@ def calculate_delay_timeseries(tropo_file, dis_file, geom_file, opera_dir): raise ValueError(f'CENTER_LINE_UTC is missing in metadata of file: {dis_file}') utc_sec = float(atr['CENTER_LINE_UTC']) - supported_dates, unsupported_dates = _split_supported_dates(date_list, min_date=OPERA_MIN_ACQ_DATE) + supported_dates, unsupported_dates = _split_supported_dates( + date_list, min_date=OPERA_MIN_ACQ_DATE + ) if len(unsupported_dates) > 0: print(f'WARNING: {len(unsupported_dates)} acquisition date(s) are earlier than ' f'OPERA availability ({OPERA_MIN_ACQ_DATE}) and will be skipped.') @@ -819,13 +834,18 @@ def calculate_delay_timeseries(tropo_file, dis_file, geom_file, opera_dir): for acq_date, model_date, model_hour in zip(supported_dates, model_dates, model_hours): fname = _get_best_local_opera_file(model_date, model_hour, opera_dir) if fname is None: - print(f'WARNING: NO OPERA file found for acquisition {acq_date} -> {model_date} {model_hour}:00') + print(f'WARNING: NO OPERA file found ' + f'for acquisition {acq_date} -> {model_date} ' + f'{model_hour}:00') continue used_dates.append(acq_date) opera_files.append(fname) if len(opera_files) == 0: - raise RuntimeError('No OPERA files found for requested acquisitions. Cannot create delay time-series.') + raise RuntimeError( + 'No OPERA files found for requested ' + 'acquisitions. Cannot create delay time-series.' + ) if _run_or_skip_opera(opera_files, used_dates, tropo_file, geom_file) == 'skip': return tropo_file @@ -894,7 +914,9 @@ def run_tropo_opera(inps): print(f'output corrected file : {inps.cor_dis_file}') date_list, utc_sec = read_inps2date_time(inps) - supported_dates, unsupported_dates = _split_supported_dates(date_list, min_date=OPERA_MIN_ACQ_DATE) + supported_dates, unsupported_dates = _split_supported_dates( + date_list, min_date=OPERA_MIN_ACQ_DATE + ) if len(unsupported_dates) > 0: print(f'WARNING: {len(unsupported_dates)} acquisition date(s) are earlier than ' f'OPERA availability ({OPERA_MIN_ACQ_DATE}) and will be skipped.') From 5faf8d4e1ad209929a443da669f6c4c76aced52e Mon Sep 17 00:00:00 2001 From: Simran S Sangha Date: Mon, 13 Apr 2026 14:37:48 -0700 Subject: [PATCH 16/24] Optimize calc_zenith_delay_from_opera_file function Refactor zenith delay calculation for memory efficiency and clarity. --- src/mintpy/tropo_opera.py | 65 +++++++++++++++++++++------------------ 1 file changed, 35 insertions(+), 30 deletions(-) diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index 1fabb5de9..2813d89a2 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -220,17 +220,7 @@ def _ensure_ascending(axis, data, dim): def calc_zenith_delay_from_opera_file( opera_file, geom_file, lat2d, lon2d, dem, pad_cells=3 ): - """Calculate 2D zenith tropospheric delay map intersected with DEM. - - Two-step interpolation (cubic lateral, linear vertical): - 1. At each OPERA height level, interpolate the 2D delay field from the - OPERA (lat, lon) grid to the pixel (lat, lon) locations using cubic - interpolation. This yields delay values at the pixel locations for - every height level — shape (nz, npixels). - 2. At each pixel, linearly interpolate along height to the DEM elevation. - - Memory usage is O(nz * npixels), where nz is typically 20–80 levels. - """ + """Calculate 2D zenith tropospheric delay map intersected with DEM.""" from scipy.ndimage import map_coordinates # derive DEM range for vertical subsetting @@ -267,33 +257,48 @@ def calc_zenith_delay_from_opera_file( # Convert geographical lat/lon into fractional array indices for the C-engine y_idx = (lat_flat - y_axis[0]) / dy x_idx = (lon_flat - x_axis[0]) / dx - coords = np.vstack([y_idx, x_idx]) - - # Step 1: cubic lateral interpolation at each height level using map_coordinates - delay_at_pixels = np.empty((nz, npixels), dtype=np.float64) - for k in range(nz): - delay_at_pixels[k, :] = map_coordinates( - data[k, :, :], - coords, - order=3, # order=3 means cubic interpolation - mode='constant', # If a pixel is outside the padded bounds... - cval=np.nan # ...fill it with NaN - ) - # Step 2: linear vertical interpolation at each pixel - # Find bracketing height indices and fractional weights + # Memory efficient logic + # Find bracketing height indices and fractional weights FIRST idx = np.searchsorted(z_axis, dem_flat, side='right') - 1 idx = np.clip(idx, 0, nz - 2) dz = z_axis[idx + 1] - z_axis[idx] dz[dz == 0] = 1.0 w = np.clip((dem_flat - z_axis[idx]) / dz, 0.0, 1.0) - # linear blend between bracketing levels - val_lo = delay_at_pixels[idx, np.arange(npixels)] - val_hi = delay_at_pixels[idx + 1, np.arange(npixels)] - ztd_flat = val_lo * (1.0 - w) + val_hi * w + # Initialize final output array (only takes ~260 MB of RAM instead of 15.6 GB) + ztd_flat = np.zeros(npixels, dtype=np.float32) + + for k in range(nz): + # Find ONLY the pixels that use layer 'k' as their lower or upper bracket + mask_lo = (idx == k) + mask_hi = (idx + 1 == k) + mask_any = mask_lo | mask_hi + + # If no pixels in the entire image live at this altitude, skip the layer entirely! + if not np.any(mask_any): + continue + + # Interpolate ONLY for the required pixels + coords_k = np.vstack([y_idx[mask_any], x_idx[mask_any]]) + delay_k = map_coordinates( + data[k, :, :], + coords_k, + order=3, + mode='constant', + cval=np.nan + ) + + # Apply the weights directly to the final ZTD array + if np.any(mask_lo): + valid_lo = mask_lo[mask_any] + ztd_flat[mask_lo] += delay_k[valid_lo] * (1.0 - w[mask_lo]) + + if np.any(mask_hi): + valid_hi = mask_hi[mask_any] + ztd_flat[mask_hi] += delay_k[valid_hi] * w[mask_hi] - ztd = ztd_flat.reshape(dem.shape).astype(np.float32) + ztd = ztd_flat.reshape(dem.shape) ztd[~np.isfinite(dem)] = np.nan return ztd, cube From c3c0a71402d19ce80d383634511c7d6ef6a6566c Mon Sep 17 00:00:00 2001 From: Simran S Sangha Date: Mon, 13 Apr 2026 15:48:50 -0700 Subject: [PATCH 17/24] Update authorship and improve docstring clarity --- src/mintpy/tropo_opera.py | 35 +++++++++++++++++++++++------------ 1 file changed, 23 insertions(+), 12 deletions(-) diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index 2813d89a2..add346d64 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -1,7 +1,7 @@ ############################################################ # Program is part of MintPy # # Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi # -# Author: David Bekaert, March 2026 # +# Author: David Bekaert, Simran Sangha, March 2026 # ############################################################ @@ -51,7 +51,8 @@ def get_geom_lat_lon_bounds(geom_file): def get_opera_crop_indices(lat, lon, geom_file, pad_cells=3): - """Get OPERA cube crop indices from geometry bounds with fixed ±0.2 deg padding. + """ + Get OPERA cube crop indices from geometry bounds with fixed ±0.2 deg padding. The pad_cells argument is kept for backward compatibility and ignored. """ @@ -89,7 +90,8 @@ def get_opera_crop_indices(lat, lon, geom_file, pad_cells=3): def get_opera_height_crop_indices(height_levels, dem_min, dem_max): - """Get height level indices spanning the DEM range. + """ + Get height level indices spanning the DEM range. The vertical interpolation is linear, so only the levels that bracket [dem_min, dem_max] are needed — no extra buffer beyond that. @@ -130,7 +132,8 @@ def get_opera_height_crop_indices(height_levels, dem_min, dem_max): def read_opera_total_delay_cube(opera_file, geom_file, dem_range=None, pad_cells=3): - """Read and crop OPERA delay cube, returning total zenith delay. + """ + Read and crop OPERA delay cube, returning total zenith delay. Parameters: opera_file - str, path to OPERA netCDF file geom_file - str, path to MintPy geometry HDF5 file @@ -181,7 +184,8 @@ def read_opera_total_delay_cube(opera_file, geom_file, dem_range=None, pad_cells def get_geom_lat_lon_dem(geom_file): - """Read DEM and latitude/longitude grids from geometry file. + """ + Read DEM and latitude/longitude grids from geometry file. Returns: lat2d/lon2d/dem in 2D np.ndarray with same shape """ @@ -306,7 +310,8 @@ def calc_zenith_delay_from_opera_file( ############################################################################ def read_inps2date_time(inps): - """Read acquisition date/time info from input arguments. + """ + Read acquisition date/time info from input arguments. Parameters: inps - Namespace for input arguments Returns: date_list - list(str), acquisition dates in YYYYMMDD @@ -338,7 +343,8 @@ def read_inps2date_time(inps): def nearest_opera_product_time(date_str, utc_sec): - """Round one acquisition date/time to the nearest OPERA model datetime. + """ + Round one acquisition date/time to the nearest OPERA model datetime. Parameters: date_str - str, acquisition date in YYYYMMDD utc_sec - float, acquisition UTC time in seconds @@ -361,7 +367,8 @@ def nearest_opera_product_time(date_str, utc_sec): def get_opera_date_time_list(date_list, utc_sec): - """Build required OPERA date/hour list for all acquisitions. + """ + Build required OPERA date/hour list for all acquisitions. Parameters: date_list - list(str), acquisition dates in YYYYMMDD utc_sec - float, acquisition UTC time in seconds @@ -379,7 +386,8 @@ def get_opera_date_time_list(date_list, utc_sec): def get_expected_opera_file_patterns(opera_date_list, opera_hour_list): - """Create expected OPERA ZTD filename patterns for model times. + """ + Create expected OPERA ZTD filename patterns for model times. The production timestamp is left as wildcard. @@ -394,7 +402,8 @@ def get_expected_opera_file_patterns(opera_date_list, opera_hour_list): def get_opera_file_status(opera_date_list, opera_hour_list, opera_dir): - """Get matched OPERA files and missing model date/hour list. + """ + Get matched OPERA files and missing model date/hour list. Returns: expected_patterns - list(str) matched_files - dict(str, list(str)) @@ -523,7 +532,8 @@ def _get_product_filename(product): def dload_opera_file_subset(product, opera_dir, geom_file, dem_range=None, session=None): - """Download a spatially and vertically subsetted OPERA file. + """ + Download a spatially and vertically subsetted OPERA file. Uses ``fsspec`` + ``h5py`` to open the remote file via byte-range HTTP requests, reads only the required spatial/vertical slice, and writes the @@ -595,7 +605,8 @@ def dload_opera_file_subset(product, opera_dir, geom_file, dem_range=None, sessi def dload_opera_files(missing_date_hour_list, opera_dir, geom_file=None, dem_range=None): - """Download missing OPERA TROPO-ZENITH files via ASF Search. + """ + Download missing OPERA TROPO-ZENITH files via ASF Search. Uses per-day ASF queries constrained by OPERA collection IDs and date windows, then down-filters by required model-time tokens. From 7ed17a343b28ce0ab6ee680a28892ab589a93d3c Mon Sep 17 00:00:00 2001 From: Simran S Sangha Date: Mon, 13 Apr 2026 15:49:21 -0700 Subject: [PATCH 18/24] Fix comment formatting in __main__.py --- src/mintpy/__main__.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/mintpy/__main__.py b/src/mintpy/__main__.py index ded28f330..0492d52bc 100644 --- a/src/mintpy/__main__.py +++ b/src/mintpy/__main__.py @@ -5,7 +5,8 @@ ############################################################ -"""Command line interface for MintPy. +""" +Command line interface for MintPy. The Miami INsar Time-series software in PYthon (MintPy as /mint pai/) is an open-source package for Interferometric Synthetic Aperture Radar From fa8b5b724201bf83646afd77bd2560564ceb3ff5 Mon Sep 17 00:00:00 2001 From: Simran S Sangha Date: Mon, 13 Apr 2026 16:02:15 -0700 Subject: [PATCH 19/24] Update author information and fix docstring format --- src/mintpy/tropo_opera.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index add346d64..c74a657d4 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -1,7 +1,7 @@ ############################################################ # Program is part of MintPy # # Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi # -# Author: David Bekaert, Simran Sangha, March 2026 # +# Author: David Bekaert, March 2026 # ############################################################ @@ -51,8 +51,7 @@ def get_geom_lat_lon_bounds(geom_file): def get_opera_crop_indices(lat, lon, geom_file, pad_cells=3): - """ - Get OPERA cube crop indices from geometry bounds with fixed ±0.2 deg padding. + """Get OPERA cube crop indices from geometry bounds with fixed ±0.2 deg padding. The pad_cells argument is kept for backward compatibility and ignored. """ From 26ac2376be6febbbaf09a5ee4ca6acbc71236eda Mon Sep 17 00:00:00 2001 From: Simran S Sangha Date: Mon, 13 Apr 2026 16:02:51 -0700 Subject: [PATCH 20/24] Fix docstring formatting in __main__.py Removed unnecessary quotes from the docstring. --- src/mintpy/__main__.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/mintpy/__main__.py b/src/mintpy/__main__.py index 0492d52bc..ded28f330 100644 --- a/src/mintpy/__main__.py +++ b/src/mintpy/__main__.py @@ -5,8 +5,7 @@ ############################################################ -""" -Command line interface for MintPy. +"""Command line interface for MintPy. The Miami INsar Time-series software in PYthon (MintPy as /mint pai/) is an open-source package for Interferometric Synthetic Aperture Radar From 6d2720a7c21532e2625daa6ae7979338cd8fa2ec Mon Sep 17 00:00:00 2001 From: David Bekaert Date: Fri, 17 Apr 2026 14:20:45 +0200 Subject: [PATCH 21/24] =?UTF-8?q?fix:=20docstring=20D212=20violations=20?= =?UTF-8?q?=E2=80=94=20summary=20on=20first=20line=20(Codacy)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/mintpy/tropo_opera.py | 35 +++++++++++------------------------ 1 file changed, 11 insertions(+), 24 deletions(-) diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index c74a657d4..867080bff 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -51,10 +51,7 @@ def get_geom_lat_lon_bounds(geom_file): def get_opera_crop_indices(lat, lon, geom_file, pad_cells=3): - """Get OPERA cube crop indices from geometry bounds with fixed ±0.2 deg padding. - - The pad_cells argument is kept for backward compatibility and ignored. - """ + """Get OPERA cube crop indices from geometry bounds with fixed ±0.2 deg padding.""" if lat.size < 2 or lon.size < 2: raise ValueError('Invalid OPERA latitude/longitude arrays: size < 2') @@ -89,8 +86,7 @@ def get_opera_crop_indices(lat, lon, geom_file, pad_cells=3): def get_opera_height_crop_indices(height_levels, dem_min, dem_max): - """ - Get height level indices spanning the DEM range. + """Get height level indices spanning the DEM range. The vertical interpolation is linear, so only the levels that bracket [dem_min, dem_max] are needed — no extra buffer beyond that. @@ -131,8 +127,7 @@ def get_opera_height_crop_indices(height_levels, dem_min, dem_max): def read_opera_total_delay_cube(opera_file, geom_file, dem_range=None, pad_cells=3): - """ - Read and crop OPERA delay cube, returning total zenith delay. + """Read and crop OPERA delay cube, returning total zenith delay. Parameters: opera_file - str, path to OPERA netCDF file geom_file - str, path to MintPy geometry HDF5 file @@ -183,8 +178,7 @@ def read_opera_total_delay_cube(opera_file, geom_file, dem_range=None, pad_cells def get_geom_lat_lon_dem(geom_file): - """ - Read DEM and latitude/longitude grids from geometry file. + """Read DEM and latitude/longitude grids from geometry file. Returns: lat2d/lon2d/dem in 2D np.ndarray with same shape """ @@ -309,8 +303,7 @@ def calc_zenith_delay_from_opera_file( ############################################################################ def read_inps2date_time(inps): - """ - Read acquisition date/time info from input arguments. + """Read acquisition date/time info from input arguments. Parameters: inps - Namespace for input arguments Returns: date_list - list(str), acquisition dates in YYYYMMDD @@ -342,8 +335,7 @@ def read_inps2date_time(inps): def nearest_opera_product_time(date_str, utc_sec): - """ - Round one acquisition date/time to the nearest OPERA model datetime. + """Round one acquisition date/time to the nearest OPERA model datetime. Parameters: date_str - str, acquisition date in YYYYMMDD utc_sec - float, acquisition UTC time in seconds @@ -366,8 +358,7 @@ def nearest_opera_product_time(date_str, utc_sec): def get_opera_date_time_list(date_list, utc_sec): - """ - Build required OPERA date/hour list for all acquisitions. + """Build required OPERA date/hour list for all acquisitions. Parameters: date_list - list(str), acquisition dates in YYYYMMDD utc_sec - float, acquisition UTC time in seconds @@ -385,8 +376,7 @@ def get_opera_date_time_list(date_list, utc_sec): def get_expected_opera_file_patterns(opera_date_list, opera_hour_list): - """ - Create expected OPERA ZTD filename patterns for model times. + """Create expected OPERA ZTD filename patterns for model times. The production timestamp is left as wildcard. @@ -401,8 +391,7 @@ def get_expected_opera_file_patterns(opera_date_list, opera_hour_list): def get_opera_file_status(opera_date_list, opera_hour_list, opera_dir): - """ - Get matched OPERA files and missing model date/hour list. + """Get matched OPERA files and missing model date/hour list. Returns: expected_patterns - list(str) matched_files - dict(str, list(str)) @@ -531,8 +520,7 @@ def _get_product_filename(product): def dload_opera_file_subset(product, opera_dir, geom_file, dem_range=None, session=None): - """ - Download a spatially and vertically subsetted OPERA file. + """Download a spatially and vertically subsetted OPERA file. Uses ``fsspec`` + ``h5py`` to open the remote file via byte-range HTTP requests, reads only the required spatial/vertical slice, and writes the @@ -604,8 +592,7 @@ def dload_opera_file_subset(product, opera_dir, geom_file, dem_range=None, sessi def dload_opera_files(missing_date_hour_list, opera_dir, geom_file=None, dem_range=None): - """ - Download missing OPERA TROPO-ZENITH files via ASF Search. + """Download missing OPERA TROPO-ZENITH files via ASF Search. Uses per-day ASF queries constrained by OPERA collection IDs and date windows, then down-filters by required model-time tokens. From 057032fddc44b314f2342adec1468be9002e4805 Mon Sep 17 00:00:00 2001 From: David Bekaert Date: Fri, 17 Apr 2026 14:28:32 +0200 Subject: [PATCH 22/24] fix: convert multi-line docstrings to single-line to satisfy both D212+D213 (Codacy) --- src/mintpy/tropo_opera.py | 103 ++++---------------------------------- 1 file changed, 10 insertions(+), 93 deletions(-) diff --git a/src/mintpy/tropo_opera.py b/src/mintpy/tropo_opera.py index 867080bff..e221344f0 100644 --- a/src/mintpy/tropo_opera.py +++ b/src/mintpy/tropo_opera.py @@ -86,21 +86,7 @@ def get_opera_crop_indices(lat, lon, geom_file, pad_cells=3): def get_opera_height_crop_indices(height_levels, dem_min, dem_max): - """Get height level indices spanning the DEM range. - - The vertical interpolation is linear, so only the levels that bracket - [dem_min, dem_max] are needed — no extra buffer beyond that. - - For example, if dem_min=-5, dem_max=10 and height levels are - [-20, -10, 0, 10, 20], the returned slice covers [-10, 0, 10] - (the first level at-or-below dem_min through the first level - at-or-above dem_max). - - Parameters: height_levels - 1D np.ndarray, OPERA height levels (metres) - dem_min - float, minimum DEM elevation (metres) - dem_max - float, maximum DEM elevation (metres) - Returns: k0, k1 - int, start/end slice indices into *height_levels* - """ + """Get height level indices spanning the DEM range.""" h = np.asarray(height_levels, dtype=np.float64) n = len(h) if n < 2: @@ -127,17 +113,7 @@ def get_opera_height_crop_indices(height_levels, dem_min, dem_max): def read_opera_total_delay_cube(opera_file, geom_file, dem_range=None, pad_cells=3): - """Read and crop OPERA delay cube, returning total zenith delay. - - Parameters: opera_file - str, path to OPERA netCDF file - geom_file - str, path to MintPy geometry HDF5 file - dem_range - tuple(float, float) or None, (dem_min, dem_max) in - metres. When provided the height dimension is also - cropped to the levels bracketing [dem_min, dem_max]. - pad_cells - int, kept for backward compatibility (ignored) - Returns a dict with keys: - latitude, longitude, height, wet_delay, hydro_delay, total_delay, crop_info - """ + """Read and crop OPERA delay cube, returning total zenith delay.""" with h5py.File(opera_file, 'r') as fobj: req_dsets = ['latitude', 'longitude', 'height', 'wet_delay', 'hydrostatic_delay'] missing = [k for k in req_dsets if k not in fobj] @@ -178,10 +154,7 @@ def read_opera_total_delay_cube(opera_file, geom_file, dem_range=None, pad_cells def get_geom_lat_lon_dem(geom_file): - """Read DEM and latitude/longitude grids from geometry file. - - Returns: lat2d/lon2d/dem in 2D np.ndarray with same shape - """ + """Read DEM and latitude/longitude grids from geometry file.""" dem = readfile.read(geom_file, datasetName='height')[0].astype(np.float32) atr = readfile.read_attribute(geom_file) @@ -303,12 +276,7 @@ def calc_zenith_delay_from_opera_file( ############################################################################ def read_inps2date_time(inps): - """Read acquisition date/time info from input arguments. - - Parameters: inps - Namespace for input arguments - Returns: date_list - list(str), acquisition dates in YYYYMMDD - utc_sec - float, acquisition UTC time in seconds - """ + """Read acquisition date/time info from input arguments.""" if not inps.dis_file: raise ValueError('input displacement file is required to derive OPERA date/time info.') @@ -335,13 +303,7 @@ def read_inps2date_time(inps): def nearest_opera_product_time(date_str, utc_sec): - """Round one acquisition date/time to the nearest OPERA model datetime. - - Parameters: date_str - str, acquisition date in YYYYMMDD - utc_sec - float, acquisition UTC time in seconds - Returns: out_date - str, OPERA model date in YYYYMMDD - out_hour - str, OPERA model hour in HH format - """ + """Round one acquisition date/time to the nearest OPERA model datetime.""" acq_dt = datetime.strptime(date_str, '%Y%m%d') + timedelta(seconds=float(utc_sec)) candidates = [] @@ -358,13 +320,7 @@ def nearest_opera_product_time(date_str, utc_sec): def get_opera_date_time_list(date_list, utc_sec): - """Build required OPERA date/hour list for all acquisitions. - - Parameters: date_list - list(str), acquisition dates in YYYYMMDD - utc_sec - float, acquisition UTC time in seconds - Returns: opera_date_list - list(str), OPERA model date in YYYYMMDD - opera_hour_list - list(str), OPERA model hour in HH format - """ + """Build required OPERA date/hour list for all acquisitions.""" opera_date_list = [] opera_hour_list = [] for date_str in date_list: @@ -376,13 +332,7 @@ def get_opera_date_time_list(date_list, utc_sec): def get_expected_opera_file_patterns(opera_date_list, opera_hour_list): - """Create expected OPERA ZTD filename patterns for model times. - - The production timestamp is left as wildcard. - - Returns: expected_patterns - list(str), e.g. - OPERA_L4_TROPO-ZENITH_YYYYMMDDTHHMMSSZ_*_HRES_v* - """ + """Create expected OPERA ZTD filename patterns for model times.""" expected_patterns = sorted({ f'OPERA_L4_TROPO-ZENITH_{_model_time_token(date_str, hour_str)}_*_HRES_v*' for date_str, hour_str in zip(opera_date_list, opera_hour_list) @@ -391,12 +341,7 @@ def get_expected_opera_file_patterns(opera_date_list, opera_hour_list): def get_opera_file_status(opera_date_list, opera_hour_list, opera_dir): - """Get matched OPERA files and missing model date/hour list. - - Returns: expected_patterns - list(str) - matched_files - dict(str, list(str)) - missing_date_hour_list - list(tuple(str, str)) in (YYYYMMDD, HH) - """ + """Get matched OPERA files and missing model date/hour list.""" expected_patterns = get_expected_opera_file_patterns(opera_date_list, opera_hour_list) matched_files = {} @@ -520,20 +465,7 @@ def _get_product_filename(product): def dload_opera_file_subset(product, opera_dir, geom_file, dem_range=None, session=None): - """Download a spatially and vertically subsetted OPERA file. - - Uses ``fsspec`` + ``h5py`` to open the remote file via byte-range HTTP - requests, reads only the required spatial/vertical slice, and writes the - subset to a local HDF5 file that is compatible with - :func:`read_opera_total_delay_cube`. - - Parameters: product - ASF search product object - opera_dir - str, local directory for output files - geom_file - str, path to MintPy geometry file (for bounding box) - dem_range - tuple(float, float) or None, (min, max) DEM in metres - session - ASF session object (used for authentication cookies) - Returns: out_file - str, path to the written local subset file, or None - """ + """Download a spatially and vertically subsetted OPERA file.""" import fsspec url = _get_product_url(product) @@ -592,22 +524,7 @@ def dload_opera_file_subset(product, opera_dir, geom_file, dem_range=None, sessi def dload_opera_files(missing_date_hour_list, opera_dir, geom_file=None, dem_range=None): - """Download missing OPERA TROPO-ZENITH files via ASF Search. - - Uses per-day ASF queries constrained by OPERA collection IDs and - date windows, then down-filters by required model-time tokens. - - When *geom_file* is provided (and ``fsspec`` is installed), files are - downloaded as spatial/vertical subsets using byte-range HTTP requests, - which significantly reduces bandwidth and disk usage. Falls back to - full-file download when subsetting is not possible. - - Parameters: missing_date_hour_list - list of (date_str, hour_str) tuples - opera_dir - str, local directory for OPERA files - geom_file - str or None, MintPy geometry file - dem_range - tuple(float, float) or None, - (min, max) DEM elevation in metres - """ + """Download missing OPERA TROPO-ZENITH files via ASF Search.""" if len(missing_date_hour_list) == 0: return From adcac571dede6f55c01bd790414cc51163d075cd Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Fri, 29 May 2026 13:39:21 +0800 Subject: [PATCH 23/24] cli/tropo_opera: minor print out msg re-org --- src/mintpy/cli/tropo_opera.py | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/src/mintpy/cli/tropo_opera.py b/src/mintpy/cli/tropo_opera.py index 1250109ff..8a24326c1 100755 --- a/src/mintpy/cli/tropo_opera.py +++ b/src/mintpy/cli/tropo_opera.py @@ -13,13 +13,8 @@ ############################################################################ REFERENCE = """references: - Bekaert, D. et al., OPERA L4 Tropospheric Zenith Delay products - derived from ECMWF HRES model (https://www.earthdata.nasa.gov/data/catalog/asf-opera-l4-tropo-zenith-v1-1). -""" - -DIR_DEMO = """--dir ./OPERA - Path to the directory for downloading and storing OPERA tropospheric delay - products. Files are automatically downloaded from ASF on demand. + Bekaert, D. et al., OPERA L4 Tropospheric Zenith Delay products derived from ECMWF HRES model + (https://www.earthdata.nasa.gov/data/catalog/asf-opera-l4-tropo-zenith-v1-1). """ EXAMPLE = """example: @@ -29,8 +24,9 @@ def create_parser(subparsers=None): - synopsis = 'Tropospheric correction using OPERA Zenith Tropospheric Delays from ECMWF HRES data (https://www.earthdata.nasa.gov/data/catalog/asf-opera-l4-tropo-zenith-v1-1)' - epilog = REFERENCE + '\n' + DIR_DEMO + '\n' + EXAMPLE + synopsis = 'Tropospheric correction using OPERA Zenith Tropospheric Delays from ECMWF HRES data ' + synopsis += '\n(https://www.earthdata.nasa.gov/data/catalog/asf-opera-l4-tropo-zenith-v1-1)' + epilog = REFERENCE + '\n' + EXAMPLE name = __name__.split('.')[-1] parser = create_argument_parser( name, synopsis=synopsis, description=synopsis, epilog=epilog, subparsers=subparsers) @@ -39,7 +35,8 @@ def create_parser(subparsers=None): parser.add_argument('-g', '--geom', dest='geom_file', required=True, help='geometry file.') parser.add_argument('--dir','--opera-dir', dest='opera_dir', default='./OPERA', - help='directory to store downloaded OPERA delays data (default: %(default)s).') + help='directory to store downloaded OPERA tropospheric delay products, \n' + + 'which are automatically downloaded from ASF on demand (default: %(default)s).') parser.add_argument('-o', dest='cor_dis_file', help='Output file name for tropospheric corrected timeseries.') From b15a570b696ef93383bf84d29c1370c371ac434e Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Fri, 29 May 2026 13:41:28 +0800 Subject: [PATCH 24/24] pre-commit fix --- src/mintpy/cli/tropo_opera.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mintpy/cli/tropo_opera.py b/src/mintpy/cli/tropo_opera.py index 8a24326c1..155b43f81 100755 --- a/src/mintpy/cli/tropo_opera.py +++ b/src/mintpy/cli/tropo_opera.py @@ -13,7 +13,7 @@ ############################################################################ REFERENCE = """references: - Bekaert, D. et al., OPERA L4 Tropospheric Zenith Delay products derived from ECMWF HRES model + Bekaert, D. et al., OPERA L4 Tropospheric Zenith Delay products derived from ECMWF HRES model (https://www.earthdata.nasa.gov/data/catalog/asf-opera-l4-tropo-zenith-v1-1). """