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) 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/requirements.txt b/requirements.txt index cce66eee5..c3243be52 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,9 +1,11 @@ pip argcomplete +asf_search>=12.0.0 cartopy cvxopt dask>=1.0 dask-jobqueue>=0.3 +fsspec[http] h5py joblib lxml 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/cli/tropo_opera.py b/src/mintpy/cli/tropo_opera.py new file mode 100755 index 000000000..155b43f81 --- /dev/null +++ b/src/mintpy/cli/tropo_opera.py @@ -0,0 +1,107 @@ +#!/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: + 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: + 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 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) + 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 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.') + + 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/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/tropo_opera.py b/src/mintpy/tropo_opera.py new file mode 100644 index 000000000..e221344f0 --- /dev/null +++ b/src/mintpy/tropo_opera.py @@ -0,0 +1,921 @@ +############################################################ +# Program is part of MintPy # +# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi # +# Author: David Bekaert, March 2026 # +############################################################ + + +import glob +import netrc +import os +import re +from datetime import datetime, timedelta + +import h5py +import numpy as np + +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 = [ + 'C3717139408-ASF', + 'C1273910987-ASF', +] + + +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' + + +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 + + +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.""" + 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 + + +def get_opera_height_crop_indices(height_levels, dem_min, dem_max): + """Get height level indices spanning the DEM range.""" + 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.""" + 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 + ) + + # 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, 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_crop, + 'wet_delay': wet, + 'hydro_delay': hydro, + 'total_delay': total, + 'crop_info': crop_info, + } + + +def get_geom_lat_lon_dem(geom_file): + """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) + + 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 _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, lat2d, lon2d, dem, pad_cells=3 +): + """Calculate 2D zenith tropospheric delay map intersected with DEM.""" + from scipy.ndimage import map_coordinates + + # 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.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) + + # 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 + + # 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) + + # 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) + ztd[~np.isfinite(dem)] = np.nan + + return ztd, cube + + +############################################################################ +def read_inps2date_time(inps): + """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.') + + 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(inps.dis_file).get_date_list() + elif ftype == '.unw': + date_list = ptime.yyyymmdd(atr['DATE12'].split('-')) + else: + raise ValueError(f'un-supported displacement file type: {ftype}') + + 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.""" + 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.""" + 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.""" + 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.""" + 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 + + +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.""" + 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.""" + 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 + + # Determine whether remote subsetting is possible + _can_subset = geom_file is not None + if _can_subset: + try: + 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}) + 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 + + 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} (done)') + + 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).""" + + def get_dataset_size(fname): + atr = readfile.read_attribute(fname) + return (atr['LENGTH'], atr['WIDTH']) + + print('update mode: ON') + print(f'output file: {tropo_file}') + flag = 'skip' + + 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('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.') + + print(f'run or skip: {flag}') + return flag + + +############################################################################ +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'] + + 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}') + + 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']) + + 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: + raise RuntimeError( + f'No acquisition dates are supported by OPERA (requires >= {OPERA_MIN_ACQ_DATE}).' + ) + + 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 ' + 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.' + ) + + if _run_or_skip_opera(opera_files, used_dates, tropo_file, geom_file) == 'skip': + return tropo_file + + 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) + + 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) + + # 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, lat2d, lon2d, dem, pad_cells=3 + ) + delay_los = _project_zenith_to_los(ztd, inc_angle) + + block = [i, i + 1, 0, length, 0, width] + writefile.write_hdf5_block( + tropo_file, + data=delay_los, + datasetName='timeseries', + block=block, + print_msg=False, + ) + + prog_bar.update(i + 1, suffix=os.path.basename(opera_files[i])) + prog_bar.close() + + return tropo_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: + # 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 as exc: + print(f' WARNING: could not read DEM for subsetting: {exc}') + + 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, + 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_delay_timeseries( + tropo_file=inps.tropo_file, + dis_file=inps.dis_file, + geom_file=inps.geom_file, + opera_dir=inps.opera_dir, + ) + + 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 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',