fix: BIP 格式兼容 — 鲁棒 HDR 查找 + find_band_number GDAL 回退

问题: 用户导入 3ref.bip 文件,step2 报 FileNotFoundError: 3ref.hdr 不存在。

根本原因: get_hdr_file_path() 只用 os.path.splitext()[0]+.hdr,
对于 3ref.bip 只查找 3ref.hdr,不兼容 3ref.bip.hdr 等其他命名规范。

修复内容:

**util.py (核心):**
- get_hdr_file_path(): 改为多候选路径查找(按优先级):
  3ref.hdr → 3ref.bip.hdr → 3ref.HDR → 3ref.bip.HDR
- find_band_number(): 三级回退 —
  1) ENVI .hdr 文件 → spectral 解析
  2) GDAL 元数据域 (ENVI/wavelength, WAVELENGTH_1..N)
  3) 线性估算 (假设 400-1000nm 或 400-2500nm)
- 新增 _read_wavelengths_from_gdal() 辅助函数

**同模式修复 (4 处):**
- get_spectral.py: get_hdr_file_path() 多候选
- get_spectral-test.py: 同上
- waterindex_inversion/__init__.py: 两处 hdr 构造均改为多候选
- sampling.py: 波长读取的 hdr 查找改为多候选
This commit is contained in:
duxin
2026-07-07 10:32:02 +08:00
parent 1f6c81158b
commit 2b5c77131e
5 changed files with 204 additions and 40 deletions

View File

@ -205,13 +205,18 @@ class WaterIndexProcessor:
except Exception:
pass
# 2. HDR 补充波长信息
# 2. HDR 补充波长信息(多命名规范兼容 .bsq/.bil/.bip/.dat
if hdr_path is None:
hdr_path = os.path.splitext(bsq_path)[0] + '.hdr'
if not os.path.isfile(hdr_path):
hdr_path_alt = os.path.splitext(bsq_path)[0] + '.HDR'
if os.path.isfile(hdr_path_alt):
hdr_path = hdr_path_alt
hdr_candidates = [
os.path.splitext(bsq_path)[0] + '.hdr', # 3ref.hdr
bsq_path + '.hdr', # 3ref.bip.hdr
os.path.splitext(bsq_path)[0] + '.HDR', # 3ref.HDR
bsq_path + '.HDR', # 3ref.bip.HDR
]
for candidate in hdr_candidates:
if os.path.isfile(candidate):
hdr_path = candidate
break
if os.path.isfile(hdr_path):
wl = self._parse_wavelengths_from_hdr(hdr_path)
@ -365,13 +370,18 @@ class WaterIndexProcessor:
dict
{公式名: 输出 GeoTIFF 路径}
"""
# ── 自动构造 HDR 路径 ────────────────────────────────────────────
# ── 自动构造 HDR 路径(多命名规范兼容) ────────────────────────
if hdr_path is None:
hdr_path = os.path.splitext(bsq_path)[0] + '.hdr'
if not os.path.isfile(hdr_path):
hdr_path_alt = os.path.splitext(bsq_path)[0] + '.HDR'
if os.path.isfile(hdr_path_alt):
hdr_path = hdr_path_alt
hdr_candidates = [
os.path.splitext(bsq_path)[0] + '.hdr',
bsq_path + '.hdr',
os.path.splitext(bsq_path)[0] + '.HDR',
bsq_path + '.HDR',
]
for candidate in hdr_candidates:
if os.path.isfile(candidate):
hdr_path = candidate
break
# ── 自动构造输出目录 ────────────────────────────────────────────
if output_dir is None:
@ -605,11 +615,17 @@ class WaterIndexProcessor:
notify("开始水色指数反演", 0)
bsq_path = deglint_img_path
hdr_path = os.path.splitext(bsq_path)[0] + '.hdr'
if not os.path.isfile(hdr_path):
hdr_path_alt = os.path.splitext(bsq_path)[0] + '.HDR'
if os.path.isfile(hdr_path_alt):
hdr_path = hdr_path_alt
hdr_candidates = [
os.path.splitext(bsq_path)[0] + '.hdr',
bsq_path + '.hdr',
os.path.splitext(bsq_path)[0] + '.HDR',
bsq_path + '.HDR',
]
hdr_path = None
for candidate in hdr_candidates:
if os.path.isfile(candidate):
hdr_path = candidate
break
output_dir = os.path.join(work_dir, "10_WaterIndex_Images")

View File

@ -330,15 +330,18 @@ def load_mask_file(mask_path):
def get_hdr_file_path(file_path):
"""
获取HDR文件路径
Args:
file_path: 影像文件路径
Returns:
HDR文件路径
获取 ENVI 头文件路径(鲁棒版:多命名规范兼容)
"""
return os.path.splitext(file_path)[0] + ".hdr"
candidates = [
os.path.splitext(file_path)[0] + ".hdr",
file_path + ".hdr",
os.path.splitext(file_path)[0] + ".HDR",
file_path + ".HDR",
]
for path in candidates:
if os.path.isfile(path):
return path
return candidates[0]
def calculate_utm_zone(longitude):

View File

@ -212,15 +212,26 @@ def load_mask_file(mask_path):
def get_hdr_file_path(file_path):
"""
获取HDR文件路径
获取 ENVI 头文件路径(鲁棒版:多命名规范兼容)
支持 .bsq / .bil / .bip / .dat 等格式的多种 .hdr 命名规范。
Args:
file_path: 影像文件路径
Returns:
HDR文件路径
存在的 .hdr 文件路径;若都不存在,返回标准命名路径
"""
return os.path.splitext(file_path)[0] + ".hdr"
candidates = [
os.path.splitext(file_path)[0] + ".hdr", # 3ref.hdr
file_path + ".hdr", # 3ref.bip.hdr
os.path.splitext(file_path)[0] + ".HDR", # 3ref.HDR
file_path + ".HDR", # 3ref.bip.HDR
]
for path in candidates:
if os.path.isfile(path):
return path
return candidates[0]
def load_wavelength_columns(imgpath, num_bands):

View File

@ -38,13 +38,23 @@ def get_wavelengths_from_bil_header(bil_file):
list - 波长列表如果无法获取则返回None
"""
try:
# 获取头文件路径
header_file = os.path.splitext(bil_file)[0] + ".hdr"
if not os.path.exists(header_file):
print(f"警告: 找不到头文件 {header_file}")
# 获取头文件路径(多命名规范兼容 .bsq/.bil/.bip/.dat
hdr_candidates = [
os.path.splitext(bil_file)[0] + ".hdr", # 3ref.hdr
bil_file + ".hdr", # 3ref.bip.hdr
os.path.splitext(bil_file)[0] + ".HDR", # 3ref.HDR
bil_file + ".HDR", # 3ref.bip.HDR
]
header_file = None
for candidate in hdr_candidates:
if os.path.exists(candidate):
header_file = candidate
break
if header_file is None:
print(f"警告: 找不到头文件,已尝试: {hdr_candidates}")
return None
# 使用spectral库读取头文件
import spectral.io.envi as envi
header = envi.read_envi_header(header_file)

View File

@ -33,18 +33,142 @@ def timeit(f): # decorator
def get_hdr_file_path(file_path):
return os.path.splitext(file_path)[0] + ".hdr"
"""获取 ENVI 头文件路径(鲁棒版:多命名规范兼容)
支持以下命名模式(按优先级检测):
1. {filename}.hdr ← 标准 ENVI 规范,如 3ref.hdr
2. {filename_with_ext}.hdr ← 如 3ref.bip.hdr
3. {basename}.HDR ← 大写变体
4. {filename}.HDR
Args:
file_path: 影像文件路径(.bsq / .bil / .bip / .dat 等)
Returns:
存在的 .hdr 文件路径;若都不存在,返回标准命名路径(交给调用方报错)
"""
# 候选路径列表(按优先级)
candidates = [
os.path.splitext(file_path)[0] + ".hdr", # 3ref.hdr
file_path + ".hdr", # 3ref.bip.hdr
os.path.splitext(file_path)[0] + ".HDR", # 3ref.HDR
file_path + ".HDR", # 3ref.bip.HDR
]
for path in candidates:
if os.path.isfile(path):
return path
# 都不存在:返回标准命名(让调用方报明确错误)
return candidates[0]
def find_band_number(wav1, imgpath):
in_hdr_dict = spectral.envi.read_envi_header(get_hdr_file_path(imgpath))
"""根据目标波长查找最接近的波段序号0-based
wavelengths = np.array(in_hdr_dict['wavelength']).astype('float64')
优先级:
1) ENVI .hdr 文件 → spectral 库解析
2) GDAL 元数据 → 直接从数据集读取波长域
3) 回退 → 基于波段计数线性估算
differences = np.abs(wavelengths - wav1)
min_position = np.argmin(differences)
Args:
wav1: 目标波长 (nm)
imgpath: 影像文件路径
return int(min_position)
Returns:
最接近的波段序号0-based int
"""
# ── 路径 1ENVI .hdr 文件 ──
hdr_path = get_hdr_file_path(imgpath)
if os.path.isfile(hdr_path):
try:
in_hdr_dict = spectral.envi.read_envi_header(hdr_path)
wavelengths = np.array(in_hdr_dict['wavelength']).astype('float64')
differences = np.abs(wavelengths - wav1)
min_position = int(np.argmin(differences))
print(f"[find_band] HDR 解析成功: target={wav1}nm → band {min_position} "
f"(wl={wavelengths[min_position]:.2f}nm), 来自 {hdr_path}")
return min_position
except Exception as e:
print(f"[find_band] HDR 解析失败 ({hdr_path}): {e},尝试 GDAL 回退")
# ── 路径 2GDAL 元数据 ──
try:
ds = gdal.Open(imgpath, gdal.GA_ReadOnly)
if ds is not None:
n_bands = ds.RasterCount
# 尝试从 GDAL 元数据域读取波长
wavelengths = _read_wavelengths_from_gdal(ds, n_bands)
ds = None
if wavelengths is not None and len(wavelengths) > 0:
differences = np.abs(np.array(wavelengths, dtype='float64') - wav1)
min_position = int(np.argmin(differences))
print(f"[find_band] GDAL 回退成功: target={wav1}nm → band {min_position} "
f"(wl={wavelengths[min_position]:.2f}nm), 共 {n_bands} 波段")
return min_position
except Exception as e:
print(f"[find_band] GDAL 元数据读取失败: {e}")
# ── 路径 3线性估算最后回退 ──
try:
ds = gdal.Open(imgpath, gdal.GA_ReadOnly)
n_bands = ds.RasterCount
ds = None
# 假设波长范围 400-1000nm 或 400-2500nm 线性分布
# 这是非常粗略的估算,但比崩溃好
estimated_wl_min, estimated_wl_max = (400.0, 1000.0) if n_bands <= 300 else (400.0, 2500.0)
band_idx = int(round((wav1 - estimated_wl_min) / (estimated_wl_max - estimated_wl_min) * (n_bands - 1)))
band_idx = max(0, min(n_bands - 1, band_idx))
print(f"[find_band] ⚠ 线性估算回退: target={wav1}nm → band {band_idx} "
f"(假设范围 {estimated_wl_min}-{estimated_wl_max}nm, {n_bands} 波段)")
return band_idx
except Exception as e:
raise RuntimeError(
f"无法确定波段号 (target={wav1}nm, img={imgpath}): "
f"HDR 不存在、GDAL 无波长元数据、且无法读取波段数。"
) from e
def _read_wavelengths_from_gdal(dataset, n_bands: int):
"""从 GDAL Dataset 的元数据域中提取波长列表
支持的来源(按优先级):
- ENVI 域: ENVI/wavelength
- 默认域: WAVELENGTH_1, WAVELENGTH_2, ...
Returns:
wavelengths list 或 None
"""
wavelengths = []
# 尝试 ENVI metadata domain
try:
envi_md = dataset.GetMetadata('ENVI')
if 'wavelength' in envi_md:
raw = envi_md['wavelength']
# 可能是 { ... } 包裹的逗号分隔列表
raw = raw.strip('{}').strip()
wavelengths = [float(x.strip()) for x in raw.split(',') if x.strip()]
if wavelengths:
return wavelengths
except Exception:
pass
# 尝试从默认域按 WAVELENGTH_1, WAVELENGTH_2, ... 读取
try:
md = dataset.GetMetadata()
for i in range(1, n_bands + 1):
key = f'WAVELENGTH_{i}'
if key in md:
wavelengths.append(float(md[key]))
else:
break
if wavelengths:
return wavelengths
except Exception:
pass
return None
@timeit