276 lines
11 KiB
Python
276 lines
11 KiB
Python
#! /usr/bin/env python
|
||
#################################################################
|
||
### This program is part of PyINT v2.1 ###
|
||
### Copy Right (c): 2017-2019, Yunmeng Cao ###
|
||
### Author: Yunmeng Cao ###
|
||
### Modified: 2026, Z. Zou - HyP3 UTM GeoTIFF output ###
|
||
### Contact : ymcmrs@gmail.com ###
|
||
#################################################################
|
||
|
||
import numpy as np
|
||
import os
|
||
import sys
|
||
import argparse
|
||
import time
|
||
from datetime import datetime
|
||
|
||
from pyint import _utils as ut
|
||
|
||
|
||
INTRODUCTION = '''
|
||
-------------------------------------------------------------------
|
||
将 geocode_gamma.py 已地理编码的 GAMMA 二进制产品转换为
|
||
ASF HyP3 兼容的 UTM GeoTIFF 格式, 用于 MintPy 时间序列分析。
|
||
|
||
本脚本 **不做任何地理编码或计算**, 仅执行:
|
||
1. data2geotiff : GAMMA 二进制 (EQA) → GeoTIFF (EQA)
|
||
2. gdalwarp : GeoTIFF (EQA) → GeoTIFF (UTM)
|
||
3. rasterio : 生成水体掩膜
|
||
|
||
所有地理编码产品由 geocode_gamma.py 提供:
|
||
amp, corr, dem, unw_phase, wrapped_phase,
|
||
los_disp, vert_disp, lv_theta, lv_phi
|
||
'''
|
||
|
||
EXAMPLE = '''
|
||
Usage:
|
||
hyp3format_gamma.py projectName ifgPair
|
||
hyp3format_gamma.py shanghaiT171F128S1A 20241105-20241117
|
||
hyp3format_gamma.py shanghaiT171F128S1A 20241105-20241117 --output-dir /path/to/output
|
||
-------------------------------------------------------------------
|
||
'''
|
||
|
||
|
||
def cmdLineParse():
|
||
parser = argparse.ArgumentParser(
|
||
description='Convert geocoded GAMMA binaries to HyP3 UTM GeoTIFF.',
|
||
formatter_class=argparse.RawTextHelpFormatter,
|
||
epilog=INTRODUCTION + '\n' + EXAMPLE)
|
||
parser.add_argument('projectName', help='项目名称')
|
||
parser.add_argument('ifgPair', help='干涉图对 (如 20241105-20241117)')
|
||
parser.add_argument('--output-dir', dest='output_dir', default=None,
|
||
help='输出目录. 默认: projectDir/Hyp3Products/')
|
||
inps = parser.parse_args()
|
||
return inps
|
||
|
||
|
||
def get_utm_epsg(lon, lat):
|
||
"""根据中心经纬度自动计算 UTM 投影的 EPSG 代码"""
|
||
zone = int((lon + 180) / 6) + 1
|
||
return 32600 + zone if lat >= 0 else 32700 + zone
|
||
|
||
|
||
def run_cmd(cmd, desc=''):
|
||
"""执行 shell 命令并检查返回值"""
|
||
if desc:
|
||
print(f" [{desc}] {cmd}")
|
||
else:
|
||
print(f" {cmd}")
|
||
ret = os.system(cmd)
|
||
if ret != 0:
|
||
print(f" WARNING: 返回非零退出码 {ret}")
|
||
return ret
|
||
|
||
|
||
def main(argv):
|
||
|
||
start_time = time.time()
|
||
inps = cmdLineParse()
|
||
projectName = inps.projectName
|
||
ifgPair = inps.ifgPair
|
||
|
||
Mdate = ifgPair.split('-')[0]
|
||
Sdate = ifgPair.split('-')[1]
|
||
|
||
# 验证日期并计算时间基线
|
||
d1 = datetime.strptime(Mdate, '%Y%m%d')
|
||
d2 = datetime.strptime(Sdate, '%Y%m%d')
|
||
interval = abs((d2 - d1).days)
|
||
|
||
# 读取模板参数
|
||
scratchDir = os.getenv('SCRATCHDIR')
|
||
templateDir = os.getenv('TEMPLATEDIR')
|
||
templateFile = templateDir + '/' + projectName + '.template'
|
||
templateDict = ut.update_template(templateFile)
|
||
|
||
rlks = templateDict['range_looks']
|
||
azlks = templateDict['azimuth_looks']
|
||
masterDate = templateDict['masterDate']
|
||
|
||
# 项目目录
|
||
projectDir = scratchDir + '/' + projectName
|
||
ifgDir = projectDir + '/ifgrams'
|
||
demDir = projectDir + '/DEM'
|
||
|
||
workDir = ifgDir + '/' + ifgPair
|
||
if not os.path.isdir(workDir):
|
||
print(f"ERROR: 干涉图目录不存在: {workDir}")
|
||
sys.exit(1)
|
||
|
||
# 输出目录
|
||
outputDir = inps.output_dir or (projectDir + '/Hyp3Products')
|
||
productDir = os.path.join(outputDir, ifgPair)
|
||
os.makedirs(productDir, exist_ok=True)
|
||
|
||
# ============================================================
|
||
# DEM 参数 (直接引用, 不复制不删除)
|
||
# ============================================================
|
||
DEMpar = demDir + '/' + masterDate + '_' + rlks + 'rlks.utm.dem.par'
|
||
|
||
nWidthDEM = ut.read_gamma_par(DEMpar, 'read', 'width')
|
||
nLineDEM = ut.read_gamma_par(DEMpar, 'read', 'nlines')
|
||
|
||
# 自动检测 UTM 投影
|
||
corner_lat = float(ut.read_gamma_par(DEMpar, 'read', 'corner_lat').split()[0])
|
||
corner_lon = float(ut.read_gamma_par(DEMpar, 'read', 'corner_lon').split()[0])
|
||
post_lat = float(ut.read_gamma_par(DEMpar, 'read', 'post_lat').split()[0])
|
||
post_lon = float(ut.read_gamma_par(DEMpar, 'read', 'post_lon').split()[0])
|
||
center_lat = corner_lat + float(nLineDEM) / 2 * post_lat
|
||
center_lon = corner_lon + float(nWidthDEM) / 2 * post_lon
|
||
utm_epsg = get_utm_epsg(center_lon, center_lat)
|
||
|
||
# HyP3 文件名前缀 (含 YYYYMMDDTHHMMSS, 兼容 MintPy 日期解析)
|
||
hyp3_prefix = f"{projectName}_{Mdate}T000000_{Sdate}T000000"
|
||
|
||
print("\n" + "=" * 70)
|
||
print(f"HyP3 格式转换: {projectName} / {ifgPair}")
|
||
print(f"时间基线: {interval} 天, DEM: {nWidthDEM}×{nLineDEM}")
|
||
print(f"中心: {center_lat:.4f}°N, {center_lon:.4f}°E → EPSG:{utm_epsg}")
|
||
print(f"输出: {productDir}")
|
||
print("=" * 70)
|
||
|
||
# ============================================================
|
||
# 已地理编码 (EQA) GAMMA 二进制文件 (由 geocode_gamma.py 生成)
|
||
# ============================================================
|
||
geo_unw = workDir + '/geo_' + ifgPair + '_' + rlks + 'rlks.diff_filt.unw'
|
||
geo_wrapped_pha = workDir + '/geo_' + ifgPair + '_' + rlks + 'rlks.diff_filt.pha'
|
||
geo_los_disp = workDir + '/geo_' + ifgPair + '_' + rlks + 'rlks.los_disp'
|
||
geo_vert_disp = workDir + '/geo_' + ifgPair + '_' + rlks + 'rlks.vert_disp'
|
||
geo_amp = workDir + '/geo_' + masterDate + '_' + rlks + 'rlks.amp'
|
||
geo_hgt = workDir + '/geo_' + masterDate + '_' + rlks + 'rlks.hgt'
|
||
lv_theta = workDir + '/lv_theta'
|
||
lv_phi = workDir + '/lv_phi'
|
||
|
||
# 相干性可能用 Pair 或 masterDate 命名
|
||
geo_cor = workDir + '/geo_' + ifgPair + '_' + rlks + 'rlks.diff_filt.cor'
|
||
if not os.path.exists(geo_cor):
|
||
geo_cor = workDir + '/geo_' + masterDate + '_' + rlks + 'rlks.diff_filt.cor'
|
||
|
||
# ============================================================
|
||
# data2geotiff (EQA binary → EQA GeoTIFF) → gdalwarp (→ UTM)
|
||
# ============================================================
|
||
print("\n--- GAMMA 二进制 → EQA GeoTIFF → UTM GeoTIFF ---")
|
||
|
||
product_map = [
|
||
('amp', geo_amp, 'bilinear'),
|
||
('corr', geo_cor, 'bilinear'),
|
||
('dem', geo_hgt, 'bilinear'),
|
||
('unw_phase', geo_unw, 'bilinear'),
|
||
('wrapped_phase', geo_wrapped_pha, 'bilinear'),
|
||
('los_disp', geo_los_disp, 'bilinear'),
|
||
('vert_disp', geo_vert_disp, 'bilinear'),
|
||
('lv_theta', lv_theta, 'bilinear'),
|
||
('lv_phi', lv_phi, 'bilinear'),
|
||
]
|
||
|
||
generated = {}
|
||
|
||
for name, geo_bin, interp in product_map:
|
||
utm_tif = os.path.join(productDir, f'{hyp3_prefix}_{name}.tif')
|
||
|
||
# 跳过已存在且非空的输出
|
||
if os.path.exists(utm_tif) and os.path.getsize(utm_tif) > 0:
|
||
print(f" {name}: 已存在，跳过")
|
||
generated[name] = utm_tif
|
||
continue
|
||
|
||
if not os.path.exists(geo_bin):
|
||
print(f" {name}: 源文件不存在 ({os.path.basename(geo_bin)})，跳过")
|
||
continue
|
||
|
||
eqa_tif = workDir + '/' + name + '_hyp3_eqa.tif'
|
||
|
||
# GAMMA 二进制 → EQA GeoTIFF
|
||
run_cmd(f'data2geotiff {DEMpar} {geo_bin} 2 {eqa_tif}', f'{name} → EQA')
|
||
|
||
if not os.path.exists(eqa_tif):
|
||
print(f" WARNING: {name} EQA GeoTIFF 生成失败")
|
||
continue
|
||
|
||
# EQA → UTM (gdalwarp)
|
||
run_cmd(f'gdalwarp -t_srs EPSG:{utm_epsg} -r {interp} '
|
||
f'-co COMPRESS=LZW -overwrite {eqa_tif} {utm_tif}',
|
||
f'{name} → UTM')
|
||
|
||
# 清理 EQA 中间文件
|
||
if os.path.exists(eqa_tif):
|
||
os.remove(eqa_tif)
|
||
|
||
if os.path.exists(utm_tif) and os.path.getsize(utm_tif) > 0:
|
||
generated[name] = utm_tif
|
||
|
||
# ============================================================
|
||
# 生成水体掩膜 (相干性阈值)
|
||
# ============================================================
|
||
print("\n--- 水体掩膜 ---")
|
||
water_mask_tif = os.path.join(productDir, f'{hyp3_prefix}_water_mask.tif')
|
||
corr_tif = generated.get('corr')
|
||
|
||
if corr_tif and os.path.exists(corr_tif) and not os.path.exists(water_mask_tif):
|
||
try:
|
||
import rasterio
|
||
with rasterio.open(corr_tif) as src:
|
||
corr_data = src.read(1)
|
||
mask = np.ones_like(corr_data, dtype=np.uint8)
|
||
mask[corr_data < 0.05] = 0
|
||
mask[np.isnan(corr_data)] = 0
|
||
|
||
profile = src.profile.copy()
|
||
profile.update(dtype=rasterio.uint8, count=1, compress='lzw', nodata=None)
|
||
with rasterio.open(water_mask_tif, 'w', **profile) as dst:
|
||
dst.write(mask, 1)
|
||
generated['water_mask'] = water_mask_tif
|
||
print(f" 水体掩膜已生成")
|
||
except Exception as e:
|
||
print(f" WARNING: 水体掩膜生成失败: {e}")
|
||
elif os.path.exists(water_mask_tif):
|
||
generated['water_mask'] = water_mask_tif
|
||
print(f" 水体掩膜已存在，跳过")
|
||
|
||
# ============================================================
|
||
# 生成入射角图 (inc_map / inc_map_ell, 复制自 lv_theta)
|
||
# ============================================================
|
||
print("\n--- 入射角图 ---")
|
||
inc_map_tif = os.path.join(productDir, f'{hyp3_prefix}_inc_map.tif')
|
||
inc_map_ell_tif = os.path.join(productDir, f'{hyp3_prefix}_inc_map_ell.tif')
|
||
lv_theta_tif = generated.get('lv_theta')
|
||
|
||
if lv_theta_tif and os.path.exists(lv_theta_tif):
|
||
if not os.path.exists(inc_map_tif):
|
||
os.system(f'cp {lv_theta_tif} {inc_map_tif}')
|
||
generated['inc_map'] = inc_map_tif
|
||
print(f" inc_map 已生成 (= lv_theta)")
|
||
if not os.path.exists(inc_map_ell_tif):
|
||
os.system(f'cp {lv_theta_tif} {inc_map_ell_tif}')
|
||
generated['inc_map_ell'] = inc_map_ell_tif
|
||
print(f" inc_map_ell 已生成 (近似 lv_theta)")
|
||
|
||
# ============================================================
|
||
# 输出摘要
|
||
# ============================================================
|
||
print("\n" + "=" * 70)
|
||
print(f"HyP3 格式转换完成: {ifgPair}")
|
||
print(f"共 {len(generated)} 个产品:")
|
||
for name, path in sorted(generated.items()):
|
||
if os.path.exists(path):
|
||
size_mb = os.path.getsize(path) / 1024 / 1024
|
||
print(f" {name:20s} ({size_mb:.1f} MB)")
|
||
print("=" * 70)
|
||
|
||
ut.print_process_time(start_time, time.time())
|
||
sys.exit(0)
|
||
|
||
|
||
if __name__ == '__main__':
|
||
main(sys.argv[:])
|