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)