diff --git a/chla-model.py b/chla-model.py new file mode 100644 index 0000000..0e1a914 --- /dev/null +++ b/chla-model.py @@ -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) \ No newline at end of file