239 lines
8.5 KiB
Python
239 lines
8.5 KiB
Python
import argparse
|
||
import pandas as pd
|
||
import numpy as np
|
||
|
||
|
||
def find_nearest_bandnumber(wavelength, wavelengths):
|
||
differences = np.abs(wavelengths - wavelength)
|
||
min_position = np.argmin(differences)
|
||
return int(min_position)
|
||
|
||
|
||
def deglint(wavelengths, reflectance, start_wave=730, end_wave=750):
|
||
"""去除水表耀光"""
|
||
s = find_nearest_bandnumber(start_wave, wavelengths)
|
||
e = find_nearest_bandnumber(end_wave, wavelengths)
|
||
|
||
glint = np.mean(reflectance[s:e + 1])
|
||
deglinted_ref = reflectance - glint
|
||
deglinted_ref[deglinted_ref < 0] = 0
|
||
|
||
return deglinted_ref
|
||
|
||
|
||
def replace_non_valid_number(x):
|
||
"""检查并替换无效数字"""
|
||
if np.isnan(x) or np.isinf(x) or x < 0 or x > 65535:
|
||
return 0.0
|
||
return x
|
||
|
||
|
||
def safe_divide(numerator, denominator, default_value=0.0):
|
||
"""安全除法，避免除零错误"""
|
||
if denominator == 0 or np.isnan(denominator) or np.isinf(denominator):
|
||
return default_value
|
||
return numerator / denominator
|
||
|
||
|
||
def chl_a(wavelengths, reflectance):
|
||
"""计算叶绿素a"""
|
||
r651 = reflectance[find_nearest_bandnumber(651, wavelengths)]
|
||
r707 = reflectance[find_nearest_bandnumber(707, wavelengths)]
|
||
r670 = reflectance[find_nearest_bandnumber(670, wavelengths)]
|
||
|
||
denominator = safe_divide(r707 - r670, 1, default_value=1)
|
||
content = 5.534 + 9.809 * (r651 - r707) / denominator
|
||
return replace_non_valid_number(content)
|
||
|
||
|
||
def pc(wavelengths, reflectance):
|
||
"""计算藻蓝素"""
|
||
r620 = reflectance[find_nearest_bandnumber(620, wavelengths)]
|
||
r665 = reflectance[find_nearest_bandnumber(665, wavelengths)]
|
||
r709 = reflectance[find_nearest_bandnumber(709, wavelengths)]
|
||
r778 = reflectance[find_nearest_bandnumber(778, wavelengths)] # 修正了波段索引
|
||
|
||
# 计算背景值bb
|
||
bb_denom = safe_divide(0.082 - 0.6 * r778, 1, default_value=1)
|
||
bb = 2.71 * 0.6 * r778 / bb_denom
|
||
|
||
# 计算叶绿素chl
|
||
chl_denom = safe_divide(r665, 1, default_value=1)
|
||
chl = (r709 / chl_denom * (0.727 + bb) - bb - 0.401) * (1 / 0.68)
|
||
|
||
# 计算藻蓝素
|
||
pc_denom = safe_divide(r620, 1, default_value=1)
|
||
content = (r709 / pc_denom * (0.727 + bb) - bb - 0.281) * (1 / 0.84) - 0.24 * chl
|
||
return replace_non_valid_number(content)
|
||
|
||
|
||
def nh3(wavelengths, reflectance):
|
||
"""计算氨氮"""
|
||
band_500 = reflectance[find_nearest_bandnumber(500, wavelengths)]
|
||
band_600 = reflectance[find_nearest_bandnumber(600, wavelengths)]
|
||
band_850 = reflectance[find_nearest_bandnumber(850, wavelengths)]
|
||
|
||
ratio1 = safe_divide(band_500, band_850, default_value=1)
|
||
ratio2 = safe_divide(band_600, band_500, default_value=1)
|
||
|
||
x13 = np.log(ratio1) if ratio1 > 0 else 0
|
||
x23 = np.exp(ratio2)
|
||
|
||
content = 0.139 * x23 - 0.508 * x13 - 0.672
|
||
return replace_non_valid_number(content)
|
||
|
||
|
||
def tn(wavelengths, reflectance):
|
||
"""计算总氮"""
|
||
band_500 = reflectance[find_nearest_bandnumber(500, wavelengths)]
|
||
band_600 = reflectance[find_nearest_bandnumber(600, wavelengths)]
|
||
band_850 = reflectance[find_nearest_bandnumber(850, wavelengths)]
|
||
|
||
ratio1 = safe_divide(band_500, band_850, default_value=1)
|
||
ratio2 = safe_divide(band_600, band_500, default_value=1)
|
||
|
||
x13 = np.log(ratio1) if ratio1 > 0 else 0
|
||
x23 = np.exp(ratio2)
|
||
|
||
content = 0.059 * x23 - 0.089 * x13 - 0.012
|
||
return replace_non_valid_number(content)
|
||
|
||
|
||
def tp(wavelengths, reflectance):
|
||
"""计算总磷"""
|
||
band_440 = reflectance[find_nearest_bandnumber(440, wavelengths)]
|
||
band_500 = reflectance[find_nearest_bandnumber(500, wavelengths)]
|
||
band_800 = reflectance[find_nearest_bandnumber(800, wavelengths)]
|
||
|
||
ratio1 = safe_divide(band_500, band_440, default_value=1)
|
||
ratio2 = safe_divide(band_500, band_800, default_value=1)
|
||
|
||
x12 = np.log(ratio1) if ratio1 > 0 else 0
|
||
x17 = np.exp(ratio2)
|
||
|
||
content = 0.066 * x12 - 0.159 * x17 - 0.012
|
||
return replace_non_valid_number(content)
|
||
|
||
|
||
def mno4(wavelengths, reflectance):
|
||
"""计算高锰酸盐"""
|
||
band_440 = reflectance[find_nearest_bandnumber(440, wavelengths)]
|
||
band_500 = reflectance[find_nearest_bandnumber(500, wavelengths)]
|
||
band_610 = reflectance[find_nearest_bandnumber(610, wavelengths)]
|
||
band_800 = reflectance[find_nearest_bandnumber(800, wavelengths)]
|
||
|
||
x3 = safe_divide(band_500, band_440, default_value=1)
|
||
x6 = safe_divide(band_610, band_800, default_value=1)
|
||
|
||
content = 1.13 * x3 - 1.29 * x6 - 5.95
|
||
return replace_non_valid_number(content)
|
||
|
||
|
||
def tss(wavelengths, reflectance):
|
||
"""计算总悬浮物"""
|
||
band_490 = reflectance[find_nearest_bandnumber(490, wavelengths)]
|
||
band_555 = reflectance[find_nearest_bandnumber(555, wavelengths)]
|
||
band_670 = reflectance[find_nearest_bandnumber(670, wavelengths)]
|
||
|
||
x1 = band_555 + band_670
|
||
x2 = safe_divide(band_490, band_555, default_value=1)
|
||
|
||
content_tmp = 0.58213 + 23.84071 * x1 - 0.48532 * x2
|
||
content = np.exp(content_tmp)
|
||
|
||
return replace_non_valid_number(content)
|
||
|
||
|
||
def CDOM(wavelengths, reflectance):
|
||
"""计算有色可溶性有机物"""
|
||
band_560 = reflectance[find_nearest_bandnumber(560, wavelengths)]
|
||
band_660 = reflectance[find_nearest_bandnumber(660, wavelengths)]
|
||
|
||
x1 = band_560 + band_660
|
||
content = 40.75 * np.exp(-2.463 * x1)
|
||
|
||
return replace_non_valid_number(content)
|
||
|
||
|
||
def SD(wavelengths, reflectance):
|
||
"""计算透明度"""
|
||
band_440 = reflectance[find_nearest_bandnumber(440, wavelengths)]
|
||
band_700 = reflectance[find_nearest_bandnumber(700, wavelengths)]
|
||
|
||
x1 = safe_divide(band_700, band_440, default_value=1)
|
||
content = 1.17 * x1 - 27.43
|
||
|
||
return replace_non_valid_number(content)
|
||
|
||
|
||
def process_spectral_file(input_csv_path, output_csv_path):
|
||
# 读取CSV文件，第一行为列名
|
||
df = pd.read_csv(input_csv_path)
|
||
|
||
# 提取波长数据（假设波长在列名中）
|
||
# 跳过前18列非光谱数据，从第19列开始是光谱数据
|
||
wavelength_columns = df.columns[1:]
|
||
|
||
# 将列名转换为浮点数（波长值）
|
||
try:
|
||
wavelengths = np.array([float(col) for col in wavelength_columns])
|
||
except ValueError:
|
||
# 如果列名不是数值，尝试从第一行获取波长值
|
||
wavelengths = df.iloc[0, 1:].astype(float).values
|
||
# 删除第一行（波长行）
|
||
df = df.iloc[1:]
|
||
|
||
# 创建结果DataFrame
|
||
results = {
|
||
'timestamp': df.iloc[:, 0], # 假设时间戳在第一列
|
||
'chl-a': [],
|
||
'pc': [],
|
||
'nh3': [],
|
||
'tn': [],
|
||
'tp': [],
|
||
'mno4': [],
|
||
'tss': [],
|
||
'CDOM': [],
|
||
'SD': []
|
||
}
|
||
|
||
# 逐行处理数据
|
||
for i, row in df.iterrows():
|
||
# 提取反射率数据（第19列到最后一列）
|
||
reflectance = row.iloc[1:].astype(float).values
|
||
|
||
# 去除水表耀光
|
||
|
||
deglinted_ref = reflectance
|
||
# deglinted_ref = deglint(wavelengths, reflectance)
|
||
|
||
# 计算各项水质参数
|
||
results['chl-a'].append(chl_a(wavelengths, deglinted_ref))
|
||
results['pc'].append(pc(wavelengths, deglinted_ref))
|
||
results['nh3'].append(nh3(wavelengths, deglinted_ref))
|
||
results['tn'].append(tn(wavelengths, deglinted_ref))
|
||
results['tp'].append(tp(wavelengths, deglinted_ref))
|
||
results['mno4'].append(mno4(wavelengths, deglinted_ref))
|
||
results['tss'].append(tss(wavelengths, deglinted_ref))
|
||
results['CDOM'].append(CDOM(wavelengths, deglinted_ref))
|
||
results['SD'].append(SD(wavelengths, deglinted_ref))
|
||
|
||
# 创建结果DataFrame并保存
|
||
result_df = pd.DataFrame(results)
|
||
result_df.to_csv(output_csv_path, index=False)
|
||
print(f"结果已保存到 {output_csv_path}")
|
||
#
|
||
|
||
input_csv_path = r"D:\WQ\zhanghuilai\hyperspectral-inversion\data\input\一代高光谱\2025-07-11_2025-07-21\筛选.csv"
|
||
output_csv_path = r"D:\WQ\zhanghuilai\hyperspectral-inversion\data\input\一代高光谱\2025-07-11_2025-07-21/不去耀斑index.csv"
|
||
process_spectral_file(input_csv_path, output_csv_path)
|
||
|
||
|
||
#
|
||
# if __name__ == "__main__":
|
||
# parser = argparse.ArgumentParser(description='高光谱水质参数反演')
|
||
# parser.add_argument('--input', type=str, required=True, help='输入CSV文件路径')
|
||
# parser.add_argument('--output', type=str, required=True, help='输出CSV文件路径')
|
||
# args = parser.parse_args()
|
||
#
|
||
# process_spectral_file(args.input, args.output) |