zhanghuilai/hyperspectral-hongshengreservoir
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

上传文件至 /

Browse Source
This commit is contained in:
zhanghuilai
2025-07-22 10:28:35 +08:00
parent bd82850ce4
commit fe1bd426fa
1 changed files with 239 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

239
chla-model.py Normal file
View File

@ -0,0 +1,239 @@
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)