diff --git a/.idea/vcs.xml b/.idea/vcs.xml new file mode 100644 index 0000000..35eb1dd --- /dev/null +++ b/.idea/vcs.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/Flexbrdf/Flex_brdf.py b/Flexbrdf/Flex_brdf.py index ea23e92..d994b26 100644 --- a/Flexbrdf/Flex_brdf.py +++ b/Flexbrdf/Flex_brdf.py @@ -16,265 +16,57 @@ from Flexbrdf.hytools.masks import mask_create warnings.filterwarnings("ignore") np.seterr(divide='ignore', invalid='ignore') -# Default configuration values -DEFAULT_CONFIG = { - 'file_type': 'envi', - 'num_cpus': 16, - 'bad_bands': [], - 'corrections': [], - 'resample': False, - 'resampler': {'type': 'cubic', 'out_waves': []}, - 'export': { - 'coeffs':True, - 'image': True, - 'masks': False, - 'subset_waves': [], - 'output_dir': './output/', - 'suffix': 'BRDF' - }, - 'brdf': { - 'type': 'flex', - 'grouped': True, - 'geometric': 'li_dense_r', - 'volume': 'ross_thick', - 'b/r': 2.5, - 'h/b': 2.0, - 'sample_perc': 0.1, - 'interp_kind': 'linear', - 'calc_mask': [ - ['water', {'band_1': 850, 'band_2': 660,"threshold": 290}], - ['kernel_finite', {}], - ['ancillary', {'name': 'sensor_zn', 'min': 0.03490658503988659, 'max': 'inf'}] - ], - 'apply_mask': [ - ['water', {'band_1': 850, 'band_2': 660,"threshold": 290}] - ], - 'bin_type': 'dynamic', - 'num_bins': 18, - 'ndvi_bin_min': 0.05, - 'ndvi_bin_max': 1.0, - 'ndvi_perc_min': 10, - 'ndvi_perc_max': 95, - 'solar_zn_type': 'scene' - }, - # 'topo': { - # 'type': 'scs+c', - # 'calc_mask': [ - # ['ndi', {'band_1': 850, 'band_2': 660, 'min': 0.05, 'max': 1.0}], - # ['kernel_finite', {}], - # ['ancillary', {'name': 'sensor_zn', 'min': 0.03490658503988659, 'max': 'inf'}] - # ], - # 'apply_mask': [ - # ['ndi', {'band_1': 850, 'band_2': 660, 'min': 0.05, 'max': 1.0}] - # ], - # 'sample_perc': 0.1, - # 'subgrouped': False, - # 'subgroup': {} - # }, - # 'glint': { - # 'type': 'hochberg', - # 'correction_band': 560, - # 'deep_water_sample': {}, - # 'calc_mask': [ - # ['ndi', {'band_1': 550, 'band_2': 2150, 'min': -1, 'max': 0}], - # ['kernel_finite', {}] - # ], - # 'apply_mask': [ - # ['ndi', {'band_1': 550, 'band_2': 2150, 'min': -1, 'max': 0}] - # ] - # } -} - -def build_config_from_args(args): - """从命令行参数构建配置字典""" - # 自动发现输入文件 - input_files = [] - if os.path.isdir(args.input_dir): - # 支持的文件扩展名 - extensions = ['.tif', '.tiff', '.envi', '.img', '.dat'] - for file in os.listdir(args.input_dir): - if any(file.lower().endswith(ext) for ext in extensions): - input_files.append(os.path.join(args.input_dir, file)) - else: - input_files = [args.input_dir] if os.path.isfile(args.input_dir) else [] - - if not input_files: - raise ValueError(f"No input files found in {args.input_dir}") - - # 自动生成anc_files - anc_files = {} - for input_file in input_files: - base_name = os.path.splitext(os.path.basename(input_file))[0] - # 根据文件名模式生成对应的angles文件路径 - # 例如: 2025_9_2_3_53_45_202592_35252_0_rad_geo_corrected_reflectance.dat - # 对应的angles文件: 2025_9_2_3_53_45_202592_35252_0_rad_rgbxyz_geo_registered_angles.bip - parts = base_name.split('_') - if len(parts) >= 9: # 确保有足够的parts - anc_base = f"{parts[0]}_{parts[1]}_{parts[2]}_{parts[3]}_{parts[4]}_{parts[5]}_{parts[6]}_{parts[7]}_{parts[8]}_rad_rgbxyz_geo_angles_registered" - # 确保路径分隔符统一使用反斜杠（Windows） - anc_dir_clean = args.anc_dir.replace('/', '\\') - anc_path = os.path.join(anc_dir_clean, anc_base + ".bip") - # 确保最终路径格式正确 - anc_path = os.path.normpath(anc_path) - - if os.path.exists(anc_path): - # 生成完整的anc_files字典结构 - anc_files[input_file] = { - "path_length": [anc_path, 0], - "sensor_az": [anc_path, 9], - "sensor_zn": [anc_path, 8], - "solar_az": [anc_path, 7], - "solar_zn": [anc_path, 6], - "phase": [anc_path, 0], - "slope": [anc_path, 0], - "aspect": [anc_path, 0], - "cosine_i": [anc_path, 0], - "utc_time": [anc_path, 0] - } - else: - # 如果找不到对应的angles文件，使用默认的obs文件 - obs_path = os.path.join(args.anc_dir, base_name + "_obs") - if os.path.exists(obs_path): - anc_files[input_file] = obs_path - - # 构建配置字典 - config_dict = DEFAULT_CONFIG.copy() - config_dict.update({ - 'input_files': input_files, - 'anc_files': anc_files, - 'file_type': 'envi', # 假设都是ENVI格式 - 'num_cpus': args.num_cpus, - 'bad_bands': args.bad_bands if args.bad_bands else [], - 'corrections': args.corrections if args.corrections else [], - 'export': { - 'coeffs': args.export_coeffs, - 'image': args.export_image, - 'masks': args.export_masks, - 'subset_waves': args.subset_waves if args.subset_waves else [], - 'output_dir': args.output_dir, - 'suffix': args.suffix - } - }) - - # 根据命令行参数更新BRDF配置 - if 'brdf' in args.corrections: - brdf_config = config_dict['brdf'].copy() - if hasattr(args, 'brdf_type') and args.brdf_type: - brdf_config['type'] = args.brdf_type - if hasattr(args, 'grouped') and args.grouped is not None: - brdf_config['grouped'] = args.grouped - if hasattr(args, 'geometric') and args.geometric: - brdf_config['geometric'] = args.geometric - if hasattr(args, 'volume') and args.volume: - brdf_config['volume'] = args.volume - if hasattr(args, 'num_bins') and args.num_bins: - brdf_config['num_bins'] = args.num_bins - config_dict['brdf'] = brdf_config - - # 根据命令行参数更新TOPO配置 - if 'topo' in args.corrections: - topo_config = config_dict['topo'].copy() - if hasattr(args, 'topo_type') and args.topo_type: - topo_config['type'] = args.topo_type - config_dict['topo'] = topo_config - - # 根据命令行参数更新Glint配置 - if 'glint' in args.corrections: - glint_config = config_dict['glint'].copy() - if hasattr(args, 'glint_type') and args.glint_type: - glint_config['type'] = args.glint_type - config_dict['glint'] = glint_config - - return config_dict def main(): - parser = argparse.ArgumentParser(description="High-resolution image correction with automatic configuration") + parser = argparse.ArgumentParser(description="FlexBRDF - 只接受JSON配置文件输入") + parser.add_argument('config', help='JSON配置文件路径 (例如: config.json)') + + args = parser.parse_args() + + # 检查配置文件是否存在 + if not os.path.exists(args.config): + print(f"错误: 配置文件不存在: {args.config}") + sys.exit(1) + + # 检查是否为JSON文件 + if not args.config.endswith('.json'): + print("错误: 只接受JSON格式的配置文件") + sys.exit(1) + + # 加载JSON配置 + try: + with open(args.config, 'r', encoding='utf-8') as f: + config_dict = json.load(f) + except json.JSONDecodeError as e: + print(f"错误: JSON配置文件格式错误: {e}") + sys.exit(1) + except Exception as e: + print(f"错误: 无法读取配置文件: {e}") + sys.exit(1) - # 必需参数 - parser.add_argument('input_dir', help='Input directory containing image files or single image file path') - parser.add_argument('anc_dir', help='Ancillary directory containing angle/obs files') - parser.add_argument('--output-dir', required=True, help='Output directory for corrected images') + # 验证必需的配置字段 + required_fields = ['input_files', 'anc_files', 'corrections', 'export'] + missing_fields = [f for f in required_fields if f not in config_dict] + if missing_fields: + print(f"错误: 配置文件缺少必需字段: {', '.join(missing_fields)}") + print("必需字段包括: input_files, anc_files, corrections, export") + sys.exit(1) - # 可选参数 - 基本设置 - parser.add_argument('--num-cpus', type=int, default=16, help='Number of CPUs to use (default: 4)') - parser.add_argument('--bad-bands', nargs='*', type=int, default=[], help='Bad band indices to exclude') - - # 校正类型 - parser.add_argument('--corrections', nargs='*', choices=['topo', 'brdf', 'glint'], - default=['brdf'], help='Correction types to apply (default: brdf)') - - # BRDF参数 - parser.add_argument('--brdf-type', choices=['universal', 'flex'], default='flex', - help='BRDF correction type (default: flex)') - parser.add_argument('--grouped', action='store_true', default=True, - help='Group images for BRDF correction (default: True)') - parser.add_argument('--no-grouped', action='store_false', dest='grouped', - help='Do not group images for BRDF correction') - parser.add_argument('--geometric', default='li_dense_r', - choices=['li_sparse', 'li_dense', 'li_dense_r', 'roujean'], - help='Geometric kernel type (default: li_dense_r)') - parser.add_argument('--volume', default='ross_thick', - choices=['ross_thin', 'ross_thick', 'hotspot', 'roujean'], - help='Volume kernel type (default: ross_thick)') - parser.add_argument('--num-bins', type=int, default=18, - help='Number of NDVI bins for FlexBRDF (default: 18)') - - # TOPO参数 - parser.add_argument('--topo-type', default='scs+c', - choices=['scs', 'scs+c', 'c', 'cosine', 'mod_minneart'], - help='TOPO correction type (default: scs+c)') - - # Glint参数 - parser.add_argument('--glint-type', default='hochberg', - choices=['hochberg', 'gao', 'hedley'], - help='Glint correction type (default: hochberg)') - - # 输出选项 - parser.add_argument('--export-coeffs', action='store_true', default=True, - help='Export correction coefficients') - parser.add_argument('--export-image', action='store_true', default=True, - help='Export corrected images (default: True)') - parser.add_argument('--export-masks', action='store_true', default=True, - help='Export correction masks (default: True)') - parser.add_argument('--subset-waves', nargs='*', type=float, default=[], - help='Subset of wavelengths to export (empty for all)') - parser.add_argument('--suffix', default='corrected', - help='Suffix for output files (default: corrected)') - parser.add_argument('--output-config', type=str, default=None, - help='Output current configuration to JSON file') - - # 向后兼容：如果只有一个参数且是JSON文件，则使用传统模式 - if len(sys.argv) == 2 and sys.argv[1].endswith('.json'): - config_file = sys.argv[1] - with open(config_file, 'r') as outfile: - config_dict = json.load(outfile) - else: - args = parser.parse_args() - config_dict = build_config_from_args(args) - - # 如果指定了输出配置选项，则输出配置到JSON文件 - if hasattr(args, 'output_config') and args.output_config: - print(f"Outputting configuration to {args.output_config}...") - # 重新排序配置字典以匹配期望的格式 - ordered_config = { - "bad_bands": config_dict.get("bad_bands", []), - "file_type": config_dict.get("file_type", "envi"), - "input_files": config_dict.get("input_files", []), - "anc_files": config_dict.get("anc_files", {}), - "num_cpus": config_dict.get("num_cpus", 4), - "export": config_dict.get("export", {}), - "corrections": config_dict.get("corrections", []), - "brdf": config_dict.get("brdf", {}), - "topo": config_dict.get("topo", {}), - "glint": config_dict.get("glint", {}), - "resample": config_dict.get("resample", False), - "resampler": config_dict.get("resampler", {"type": "cubic", "out_waves": []}) - } - with open(args.output_config, 'w', encoding='utf-8') as f: - json.dump(ordered_config, f, indent=2, ensure_ascii=False) - print("Configuration output complete.") - return # 如果只是输出配置，则退出程序 + # 设置默认值 + config_dict.setdefault('file_type', 'envi') + config_dict.setdefault('num_cpus', 16) + config_dict.setdefault('bad_bands', []) + config_dict.setdefault('resample', False) + config_dict.setdefault('resampler', {'type': 'cubic', 'out_waves': []}) + print(f"=" * 60) + print(f"FlexBRDF 图像校正工具") + print(f"=" * 60) + print(f"配置文件: {args.config}") + print(f"输入文件数量: {len(config_dict['input_files'])}") + print(f"校正类型: {', '.join(config_dict['corrections'])}") + print(f"输出目录: {config_dict['export']['output_dir']}") + print(f"=" * 60) print("Starting image correction process...") total_start = time.time() diff --git a/Flexbrdf/examples/configs/1.json b/Flexbrdf/examples/configs/1.json new file mode 100644 index 0000000..1c1387a --- /dev/null +++ b/Flexbrdf/examples/configs/1.json @@ -0,0 +1,167 @@ +{ + "input_files": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/rad2geo/R/2025_9_2_3_53_45_202592_35252_0_rad_geo_corrected_reflectance.dat", + "D:/BaiduNetdiskDownload/20250902/_3_52_52/rad2geo/R/2025_9_2_3_53_45_202592_35252_1_rad_geo_corrected_reflectance.dat" + ], + "file_type": "envi", + "bad_bands": [], + "num_cpus": 2, + "anc_files": { + "D:/BaiduNetdiskDownload/20250902/_3_52_52/rad2geo/R/2025_9_2_3_53_45_202592_35252_0_rad_geo_corrected_reflectance.dat": { + "path_length": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_0_rad_rgbxyz_geo.bip_with_angles.bip", + 0 + ], + "sensor_az": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_0_rad_rgbxyz_geo.bip_with_angles.bip", + 9 + ], + "sensor_zn": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_0_rad_rgbxyz_geo.bip_with_angles.bip", + 8 + ], + "solar_az": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_0_rad_rgbxyz_geo.bip_with_angles.bip", + 7 + ], + "solar_zn": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_0_rad_rgbxyz_geo.bip_with_angles.bip", + 6 + ], + "phase": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_0_rad_rgbxyz_geo.bip_with_angles.bip", + 0 + ], + "slope": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_0_rad_rgbxyz_geo.bip_with_angles.bip", + 0 + ], + "aspect": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_0_rad_rgbxyz_geo.bip_with_angles.bip", + 0 + ], + "cosine_i": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_0_rad_rgbxyz_geo.bip_with_angles.bip", + 0 + ], + "utc_time": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_0_rad_rgbxyz_geo.bip_with_angles.bip", + 0 + ] + }, + "D:/BaiduNetdiskDownload/20250902/_3_52_52/rad2geo/R/2025_9_2_3_53_45_202592_35252_1_rad_geo_corrected_reflectance.dat": { + "path_length": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_1_rad_rgbxyz_geo.bip_angles.bip", + 0 + ], + "sensor_az": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_1_rad_rgbxyz_geo.bip_angles.bip", + 9 + ], + "sensor_zn": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_1_rad_rgbxyz_geo.bip_angles.bip", + 8 + ], + "solar_az": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_1_rad_rgbxyz_geo.bip_angles.bip", + 7 + ], + "solar_zn": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_1_rad_rgbxyz_geo.bip_angles.bip", + 6 + ], + "phase": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_1_rad_rgbxyz_geo.bip_angles.bip", + 0 + ], + "slope": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_1_rad_rgbxyz_geo.bip_angles.bip", + 0 + ], + "aspect": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_1_rad_rgbxyz_geo.bip_angles.bip", + 0 + ], + "cosine_i": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_1_rad_rgbxyz_geo.bip_angles.bip", + 0 + ], + "utc_time": [ + "D:/BaiduNetdiskDownload/20250902/_3_52_52/BRDF/angle/test/2025_9_2_3_53_45_202592_35252_1_rad_rgbxyz_geo.bip_angles.bip", + 0 + ] + } + }, + "corrections": [ + "brdf" + ], + "topo": { + "type": "scs+c", + "calc_mask": [], + "apply_mask": [], + "c_fit_type": "nnls" + }, + "brdf": { + "type": "flex", + "grouped": true, + "geometric": "li_dense_r", + "volume": "ross_thick", + "b/r": 2.5, + "h/b": 2.0, + "sample_perc": 0.1, + "interp_kind": "linear", + "calc_mask": [ + [ + "water", + { + "band_1": 850, + "band_2": 660, + "threshold": 290.0 + } + ], + [ + "kernel_finite", + {} + ], + [ + "ancillary", + { + "name": "sensor_zn", + "min": 0.034907, + "max": 1000000.0 + } + ] + ], + "apply_mask": [ + [ + "water", + { + "band_1": 850, + "band_2": 660, + "threshold": 290.0 + } + ] + ], + "bin_type": "dynamic", + "num_bins": 18, + "ndvi_bin_min": 0.05, + "ndvi_bin_max": 1.0, + "ndvi_perc_min": 10, + "ndvi_perc_max": 95, + "solar_zn_type": "scene" + }, + "glint": { + "type": "hochberg", + "correction_wave": 2150, + "apply_mask": [] + }, + "export": { + "coeffs": true, + "image": true, + "masks": true, + "subset_waves": [], + "output_dir": "E:/code/hytools-master/hytools-master/data/output", + "suffix": "brdf_glint" + }, + "resample": false +} \ No newline at end of file diff --git a/Flexbrdf/examples/img/config_json_gen_gui.jpg b/Flexbrdf/examples/img/config_json_gen_gui.jpg deleted file mode 100644 index dcac7b7..0000000 Binary files a/Flexbrdf/examples/img/config_json_gen_gui.jpg and /dev/null differ diff --git a/Flexbrdf/examples/img/emit_001_20231101T024133_2330502_014_rgb.jpg b/Flexbrdf/examples/img/emit_001_20231101T024133_2330502_014_rgb.jpg deleted file mode 100644 index c2e7da4..0000000 Binary files a/Flexbrdf/examples/img/emit_001_20231101T024133_2330502_014_rgb.jpg and /dev/null differ diff --git a/Flexbrdf/examples/img/emit_001_20231101T024133_2330502_014_rgb_warp.jpg b/Flexbrdf/examples/img/emit_001_20231101T024133_2330502_014_rgb_warp.jpg deleted file mode 100644 index f25dfc4..0000000 Binary files a/Flexbrdf/examples/img/emit_001_20231101T024133_2330502_014_rgb_warp.jpg and /dev/null differ diff --git a/Flexbrdf/examples/img/emit_aus_subset1_topo.png b/Flexbrdf/examples/img/emit_aus_subset1_topo.png deleted file mode 100644 index 11ec119..0000000 Binary files a/Flexbrdf/examples/img/emit_aus_subset1_topo.png and /dev/null differ diff --git a/Flexbrdf/examples/img/emit_aus_subset1_uncorr.png b/Flexbrdf/examples/img/emit_aus_subset1_uncorr.png deleted file mode 100644 index bc62dca..0000000 Binary files a/Flexbrdf/examples/img/emit_aus_subset1_uncorr.png and /dev/null differ diff --git a/Flexbrdf/hytools/topo/原理.md b/Flexbrdf/hytools/topo/原理.md new file mode 100644 index 0000000..60c55ec --- /dev/null +++ b/Flexbrdf/hytools/topo/原理.md @@ -0,0 +1,186 @@ +基于 Flexbrdf/hytools/topo 目录下的代码，我来详细解释地形校正的原理。 + +地形校正原理详解 +地形校正（Topographic Correction）的目的是消除地形起伏（坡度、坡向）对高光谱图像反射率的影响。由于太阳光入射角度在不同坡度和坡向上的变化，同一地物在不同地形位置会呈现不同的辐射亮度，地形校正就是要将这种差异归一化到水平地表的条件。 + +一、核心概念：入射角余弦（cosine i） +地形校正的核心是计算入射角余弦（cosine of incidence angle），即太阳光线与地表法线之间夹角的余弦值。 + +计算公式（在 topo.py:53）： + +cosine_i = np.cos(solar_zn)*np.cos(slope) + np.sin(solar_zn)*np.sin(slope)*np.cos(relative_az) +其中： + +solar_zn：太阳天顶角（太阳光线与垂直方向的夹角） +solar_az：太阳方位角 +slope：地面坡度 +aspect：地面坡向 +relative_az = aspect - solar_az：相对方位角 +物理意义： + +当 cosine_i = 1 时，太阳光垂直照射地面（正午水平地表） +当 cosine_i < 1 时，太阳光斜射，照度降低 +背阴坡 cosine_i 较小，向阳坡 cosine_i 较大 +二、支持的5种地形校正方法 +1. Cosine 校正（余弦校正） +原理：最简单的地形校正方法，假设地表为朗伯体（各向同性反射）。 + +校正公式（cosine.py:50）： + +c_factor = cos_solar_zn / cos_i +cos_solar_zn：水平地表的太阳天顶角余弦 +cos_i：实际地表的入射角余弦 +校正效果： + +corrected_reflectance = original_reflectance × (cos_solar_zn / cos_i) +特点： + +简单快速 +假设地表是理想朗伯体，对非朗伯体（如植被）校正过度 +2. SCS 校正（Sun-Canopy-Sensor） +原理：考虑森林冠层对地形影响的简化模型，基于 "cosine i" 但增加坡度因子。 + +校正公式（scs.py:51）： + +c_factor = (cos_slope × cos_solar_zn) / cos_i +与 Cosine 方法的区别： + +增加了 cos_slope 项（坡度余弦） +考虑了坡度对观测路径的影响 +适用场景：森林覆盖区域 + +3. C 校正（C-Correction） +原理：在 Cosine 校正基础上增加一个经验系数 C，用于调整非朗伯体效应。 + +校正公式（c.py:177）： + +correction_factor = (cos_sz + C) / (cosine_i + C) +其中 C 系数的计算过程（c.py:42-76）： + +对每个波段，建立反射率与 cosine_i 的线性关系： + +ρ = slope × cos_i + intercept +通过最小二乘拟合（OLS）或非负最小二乘（NNLS）求解 + +计算 C 系数（Soenen et al. 2005, Eq. 8）： + +C = intercept / slope +物理意义：C 值越大，表示地表越接近朗伯体，校正强度越接近 Cosine 方法；C 值越小，校正越保守。 + +拟合方法选择（c_fit_type）： + +ols：普通最小二乘（Ordinary Least Squares） +nnls：非负最小二乘（Non-Negative Least Squares） +wls：加权最小二乘（Weighted Least Squares） +4. SCS+C 校正（推荐用于森林） +原理：结合 SCS 和 C 校正的优势，是目前森林地形校正的主流方法（Soenen et al. 2005, IEEE TGRS）。 + +计算步骤（scsc.py）： + +计算 c1 因子： + +c1 = cos(solar_zn) × cos(slope) +计算 C 系数（与 C 校正相同，通过线性回归） + +综合校正因子（scsc.py:115-125）： + +correction_factor = (c1 + C × cos_slope) / (cos_i + C) +为什么有效： + +SCS 部分处理冠层-地形几何关系 +C 部分调整非朗伯体效应 +相比单独使用 SCS 或 C 校正，能更好地保持植被反射率的一致性 +5. Modified Minnaert 校正（改进的 Minnaert） +原理：基于 Minnaert 指数的半经验方法，区分植被和非植被区域进行不同强度的校正（Richter et al. 2009）。 + +算法流程（modminn.py:46-83）： + +计算调整后的太阳天顶角： + +if solar_zn < 45°: solar_zn_t = solar_zn + 20° +if 45° ≤ solar_zn ≤ 55°: solar_zn_t = solar_zn + 15° +if solar_zn > 55°: solar_zn_t = solar_zn + 10° +NDVI 掩膜分类： + +植被区（NDVI > 0.2） +非植被区（NDVI ≤ 0.2） +计算基础校正因子： + +factor = cos_i / cos(solar_zn_t) +分区应用不同指数： + +区域 短波（<720nm） 长波（≥720nm） +非植被 +factor^(1/2) +factor^(1/2) +植被 +factor^(3/4) +factor^(1/3) +限幅处理： + +factor = np.clip(factor, 0.25, 1.0) +阈值控制：当入射角小于调整后的太阳天顶角时不进行校正 + +特点： + +自适应不同地物类型 +区分短波和长波有不同的散射特性 +通过指数控制校正强度，避免过度校正 +三、地形校正的工作流程 +从 topo.py:104-154 可以看出地形校正的完整流程： + +1. 设置地形参数（set_topo） + ↓ +2. 生成计算掩膜（calc_mask） + ↓ +3. 计算校正系数（calc_*_coeffs） + - 遍历每个波段（跳过坏波段） + - 根据选定的算法类型计算校正因子 + ↓ +4. 应用校正（apply_*） + - 支持多种数据切分方式：line/column/band/chunk/pixels +四、关键配置参数 +参数 类型 说明 适用算法 +type +str +校正算法类型 +所有 +calc_mask +list +计算校正系数的掩膜条件 +所有 +apply_mask +list +应用校正的掩膜条件 +所有 +c_fit_type +str +C系数拟合方法：ols/nnls/wls +c, scs+c +sample_perc +float +采样比例（用于分组计算） +分组计算时 +五、算法选择建议 +场景 推荐算法 理由 +裸土/岩石 +Cosine 或 C +地表接近朗伯体 +森林植被 +SCS+C +考虑冠层几何，校正效果最好 +草地/灌木 +SCS +简单有效 +混合地物 +Modified Minnaert +自适应不同地物 +需要精确结果 +C 或 SCS+C + 分组计算 +按NDVI分组计算C系数 +六、参考文献 +SCS+C 方法：Soenen, S. A., Peddle, D. R., & Coburn, C. A. (2005). SCS+C: A Modified Sun-Canopy-Sensor Topographic Correction in Forested Terrain. IEEE TGRS, 43(9), 2148-2159. + +方法比较：Richter, R., Kellenberger, T., & Kaufmann, H. (2009). Comparison of topographic correction methods. Remote Sensing, 1(3), 184-196. + +C 校正：Teillet, P. M., Guindon, B., & Goodenough, D. G. (1982). On the slope-aspect correction of multispectral scanner data. Canadian Journal of Remote Sensing, 8(2), 84-106. \ No newline at end of file diff --git a/GUI/brdf_gui.py b/GUI/brdf_gui.py index 6c05fc0..e72260f 100644 --- a/GUI/brdf_gui.py +++ b/GUI/brdf_gui.py @@ -1,17 +1,43 @@ # -*- coding: utf-8 -*- +""" +BRDF校正工具 - 统一GUI应用 +支持陆地植被BRDF校正（FlexBRDF）和水体BRDF校正（ocbrdf）两个模块 +""" import sys import os import json +import subprocess +import threading from pathlib import Path from PyQt5.QtWidgets import ( QApplication, QMainWindow, QWidget, QVBoxLayout, QHBoxLayout, QGroupBox, QFormLayout, QLabel, QLineEdit, QPushButton, QFileDialog, QComboBox, - QCheckBox, QSpinBox, QDoubleSpinBox, QTabWidget, QMessageBox, QScrollArea, + QCheckBox, QRadioButton, QSpinBox, QDoubleSpinBox, QTabWidget, QMessageBox, QScrollArea, QStackedWidget, QGridLayout, QDialog, QDialogButtonBox, - QListWidget, QListWidgetItem # 添加这两个 + QListWidget, QListWidgetItem, QTextEdit, QProgressBar, QSplitter ) -from PyQt5.QtCore import Qt +from PyQt5.QtCore import Qt, QThread, pyqtSignal, QTimer, QItemSelectionModel +from PyQt5.QtGui import QFont, QIcon, QPalette, QWheelEvent + + +# ==================== 禁用滚轮的自定义控件 ==================== +class NoWheelSpinBox(QSpinBox): + """禁用鼠标滚轮的SpinBox""" + def wheelEvent(self, event: QWheelEvent): + event.ignore() + + +class NoWheelDoubleSpinBox(QDoubleSpinBox): + """禁用鼠标滚轮的DoubleSpinBox""" + def wheelEvent(self, event: QWheelEvent): + event.ignore() + + +class NoWheelComboBox(QComboBox): + """禁用鼠标滚轮的ComboBox""" + def wheelEvent(self, event: QWheelEvent): + event.ignore() # ==================== 掩膜条件编辑器 ==================== @@ -80,7 +106,7 @@ class MaskConditionEditor(QWidget): # 类型选择行 type_layout = QHBoxLayout() - self.type_combo = QComboBox() + self.type_combo = NoWheelComboBox() self.type_combo.addItems(list(self.MASK_TYPES.keys())) self.type_combo.currentTextChanged.connect(self.on_type_changed) type_layout.addWidget(QLabel("类型:")) @@ -110,33 +136,36 @@ class MaskConditionEditor(QWidget): widget = QWidget() layout = QFormLayout(widget) layout.setContentsMargins(10, 5, 10, 5) + layout.setFieldGrowthPolicy(QFormLayout.ExpandingFieldsGrow) + layout.setLabelAlignment(Qt.AlignLeft) + layout.setFormAlignment(Qt.AlignLeft) for param_name, param_type, default_value in params: if param_type == int: if "band" in param_name or "radius" in param_name: - spin = QSpinBox() + spin = NoWheelSpinBox() spin.setRange(0, 3000) spin.setValue(default_value) param_widget = spin else: - spin = QSpinBox() + spin = NoWheelSpinBox() spin.setRange(-10000, 10000) spin.setValue(default_value) param_widget = spin elif param_type == float: - dspin = QDoubleSpinBox() + dspin = NoWheelDoubleSpinBox() dspin.setRange(-1e6, 1e6) dspin.setDecimals(6) dspin.setValue(default_value) param_widget = dspin else: # str if param_name == "name" and "name_choices" in self.MASK_TYPES[mask_type]: - combo = QComboBox() + combo = NoWheelComboBox() combo.addItems(self.MASK_TYPES[mask_type]["name_choices"]) combo.setCurrentText(default_value) param_widget = combo elif param_name == "method" and "method_choices" in self.MASK_TYPES[mask_type]: - combo = QComboBox() + combo = NoWheelComboBox() combo.addItems(self.MASK_TYPES[mask_type]["method_choices"]) combo.setCurrentText(default_value) param_widget = combo @@ -281,19 +310,331 @@ class MaskListWidget(QWidget): self.add_condition(cond) -# ==================== 主窗口 ==================== +# ==================== 命令执行线程 ==================== +class CommandThread(QThread): + """后台执行命令行程序的线程""" + output_updated = pyqtSignal(str) # 输出更新信号 + progress_updated = pyqtSignal(int) # 进度更新信号 + finished_signal = pyqtSignal(bool, str) # 完成信号(成功/失败, 消息) + + def __init__(self, command, working_dir=None, env=None): + super().__init__() + self.command = command + self.working_dir = working_dir + self.env = env # 环境变量字典 + self.process = None + + def run(self): + try: + # 合并环境变量 + run_env = os.environ.copy() + if self.env: + run_env.update(self.env) + + # 设置Python环境变量以确保UTF-8编码 + run_env['PYTHONIOENCODING'] = 'utf-8' + + self.process = subprocess.Popen( + self.command, + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, + bufsize=1, + cwd=self.working_dir, + shell=True, + env=run_env, + encoding='utf-8', + errors='replace' # 无法解码的字符用替换 + ) + + while True: + output = self.process.stdout.readline() + if output == '' and self.process.poll() is not None: + break + if output: + self.output_updated.emit(output.strip()) + + return_code = self.process.poll() + if return_code == 0: + self.finished_signal.emit(True, "处理完成！") + else: + self.finished_signal.emit(False, f"处理失败，退出代码：{return_code}") + + except Exception as e: + self.finished_signal.emit(False, f"执行错误：{str(e)}") + + def stop(self): + if self.process: + self.process.terminate() + + +# ==================== 模块选择主窗口 ==================== class MainWindow(QMainWindow): def __init__(self): super().__init__() - self.setWindowTitle("BRDF / Topo / Glint 校正配置工具") - self.setMinimumSize(1100, 800) + self.setWindowTitle("BRDF 校正工具 - 高光谱图像处理套件") + self.setMinimumSize(1000, 700) + + # 设置应用图标和样式 + self.setup_styles() + + # 创建主界面 + self.setup_main_interface() + + def setup_styles(self): + """设置NVIDIA风格的样式""" + self.setStyleSheet(""" + QMainWindow { + background-color: #000000; + } + QWidget { + background-color: #000000; + color: #ffffff; + font-family: 'Segoe UI', Arial, sans-serif; + font-size: 14px; + } + QLabel { + color: #ffffff; + font-size: 14px; + } + QPushButton { + background-color: transparent; + color: #ffffff; + border: 2px solid #76b900; + border-radius: 8px; + padding: 12px 24px; + font-weight: bold; + font-size: 16px; + min-height: 40px; + } + QPushButton:hover { + background-color: #76b900; + color: #000000; + } + QPushButton:pressed { + background-color: #5e9400; + border-color: #5e9400; + } + QGroupBox { + background-color: #1a1a1a; + border: 2px solid #76b900; + border-radius: 8px; + margin-top: 20px; + padding-top: 20px; + font-weight: bold; + font-size: 16px; + } + QGroupBox::title { + subcontrol-origin: margin; + left: 15px; + padding: 5px 15px; + color: #000000; + background-color: #76b900; + border-radius: 4px; + font-weight: bold; + } + """) + + def setup_main_interface(self): + """创建模块选择主界面""" + central_widget = QWidget() + self.setCentralWidget(central_widget) + + main_layout = QVBoxLayout(central_widget) + main_layout.setContentsMargins(40, 40, 40, 40) + main_layout.setSpacing(30) + + # 标题区域 + title_label = QLabel("BRDF 校正工具") + title_label.setAlignment(Qt.AlignCenter) + title_label.setStyleSheet(""" + font-size: 36px; + font-weight: bold; + color: #76b900; + margin-bottom: 20px; + """) + main_layout.addWidget(title_label) + + subtitle_label = QLabel("高光谱图像双向反射分布函数校正 - 支持陆地植被和水体两种场景") + subtitle_label.setAlignment(Qt.AlignCenter) + subtitle_label.setStyleSheet(""" + font-size: 16px; + color: #cccccc; + margin-bottom: 40px; + """) + main_layout.addWidget(subtitle_label) + + # 模块选择区域 + modules_layout = QHBoxLayout() + modules_layout.setSpacing(40) + + # 陆地植被BRDF模块 + land_group = self.create_module_card( + title="陆地植被 BRDF 校正", + subtitle="FlexBRDF算法", + description="适用于陆地植被和地表的BRDF、地形和镜面反射校正\n支持AVIRIS、AVIRIS-NG等机载高光谱数据", + features=[ + "✓ FlexBRDF 算法", + "✓ 地形校正 (SCS, SCS+C, Cosine等)", + "✓ 镜面反射校正", + "✓ 多种掩膜条件", + "✓ 并行处理支持" + ], + color="#76b900", + button_text="进入陆地模块", + callback=self.open_land_module + ) + modules_layout.addWidget(land_group) + + # 水体BRDF模块 + water_group = self.create_module_card( + title="水体 BRDF 校正", + subtitle="Ocean BRDF Correction", + description="专用于水体和海洋的BRDF校正\n支持多种BRDF模型和不确定性估算", + features=[ + "✓ L11, M02, O25 BRDF模型", + "✓ 不确定性估算", + "✓ 水体掩膜自动识别", + "✓ 分块处理支持", + "✓ 多种输出格式" + ], + color="#0080ff", + button_text="进入水体模块", + callback=self.open_water_module + ) + modules_layout.addWidget(water_group) + + main_layout.addLayout(modules_layout) + main_layout.addStretch() + + def create_module_card(self, title, subtitle, description, features, color, button_text, callback): + """创建模块卡片""" + group = QGroupBox() + group.setMinimumHeight(400) + group.setStyleSheet(f""" + QGroupBox {{ + background-color: #1a1a1a; + border: 2px solid {color}; + border-radius: 12px; + margin: 10px; + padding: 20px; + }} + """) + + layout = QVBoxLayout(group) + layout.setSpacing(20) + + # 标题 + title_label = QLabel(title) + title_label.setAlignment(Qt.AlignCenter) + title_label.setStyleSheet(f""" + font-size: 24px; + font-weight: bold; + color: {color}; + margin-bottom: 5px; + """) + layout.addWidget(title_label) + + # 副标题 + subtitle_label = QLabel(subtitle) + subtitle_label.setAlignment(Qt.AlignCenter) + subtitle_label.setStyleSheet(""" + font-size: 14px; + color: #aaaaaa; + margin-bottom: 15px; + """) + layout.addWidget(subtitle_label) + + # 描述 + desc_label = QLabel(description) + desc_label.setAlignment(Qt.AlignCenter) + desc_label.setWordWrap(True) + desc_label.setStyleSheet(""" + font-size: 13px; + color: #dddddd; + margin-bottom: 20px; + line-height: 1.4; + """) + layout.addWidget(desc_label) + + # 功能特性列表 + features_label = QLabel("\n".join(features)) + features_label.setAlignment(Qt.AlignLeft) + features_label.setStyleSheet(""" + font-size: 12px; + color: #cccccc; + margin-bottom: 20px; + background-color: #2a2a2a; + padding: 15px; + border-radius: 6px; + """) + layout.addWidget(features_label) + + layout.addStretch() + + # 操作按钮 + button = QPushButton(button_text) + button.setStyleSheet(f""" + QPushButton {{ + background-color: transparent; + color: {color}; + border: 2px solid {color}; + border-radius: 8px; + padding: 12px 24px; + font-weight: bold; + font-size: 16px; + min-height: 40px; + }} + QPushButton:hover {{ + background-color: {color}; + color: #000000; + }} + QPushButton:pressed {{ + background-color: {color}; + opacity: 0.8; + }} + """) + button.clicked.connect(callback) + layout.addWidget(button) + + return group + + def open_land_module(self): + """打开陆地植被BRDF模块""" + self.land_window = LandBRDFWindow() + self.land_window.show() + # 可选：隐藏主窗口 + # self.hide() + + def open_water_module(self): + """打开水体BRDF模块""" + self.water_window = WaterBRDFWindow() + self.water_window.show() + # 可选：隐藏主窗口 + # self.hide() + + +# ==================== 陆地植被BRDF模块窗口 ==================== +class LandBRDFWindow(QMainWindow): + def __init__(self): + super().__init__() + self.setWindowTitle("陆地植被 BRDF 校正 - FlexBRDF") + self.setMinimumSize(1200, 800) + self.resize(1400, 900) # 设置默认大小 # 存储匹配后的文件对列表 self.file_pairs = [] # [(rfl_path, obs_path), ...] # 存储每个输入文件对应的辅助文件波段映射 self.band_mappings = {} # {rfl_path: {param_name: band_index}} - - # 设置全局 NVIDIA 风格样式表 + + # 命令执行相关 + self.command_thread = None + + # 设置样式和界面 + self.setup_land_styles() + self.setup_land_interface() + + def setup_land_styles(self): + """设置陆地植被BRDF模块的样式""" self.setStyleSheet(""" QMainWindow { background-color: #000000; @@ -396,17 +737,46 @@ class MainWindow(QMainWindow): border: none; background-color: #000000; } + QTextEdit { + background-color: #1a1a1a; + color: #ffffff; + border: 1px solid #5e5e5e; + border-radius: 2px; + font-family: 'Consolas', 'Monaco', monospace; + font-size: 12px; + } + QProgressBar { + background-color: #1a1a1a; + border: 1px solid #5e5e5e; + border-radius: 2px; + text-align: center; + color: #ffffff; + } + QProgressBar::chunk { + background-color: #76b900; + border-radius: 2px; + } """) + def setup_land_interface(self): + """创建陆地植被BRDF模块界面""" central_widget = QWidget() self.setCentralWidget(central_widget) main_layout = QVBoxLayout(central_widget) main_layout.setSpacing(15) main_layout.setContentsMargins(20, 20, 20, 20) + # 创建分割器：左侧配置，右侧日志 + splitter = QSplitter(Qt.Horizontal) + + # 左侧：配置选项卡 + left_widget = QWidget() + left_layout = QVBoxLayout(left_widget) + tabs = QTabWidget() tabs.setDocumentMode(True) - main_layout.addWidget(tabs) + # tabs.setMaximumHeight(500) # 限制标签页高度，确保按钮可见 + left_layout.addWidget(tabs, 1) # stretch factor = 1 # 1. 文件与波段映射页 self.file_tab = QWidget() @@ -423,21 +793,345 @@ class MainWindow(QMainWindow): tabs.addTab(self.output_tab, "输出设置") self.setup_output_tab() - # 底部按钮 - btn_layout = QHBoxLayout() - load_btn = QPushButton("加载 JSON 配置") + # 底部按钮区域 - 使用独立的分组框使其更显眼 + btn_group = QGroupBox("操作") + btn_group.setMinimumHeight(90) # 确保按钮组有足够高度 + btn_group.setStyleSheet(""" + QGroupBox { + background-color: #1a1a1a; + border: 2px solid #76b900; + border-radius: 4px; + margin-top: 5px; + padding: 10px; + } + QGroupBox::title { + subcontrol-origin: margin; + left: 10px; + padding: 2px 10px; + color: #76b900; + font-weight: bold; + font-size: 13px; + } + """) + btn_layout = QHBoxLayout(btn_group) + btn_layout.setSpacing(15) + + load_btn = QPushButton("📂 加载 JSON 配置") + load_btn.setMinimumHeight(40) load_btn.clicked.connect(self.load_config) - save_btn = QPushButton("保存 JSON 配置") + + save_btn = QPushButton("💾 保存 JSON 配置") + save_btn.setMinimumHeight(40) save_btn.clicked.connect(self.save_config) + + generate_run_btn = QPushButton("▶ 生成配置并运行") + generate_run_btn.setMinimumHeight(55) + generate_run_btn.setMinimumWidth(220) + generate_run_btn.setStyleSheet(""" + QPushButton { + background-color: #76b900; + color: #000000; + border: 3px solid #76b900; + border-radius: 8px; + font-weight: bold; + font-size: 16px; + padding: 15px 30px; + } + QPushButton:hover { + background-color: #8fce00; + border-color: #8fce00; + } + QPushButton:pressed { + background-color: #5e9400; + border-color: #5e9400; + } + """) + generate_run_btn.clicked.connect(self.generate_and_run) + btn_layout.addWidget(load_btn) btn_layout.addWidget(save_btn) - main_layout.addLayout(btn_layout) + btn_layout.addStretch() + btn_layout.addWidget(generate_run_btn) + left_layout.addWidget(btn_group, 0, Qt.AlignBottom) # 固定高度，底部对齐 + + # 右侧：日志和进度显示 + right_widget = QWidget() + right_layout = QVBoxLayout(right_widget) + + log_label = QLabel("执行日志") + log_label.setStyleSheet("font-weight: bold; color: #76b900; font-size: 16px;") + right_layout.addWidget(log_label) + + self.log_text = QTextEdit() + self.log_text.setMaximumWidth(400) + self.log_text.setPlainText("等待执行命令...\n") + right_layout.addWidget(self.log_text) + + self.progress_bar = QProgressBar() + self.progress_bar.setVisible(False) + right_layout.addWidget(self.progress_bar) + + # 添加到分割器 + splitter.addWidget(left_widget) + splitter.addWidget(right_widget) + splitter.setStretchFactor(0, 2) # 左侧占2/3 + splitter.setStretchFactor(1, 1) # 右侧占1/3 + + main_layout.addWidget(splitter) + + def generate_and_run(self): + """生成JSON配置并运行FlexBRDF""" + if not self.file_pairs: + QMessageBox.warning(self, "警告", "没有有效的文件对，请先匹配文件。") + return + + # 生成JSON配置 + config = self.build_config() + + # 保存临时JSON文件 + temp_json_path = os.path.join(os.path.dirname(__file__), "temp_flexbrdf_config.json") + try: + with open(temp_json_path, 'w', encoding='utf-8') as f: + json.dump(config, f, indent=2, ensure_ascii=False) + except Exception as e: + QMessageBox.critical(self, "错误", f"保存配置文件失败: {e}") + return + + # 构建命令 + flexbrdf_path = os.path.join(os.path.dirname(__file__), "..", "Flexbrdf", "Flex_brdf.py") + flexbrdf_path = os.path.abspath(flexbrdf_path) + flexbrdf_dir = os.path.dirname(flexbrdf_path) + + if not os.path.exists(flexbrdf_path): + QMessageBox.critical(self, "错误", f"找不到 Flex_brdf.py 文件:\n{flexbrdf_path}") + return + + # 使用当前Python解释器 + python_exe = sys.executable + command = f'"{python_exe}" "{flexbrdf_path}" "{temp_json_path}"' + + # 显示进度条和日志 + self.progress_bar.setVisible(True) + self.progress_bar.setRange(0, 0) # 不确定进度 + self.log_text.clear() + self.log_text.append("开始执行 FlexBRDF 校正...") + self.log_text.append(f"Python: {python_exe}") + self.log_text.append(f"命令: {command}\n") + + # 启动命令执行线程 + self.command_thread = CommandThread(command, flexbrdf_dir) + self.command_thread.output_updated.connect(self.append_log) + self.command_thread.finished_signal.connect(self.command_finished) + self.command_thread.start() + + def append_log(self, text): + """添加日志文本""" + self.log_text.append(text) + # 自动滚动到底部 + cursor = self.log_text.textCursor() + cursor.movePosition(cursor.End) + self.log_text.setTextCursor(cursor) + + def command_finished(self, success, message): + """命令执行完成""" + self.progress_bar.setVisible(False) + self.log_text.append(f"\n{message}") + + if success: + QMessageBox.information(self, "成功", "FlexBRDF 校正完成！") + else: + QMessageBox.warning(self, "警告", f"FlexBRDF 校正失败:\n{message}") + + def build_config(self): + """构建FlexBRDF配置字典""" + config = {} + config["input_files"] = [rfl for rfl, _ in self.file_pairs] + config["file_type"] = "envi" + config["bad_bands"] = [] + config["num_cpus"] = self.num_cpus_spin.value() + + anc_files_dict = {} + for rfl_path, obs_path in self.file_pairs: + mapping = self.band_mappings.get(rfl_path, {}) + anc_entry = {} + for param, spin in self.param_spinboxes.items(): + band_idx = mapping.get(param, spin.value()) + anc_entry[param] = [obs_path, band_idx] + anc_files_dict[rfl_path] = anc_entry + config["anc_files"] = anc_files_dict + + corrections = [] + if self.corr_topo_cb.isChecked(): + corrections.append("topo") + if self.corr_brdf_cb.isChecked(): + corrections.append("brdf") + if self.corr_glint_cb.isChecked(): + corrections.append("glint") + config["corrections"] = corrections + + config["topo"] = { + "type": self.topo_type.currentText(), + "calc_mask": self.topo_calc_mask_list.get_conditions(), + "apply_mask": self.topo_apply_mask_list.get_conditions(), + "c_fit_type": self.c_fit_type.currentText() + } + + config["brdf"] = { + "type": self.brdf_type.currentText(), + "grouped": self.brdf_grouped.isChecked(), + "geometric": self.geometric_kernel.currentText(), + "volume": self.volume_kernel.currentText(), + "b/r": self.br_b_ratio.value(), + "h/b": self.hb_ratio.value(), + "sample_perc": self.sample_perc.value(), + "interp_kind": self.interp_kind.currentText(), + "calc_mask": self.brdf_calc_mask_list.get_conditions(), + "apply_mask": self.brdf_apply_mask_list.get_conditions(), + "bin_type": self.bin_type.currentText(), + "num_bins": self.num_bins.value(), + "ndvi_bin_min": self.ndvi_bin_min.value(), + "ndvi_bin_max": self.ndvi_bin_max.value(), + "ndvi_perc_min": self.ndvi_perc_min.value(), + "ndvi_perc_max": self.ndvi_perc_max.value(), + "solar_zn_type": self.solar_zn_type.currentText() + } + + config["glint"] = { + "type": self.glint_type.currentText(), + "correction_wave": self.correction_wave.value(), + "apply_mask": self.glint_apply_mask_list.get_conditions() + } + + subset_waves = [] + if self.subset_waves_edit.text().strip(): + try: + subset_waves = [int(x.strip()) for x in self.subset_waves_edit.text().split(",") if x.strip()] + except: + pass + config["export"] = { + "coeffs": self.coeffs_cb.isChecked(), + "image": self.image_cb.isChecked(), + "masks": self.masks_cb.isChecked(), + "subset_waves": subset_waves, + "output_dir": self.output_dir_edit.text(), + "suffix": self.suffix_edit.text() + } + config["resample"] = self.resample_cb.isChecked() + + return config # ==================== 文件与波段映射页 ==================== def setup_file_tab(self): - layout = QVBoxLayout(self.file_tab) + # 创建滚动区域作为主布局 + scroll = QScrollArea() + scroll.setWidgetResizable(True) + scroll.setStyleSheet(""" + QScrollArea { + border: none; + background-color: #000000; + } + QScrollBar:vertical { + background-color: #1a1a1a; + width: 12px; + border-radius: 6px; + } + QScrollBar::handle:vertical { + background-color: #5e5e5e; + border-radius: 6px; + min-height: 20px; + } + QScrollBar::handle:vertical:hover { + background-color: #76b900; + } + QScrollBar::add-line:vertical, QScrollBar::sub-line:vertical { + height: 0px; + } + QScrollBar:horizontal { + background-color: #1a1a1a; + height: 12px; + border-radius: 6px; + } + QScrollBar::handle:horizontal { + background-color: #5e5e5e; + border-radius: 6px; + min-width: 20px; + } + QScrollBar::handle:horizontal:hover { + background-color: #76b900; + } + QScrollBar::add-line:horizontal, QScrollBar::sub-line:horizontal { + width: 0px; + } + """) + + # 创建内容容器 + content_widget = QWidget() + layout = QVBoxLayout(content_widget) layout.setSpacing(15) layout.setContentsMargins(15, 15, 15, 15) + + # ===== 模式选择 ===== + mode_group = QGroupBox("文件选择模式") + mode_group.setStyleSheet(""" + QGroupBox { + background-color: #1a1a1a; + border: 2px solid #76b900; + border-radius: 2px; + color: #ffffff; + font-weight: 700; + font-size: 14px; + padding: 10px; + margin-top: 10px; + } + QGroupBox::title { + subcontrol-origin: margin; + left: 10px; + padding: 5px 10px; + color: #000000; + background-color: #76b900; + border-radius: 2px; + font-weight: 700; + } + QRadioButton { + color: #ffffff; + font-size: 14px; + spacing: 8px; + } + QRadioButton::indicator { + width: 18px; + height: 18px; + border: 2px solid #76b900; + border-radius: 9px; + background-color: transparent; + } + QRadioButton::indicator:checked { + background-color: #76b900; + } + """) + mode_layout = QHBoxLayout() + self.mode_auto_radio = QRadioButton("自动匹配文件夹") + self.mode_auto_radio.setChecked(True) + self.mode_auto_radio.toggled.connect(self.on_mode_changed) + self.mode_manual_radio = QRadioButton("手动指定文件对") + self.mode_manual_radio.setChecked(False) + self.mode_manual_radio.toggled.connect(self.on_mode_changed) + mode_layout.addWidget(self.mode_auto_radio) + mode_layout.addWidget(self.mode_manual_radio) + mode_layout.addStretch() + mode_group.setLayout(mode_layout) + layout.addWidget(mode_group) + + # 提示标签 + hint_label = QLabel("提示：自动匹配会从文件夹批量导入文件对，手动模式可逐个指定。文件对列表中会标记 [自动] 或 [手动]。") + hint_label.setStyleSheet("color: #aaaaaa; font-size: 11px; margin-top: 5px;") + hint_label.setWordWrap(True) + layout.addWidget(hint_label) + + # ===== 自动匹配区域 ===== + self.auto_match_widget = QWidget() + auto_layout = QVBoxLayout(self.auto_match_widget) + auto_layout.setContentsMargins(0, 0, 0, 0) + auto_layout.setSpacing(10) # 文件夹设置组 - 使用 NVIDIA 绿色主题 folder_group = QGroupBox("文件夹设置 (Folder Settings)") @@ -470,6 +1164,9 @@ class MainWindow(QMainWindow): folder_layout = QFormLayout() folder_layout.setVerticalSpacing(10) folder_layout.setHorizontalSpacing(10) + folder_layout.setFieldGrowthPolicy(QFormLayout.ExpandingFieldsGrow) + folder_layout.setLabelAlignment(Qt.AlignLeft) + folder_layout.setFormAlignment(Qt.AlignLeft) self.input_folder_edit = QLineEdit() self.input_folder_edit.setPlaceholderText("选择反射率数据文件夹...") self.anc_folder_edit = QLineEdit() @@ -481,7 +1178,7 @@ class MainWindow(QMainWindow): folder_layout.addRow("输入文件夹 (反射率):", self._horizontal_widget(self.input_folder_edit, btn_input)) folder_layout.addRow("辅助文件夹 (观测几何):", self._horizontal_widget(self.anc_folder_edit, btn_anc)) folder_group.setLayout(folder_layout) - layout.addWidget(folder_group) + auto_layout.addWidget(folder_group) # 匹配按钮 match_btn = QPushButton("自动匹配文件对 (Auto Match File Pairs)") @@ -501,16 +1198,145 @@ class MainWindow(QMainWindow): } """) match_btn.clicked.connect(self.match_files) - layout.addWidget(match_btn) + auto_layout.addWidget(match_btn) + + layout.addWidget(self.auto_match_widget) - # 文件列表 - list_label = QLabel("匹配的文件对 (点击编辑波段映射):") - list_label.setStyleSheet("color: #a7a7a7; font-size: 13px; margin-top: 10px;") - layout.addWidget(list_label) + # ===== 手动指定区域 ===== + self.manual_widget = QWidget() + manual_layout = QVBoxLayout(self.manual_widget) + manual_layout.setContentsMargins(0, 0, 0, 0) + manual_layout.setSpacing(10) + manual_group = QGroupBox("手动添加文件对") + manual_group.setStyleSheet(""" + QGroupBox { + background-color: #1a1a1a; + border: 2px solid #0080ff; + border-radius: 2px; + color: #ffffff; + font-weight: 700; + font-size: 14px; + padding: 10px; + margin-top: 10px; + } + QGroupBox::title { + subcontrol-origin: margin; + left: 10px; + padding: 5px 10px; + color: #ffffff; + background-color: #0080ff; + border-radius: 2px; + font-weight: 700; + } + QLabel { + color: #ffffff; + font-size: 14px; + font-weight: 700; + } + """) + manual_form = QFormLayout() + manual_form.setVerticalSpacing(10) + manual_form.setFieldGrowthPolicy(QFormLayout.ExpandingFieldsGrow) + manual_form.setLabelAlignment(Qt.AlignLeft) + manual_form.setFormAlignment(Qt.AlignLeft) + + self.manual_rfl_edit = QLineEdit() + self.manual_rfl_edit.setPlaceholderText("选择反射率文件...") + btn_rfl = QPushButton("浏览...") + btn_rfl.clicked.connect(self.select_rfl_file) + manual_form.addRow("反射率文件:", self._horizontal_widget(self.manual_rfl_edit, btn_rfl)) + + self.manual_obs_edit = QLineEdit() + self.manual_obs_edit.setPlaceholderText("选择观测几何文件...") + btn_obs = QPushButton("浏览...") + btn_obs.clicked.connect(self.select_obs_file) + manual_form.addRow("观测几何文件:", self._horizontal_widget(self.manual_obs_edit, btn_obs)) + + manual_group.setLayout(manual_form) + manual_layout.addWidget(manual_group) + + # 添加文件对按钮 + add_pair_btn = QPushButton("添加文件对 (Add File Pair)") + add_pair_btn.setStyleSheet(""" + QPushButton { + background-color: transparent; + color: #ffffff; + border: 2px solid #0080ff; + border-radius: 2px; + padding: 10px 20px; + font-weight: 700; + font-size: 14px; + margin-top: 5px; + } + QPushButton:hover { + background-color: #1eaedb; + } + """) + add_pair_btn.clicked.connect(self.add_manual_file_pair) + manual_layout.addWidget(add_pair_btn) + + self.manual_widget.setVisible(False) # 默认隐藏 + layout.addWidget(self.manual_widget) + + # ===== 文件对列表 ===== + list_group = QGroupBox("文件对列表 (点击编辑波段映射)") + list_group.setStyleSheet(""" + QGroupBox { + background-color: #1a1a1a; + border: 1px solid #5e5e5e; + border-radius: 2px; + color: #ffffff; + font-weight: 700; + font-size: 14px; + padding: 10px; + margin-top: 10px; + } + QGroupBox::title { + subcontrol-origin: margin; + left: 10px; + padding: 0 5px; + color: #ffffff; + } + """) + list_layout = QVBoxLayout(list_group) + + # 文件列表（带滚动区域） + list_scroll = QScrollArea() + list_scroll.setWidgetResizable(True) + list_scroll.setMaximumHeight(200) # 限制最大高度，避免挤压 + list_scroll.setStyleSheet(""" + QScrollArea { + border: none; + background-color: #1a1a1a; + } + QScrollBar:vertical { + background-color: #1a1a1a; + width: 12px; + border-radius: 6px; + } + QScrollBar::handle:vertical { + background-color: #5e5e5e; + border-radius: 6px; + min-height: 20px; + } + QScrollBar::handle:vertical:hover { + background-color: #76b900; + } + QScrollBar::add-line:vertical, QScrollBar::sub-line:vertical { + height: 0px; + } + """) + + list_container = QWidget() + list_container_layout = QVBoxLayout(list_container) + list_container_layout.setContentsMargins(0, 0, 0, 0) + list_container_layout.setSpacing(0) + self.file_list_widget = QListWidget() self.file_list_widget.setSelectionMode(QListWidget.SingleSelection) self.file_list_widget.itemClicked.connect(self.on_file_selected) + self.file_list_widget.setMaximumHeight(400) # 内部最大高度 self.file_list_widget.setStyleSheet(""" QListWidget { background-color: #1a1a1a; @@ -531,7 +1357,31 @@ class MainWindow(QMainWindow): background-color: #3860be; } """) - layout.addWidget(self.file_list_widget) + list_container_layout.addWidget(self.file_list_widget) + list_scroll.setWidget(list_container) + list_layout.addWidget(list_scroll) + + # 删除按钮 + delete_btn = QPushButton("删除选中文件对 (Delete Selected)") + delete_btn.setStyleSheet(""" + QPushButton { + background-color: transparent; + color: #ff4444; + border: 2px solid #ff4444; + border-radius: 2px; + padding: 8px 16px; + font-weight: 700; + font-size: 13px; + } + QPushButton:hover { + background-color: #ff4444; + color: #000000; + } + """) + delete_btn.clicked.connect(self.delete_selected_file_pair) + list_layout.addWidget(delete_btn) + + layout.addWidget(list_group) # 波段映射组 - 使用蓝色主题 mapping_group = QGroupBox("当前选中文件的辅助文件波段索引 (0-based)") @@ -569,29 +1419,69 @@ class MainWindow(QMainWindow): min-width: 80px; } """) - mapping_layout = QFormLayout() - mapping_layout.setVerticalSpacing(8) + mapping_layout = QGridLayout() + mapping_layout.setVerticalSpacing(10) + mapping_layout.setHorizontalSpacing(15) self.param_spinboxes = {} - param_names = [ - "path_length", "sensor_az", "sensor_zn", "solar_az", "solar_zn", - "phase", "slope", "aspect", "cosine_i", "utc_time" + + # 参数名称映射：内部名称 -> 显示名称 + param_names_display = [ + ("path_length", "路径长度 (path_length)"), + ("sensor_az", "传感器方位角 (sensor_az)"), + ("sensor_zn", "传感器天顶角 (sensor_zn)"), + ("solar_az", "太阳方位角 (solar_az)"), + ("solar_zn", "太阳天顶角 (solar_zn)"), + ("phase", "相位角 (phase)"), + ("slope", "坡度 (slope)"), + ("aspect", "坡向 (aspect)"), + ("cosine_i", "入射角余弦 (cosine_i)"), + ("utc_time", "UTC时间 (utc_time)") ] - for name in param_names: - spin = QSpinBox() + + # 使用两列网格布局 + for i, (name, display_name) in enumerate(param_names_display): + row = i // 2 + col = (i % 2) * 2 # 0, 2, 4... for labels + + label = QLabel(f"{display_name}:") + label.setStyleSheet("color: #ffffff; font-size: 12px; min-width: 180px;") + mapping_layout.addWidget(label, row, col) + + spin = NoWheelSpinBox() spin.setRange(0, 1000) spin.setValue(0) + spin.setFixedWidth(80) + spin.setStyleSheet(""" + QSpinBox { + background-color: #000000; + color: #ffffff; + border: 1px solid #0046a4; + border-radius: 2px; + padding: 3px; + } + """) self.param_spinboxes[name] = spin - mapping_layout.addRow(f"{name}:", spin) + mapping_layout.addWidget(spin, row, col + 1) + mapping_group.setLayout(mapping_layout) layout.addWidget(mapping_group) self.current_rfl_path = None + + # 设置滚动区域的内容 + scroll.setWidget(content_widget) + + # 将滚动区域添加到文件标签页 + tab_layout = QVBoxLayout(self.file_tab) + tab_layout.setContentsMargins(0, 0, 0, 0) + tab_layout.addWidget(scroll) def _horizontal_widget(self, line_edit, button): w = QWidget() hbox = QHBoxLayout(w) hbox.setContentsMargins(0, 0, 0, 0) - hbox.addWidget(line_edit) + hbox.setSpacing(10) + hbox.addWidget(line_edit, 1) hbox.addWidget(button) return w @@ -624,28 +1514,156 @@ class MainWindow(QMainWindow): key = basename[:-4] if basename.endswith("_rfl") else basename if key in anc_files: self.file_pairs.append((rfl_path, anc_files[key])) - self.file_list_widget.addItem(f"{basename} -> {os.path.basename(anc_files[key])}") + self.file_list_widget.addItem(f"[自动] {basename} -> {os.path.basename(anc_files[key])}") else: - self.file_list_widget.addItem(f"{basename} -> (未找到匹配的 obs 文件)") + self.file_list_widget.addItem(f"[自动] {basename} -> (未找到匹配的 obs 文件)") if self.file_pairs: self.file_list_widget.setCurrentRow(0) self.on_file_selected(self.file_list_widget.item(0)) + def on_mode_changed(self): + """切换文件选择模式 - 互斥单选""" + auto_mode = self.mode_auto_radio.isChecked() + + # 显示/隐藏对应区域 + self.auto_match_widget.setVisible(auto_mode) + self.manual_widget.setVisible(not auto_mode) + + def select_rfl_file(self): + """手动选择反射率文件""" + file_path, _ = QFileDialog.getOpenFileName( + self, + "选择反射率文件", + "", + "ENVI Files (*.dat *.img *.hdr);;All Files (*)" + ) + if file_path: + self.manual_rfl_edit.setText(file_path) + + def select_obs_file(self): + """手动选择观测几何文件""" + file_path, _ = QFileDialog.getOpenFileName( + self, + "选择观测几何文件", + "", + "ENVI Files (*.dat *.img *.hdr *.bip *.bsq *.bil);;All Files (*)" + ) + if file_path: + self.manual_obs_edit.setText(file_path) + + def add_manual_file_pair(self): + """添加手动指定的文件对""" + rfl_path = self.manual_rfl_edit.text().strip() + obs_path = self.manual_obs_edit.text().strip() + + # 验证输入 + if not rfl_path: + QMessageBox.warning(self, "警告", "请选择反射率文件") + return + if not obs_path: + QMessageBox.warning(self, "警告", "请选择观测几何文件") + return + + # 检查文件是否存在 + if not os.path.exists(rfl_path): + QMessageBox.warning(self, "警告", f"反射率文件不存在：\n{rfl_path}") + return + if not os.path.exists(obs_path): + QMessageBox.warning(self, "警告", f"观测几何文件不存在：\n{obs_path}") + return + + # 检查是否已存在相同的文件对 + for existing_rfl, existing_obs in self.file_pairs: + if existing_rfl == rfl_path and existing_obs == obs_path: + QMessageBox.information(self, "提示", "该文件对已存在") + return + + # 添加文件对 + self.file_pairs.append((rfl_path, obs_path)) + rfl_basename = os.path.basename(rfl_path) + obs_basename = os.path.basename(obs_path) + self.file_list_widget.addItem(f"[手动] {rfl_basename} -> {obs_basename}") + + # 为新文件初始化默认波段映射 + if rfl_path not in self.band_mappings: + self.band_mappings[rfl_path] = { + "path_length": 0, "sensor_az": 9, "sensor_zn": 8, + "solar_az": 7, "solar_zn": 6, "phase": 0, + "slope": 0, "aspect": 0, "cosine_i": 0, "utc_time": 0 + } + + # 选中新添加的文件对 + last_idx = len(self.file_pairs) - 1 + self.file_list_widget.setCurrentRow(last_idx) + self.on_file_selected(self.file_list_widget.item(last_idx)) + + # 清空输入框 + self.manual_rfl_edit.clear() + self.manual_obs_edit.clear() + + QMessageBox.information(self, "成功", "文件对已添加") + + def delete_selected_file_pair(self): + """删除选中的文件对""" + idx = self.file_list_widget.currentRow() + if idx < 0 or idx >= len(self.file_pairs): + QMessageBox.warning(self, "警告", "请先选择一个文件对") + return + + rfl_path, obs_path = self.file_pairs[idx] + rfl_name = os.path.basename(rfl_path) + + # 确认删除 + reply = QMessageBox.question( + self, "确认删除", + f"确定要删除文件对吗？\n\n{rfl_name}", + QMessageBox.Yes | QMessageBox.No, + QMessageBox.No + ) + + if reply == QMessageBox.Yes: + # 从列表和数据结构中删除 + self.file_pairs.pop(idx) + self.file_list_widget.takeItem(idx) + + # 清除对应的波段映射 + if rfl_path in self.band_mappings: + del self.band_mappings[rfl_path] + + # 清除当前选中 + self.current_rfl_path = None + for spin in self.param_spinboxes.values(): + spin.setValue(0) + + # 如果有剩余文件对，选中第一个 + if self.file_pairs: + self.file_list_widget.setCurrentRow(0) + self.on_file_selected(self.file_list_widget.item(0)) + + QMessageBox.information(self, "成功", "文件对已删除") + def on_file_selected(self, item): idx = self.file_list_widget.row(item) if idx >= len(self.file_pairs): return rfl_path, obs_path = self.file_pairs[idx] self.current_rfl_path = rfl_path + + # 更新波段映射spinbox值 if rfl_path in self.band_mappings: mapping = self.band_mappings[rfl_path] for param, spin in self.param_spinboxes.items(): - spin.setValue(mapping.get(param, 0)) + value = mapping.get(param, 0) + spin.setValue(value) else: for spin in self.param_spinboxes.values(): spin.setValue(0) + + # 强制刷新UI + QApplication.processEvents() + # 重新连接信号 for param, spin in self.param_spinboxes.items(): try: spin.valueChanged.disconnect() @@ -712,8 +1730,11 @@ class MainWindow(QMainWindow): """) corr_layout = QHBoxLayout() self.corr_topo_cb = QCheckBox("Topo (地形校正)") + self.corr_topo_cb.stateChanged.connect(self.on_correction_changed) self.corr_brdf_cb = QCheckBox("BRDF (双向反射分布函数)") + self.corr_brdf_cb.stateChanged.connect(self.on_correction_changed) self.corr_glint_cb = QCheckBox("Glint (耀斑校正)") + self.corr_glint_cb.stateChanged.connect(self.on_correction_changed) corr_layout.addWidget(self.corr_topo_cb) corr_layout.addWidget(self.corr_brdf_cb) corr_layout.addWidget(self.corr_glint_cb) @@ -721,8 +1742,9 @@ class MainWindow(QMainWindow): layout.addWidget(corr_group) # 2. Topo 参数 - 蓝色主题 (#0046a4) - topo_group = QGroupBox("地形校正 (Topo Correction)") - topo_group.setStyleSheet(""" + self.topo_group = QGroupBox("地形校正 (Topo Correction)") + self.topo_group.setVisible(False) # 默认隐藏 + self.topo_group.setStyleSheet(""" QGroupBox { background-color: #1a1a1a; border: 2px solid #0046a4; @@ -778,19 +1800,23 @@ class MainWindow(QMainWindow): topo_layout.setSpacing(10) # type type_layout = QHBoxLayout() + type_layout.setSpacing(10) type_label = QLabel("类型 (type):") type_layout.addWidget(type_label) - self.topo_type = QComboBox() + self.topo_type = NoWheelComboBox() self.topo_type.addItems(["mod_minneart", "scs+c", "cosine", "c", "scs"]) + self.topo_type.setMinimumWidth(250) type_layout.addWidget(self.topo_type) type_layout.addStretch() topo_layout.addLayout(type_layout) # c_fit_type fit_layout = QHBoxLayout() + fit_layout.setSpacing(10) fit_label = QLabel("拟合类型 (c_fit_type):") fit_layout.addWidget(fit_label) - self.c_fit_type = QComboBox() + self.c_fit_type = NoWheelComboBox() self.c_fit_type.addItems(["nnls", "ols", "wls"]) + self.c_fit_type.setMinimumWidth(250) fit_layout.addWidget(self.c_fit_type) fit_layout.addStretch() topo_layout.addLayout(fit_layout) @@ -805,12 +1831,13 @@ class MainWindow(QMainWindow): topo_layout.addWidget(apply_mask_label) self.topo_apply_mask_list = MaskListWidget() topo_layout.addWidget(self.topo_apply_mask_list) - topo_group.setLayout(topo_layout) - layout.addWidget(topo_group) + self.topo_group.setLayout(topo_layout) + layout.addWidget(self.topo_group) # 3. BRDF 参数 - NVIDIA 绿色主题 (#76b900) - brdf_group = QGroupBox("BRDF 校正 (BRDF Correction)") - brdf_group.setStyleSheet(""" + self.brdf_group = QGroupBox("BRDF 校正 (BRDF Correction)") + self.brdf_group.setVisible(False) # 默认隐藏 + self.brdf_group.setStyleSheet(""" QGroupBox { background-color: #1a1a1a; border: 2px solid #76b900; @@ -888,35 +1915,46 @@ class MainWindow(QMainWindow): """) brdf_layout = QFormLayout() brdf_layout.setVerticalSpacing(10) - self.brdf_type = QComboBox() + brdf_layout.setFieldGrowthPolicy(QFormLayout.ExpandingFieldsGrow) + brdf_layout.setLabelAlignment(Qt.AlignLeft) + brdf_layout.setFormAlignment(Qt.AlignLeft) + + self.brdf_type = NoWheelComboBox() self.brdf_type.addItems(["flex", "universal"]) + self.brdf_type.setMinimumWidth(300) brdf_layout.addRow("算法类型 (type):", self.brdf_type) self.brdf_grouped = QCheckBox("启用分组 (grouped)") self.brdf_grouped.setChecked(True) brdf_layout.addRow("", self.brdf_grouped) - self.geometric_kernel = QComboBox() + self.geometric_kernel = NoWheelComboBox() self.geometric_kernel.addItems(["li_sparse", "li_dense", "li_sparse_r", "li_dense_r", "roujean"]) + self.geometric_kernel.setMinimumWidth(300) brdf_layout.addRow("几何核 (geometric):", self.geometric_kernel) - self.volume_kernel = QComboBox() + self.volume_kernel = NoWheelComboBox() self.volume_kernel.addItems(["ross_thin", "ross_thick", "hotspot", "roujean"]) + self.volume_kernel.setMinimumWidth(300) brdf_layout.addRow("体积核 (volume):", self.volume_kernel) - self.br_b_ratio = QDoubleSpinBox() + self.br_b_ratio = NoWheelDoubleSpinBox() self.br_b_ratio.setRange(0.1, 10.0) self.br_b_ratio.setValue(2.5) self.br_b_ratio.setDecimals(2) + self.br_b_ratio.setMinimumWidth(300) brdf_layout.addRow("b/r 比率:", self.br_b_ratio) - self.hb_ratio = QDoubleSpinBox() + self.hb_ratio = NoWheelDoubleSpinBox() self.hb_ratio.setRange(0.1, 5.0) self.hb_ratio.setValue(2.0) self.hb_ratio.setDecimals(2) + self.hb_ratio.setMinimumWidth(300) brdf_layout.addRow("h/b 比率:", self.hb_ratio) - self.sample_perc = QDoubleSpinBox() + self.sample_perc = NoWheelDoubleSpinBox() self.sample_perc.setRange(0.01, 1.0) self.sample_perc.setValue(0.1) self.sample_perc.setDecimals(2) + self.sample_perc.setMinimumWidth(300) brdf_layout.addRow("采样比例 (sample_perc):", self.sample_perc) - self.interp_kind = QComboBox() + self.interp_kind = NoWheelComboBox() self.interp_kind.addItems(["linear", "nearest", "cubic"]) + self.interp_kind.setMinimumWidth(300) brdf_layout.addRow("插值方法 (interp_kind):", self.interp_kind) # BRDF 掩膜 calc_mask_label_brdf = QLabel("计算掩膜条件 (calc_mask):") @@ -930,40 +1968,48 @@ class MainWindow(QMainWindow): self.brdf_apply_mask_list = MaskListWidget() brdf_layout.addRow(self.brdf_apply_mask_list) # bin_type 等 - self.bin_type = QComboBox() + self.bin_type = NoWheelComboBox() self.bin_type.addItems(["dynamic", "fixed"]) + self.bin_type.setMinimumWidth(300) brdf_layout.addRow("分箱类型 (bin_type):", self.bin_type) - self.num_bins = QSpinBox() + self.num_bins = NoWheelSpinBox() self.num_bins.setRange(1, 100) self.num_bins.setValue(18) + self.num_bins.setMinimumWidth(300) brdf_layout.addRow("分箱数量 (num_bins):", self.num_bins) - self.ndvi_bin_min = QDoubleSpinBox() + self.ndvi_bin_min = NoWheelDoubleSpinBox() self.ndvi_bin_min.setRange(0.0, 1.0) self.ndvi_bin_min.setValue(0.05) self.ndvi_bin_min.setDecimals(2) + self.ndvi_bin_min.setMinimumWidth(300) brdf_layout.addRow("NDVI 最小值 (ndvi_bin_min):", self.ndvi_bin_min) - self.ndvi_bin_max = QDoubleSpinBox() + self.ndvi_bin_max = NoWheelDoubleSpinBox() self.ndvi_bin_max.setRange(0.0, 1.0) self.ndvi_bin_max.setValue(1.0) self.ndvi_bin_max.setDecimals(2) + self.ndvi_bin_max.setMinimumWidth(300) brdf_layout.addRow("NDVI 最大值 (ndvi_bin_max):", self.ndvi_bin_max) - self.ndvi_perc_min = QSpinBox() + self.ndvi_perc_min = NoWheelSpinBox() self.ndvi_perc_min.setRange(0, 100) self.ndvi_perc_min.setValue(10) + self.ndvi_perc_min.setMinimumWidth(300) brdf_layout.addRow("NDVI 百分比最小值 (ndvi_perc_min):", self.ndvi_perc_min) - self.ndvi_perc_max = QSpinBox() + self.ndvi_perc_max = NoWheelSpinBox() self.ndvi_perc_max.setRange(0, 100) self.ndvi_perc_max.setValue(95) + self.ndvi_perc_max.setMinimumWidth(300) brdf_layout.addRow("NDVI 百分比最大值 (ndvi_perc_max):", self.ndvi_perc_max) - self.solar_zn_type = QComboBox() + self.solar_zn_type = NoWheelComboBox() self.solar_zn_type.addItems(["scene", "pixel"]) + self.solar_zn_type.setMinimumWidth(300) brdf_layout.addRow("太阳天顶角类型 (solar_zn_type):", self.solar_zn_type) - brdf_group.setLayout(brdf_layout) - layout.addWidget(brdf_group) + self.brdf_group.setLayout(brdf_layout) + layout.addWidget(self.brdf_group) # 4. Glint 参数 - 橙色主题 (#df6500) - glint_group = QGroupBox("耀斑校正 (Glint Correction)") - glint_group.setStyleSheet(""" + self.glint_group = QGroupBox("耀斑校正 (Glint Correction)") + self.glint_group.setVisible(False) # 默认隐藏 + self.glint_group.setStyleSheet(""" QGroupBox { background-color: #1a1a1a; border: 2px solid #df6500; @@ -1024,20 +2070,26 @@ class MainWindow(QMainWindow): """) glint_layout = QFormLayout() glint_layout.setVerticalSpacing(10) - self.glint_type = QComboBox() + glint_layout.setFieldGrowthPolicy(QFormLayout.ExpandingFieldsGrow) + glint_layout.setLabelAlignment(Qt.AlignLeft) + glint_layout.setFormAlignment(Qt.AlignLeft) + + self.glint_type = NoWheelComboBox() self.glint_type.addItems(["hochberg", "gao", "hedley"]) + self.glint_type.setMinimumWidth(300) glint_layout.addRow("校正算法 (type):", self.glint_type) - self.correction_wave = QSpinBox() + self.correction_wave = NoWheelSpinBox() self.correction_wave.setRange(1000, 2500) self.correction_wave.setValue(2150) + self.correction_wave.setMinimumWidth(300) glint_layout.addRow("校正波长 (nm):", self.correction_wave) apply_mask_label_glint = QLabel("应用掩膜条件 (apply_mask):") apply_mask_label_glint.setStyleSheet("color: #a7a7a7; font-size: 13px;") glint_layout.addRow(apply_mask_label_glint) self.glint_apply_mask_list = MaskListWidget() glint_layout.addRow(self.glint_apply_mask_list) - glint_group.setLayout(glint_layout) - layout.addWidget(glint_group) + self.glint_group.setLayout(glint_layout) + layout.addWidget(self.glint_group) layout.addStretch() self.corr_tab.setLayout(QVBoxLayout()) @@ -1045,6 +2097,19 @@ class MainWindow(QMainWindow): # 设置容器背景色 container.setStyleSheet("background-color: #000000;") + + # 初始化校正参数组显示状态 + self.on_correction_changed() + + def on_correction_changed(self): + """根据选择的校正类型显示/隐藏对应的参数组""" + # 显示/隐藏对应的参数组 + if hasattr(self, 'topo_group'): + self.topo_group.setVisible(self.corr_topo_cb.isChecked()) + if hasattr(self, 'brdf_group'): + self.brdf_group.setVisible(self.corr_brdf_cb.isChecked()) + if hasattr(self, 'glint_group'): + self.glint_group.setVisible(self.corr_glint_cb.isChecked()) # ==================== 输出设置页 ==================== def setup_output_tab(self): @@ -1052,6 +2117,9 @@ class MainWindow(QMainWindow): layout.setVerticalSpacing(12) layout.setHorizontalSpacing(10) layout.setContentsMargins(15, 15, 15, 15) + layout.setFieldGrowthPolicy(QFormLayout.ExpandingFieldsGrow) + layout.setLabelAlignment(Qt.AlignLeft) + layout.setFormAlignment(Qt.AlignLeft) # 设置输出设置页的样式 self.output_tab.setStyleSheet(""" @@ -1071,6 +2139,17 @@ class MainWindow(QMainWindow): QLineEdit:focus { border: 2px solid #76b900; } + QSpinBox { + background-color: #1a1a1a; + color: #ffffff; + border: 1px solid #76b900; + border-radius: 2px; + padding: 5px; + min-width: 300px; + } + QSpinBox:focus { + border: 2px solid #76b900; + } QCheckBox { color: #ffffff; font-size: 14px; @@ -1121,6 +2200,13 @@ class MainWindow(QMainWindow): self.subset_waves_edit.setPlaceholderText("例如: 450, 550, 650 (逗号分隔，留空导出全部)") layout.addRow("子集波段 (Subset Waves):", self.subset_waves_edit) + # CPU核心数 + self.num_cpus_spin = NoWheelSpinBox() + self.num_cpus_spin.setRange(1, 64) + self.num_cpus_spin.setValue(2) + self.num_cpus_spin.setMinimumWidth(300) + layout.addRow("CPU核心数 (num_cpus):", self.num_cpus_spin) + # 重采样 self.resample_cb = QCheckBox("启用重采样 (Resample)") layout.addRow("", self.resample_cb) @@ -1137,24 +2223,56 @@ class MainWindow(QMainWindow): QMessageBox.critical(self, "错误", f"无法加载 JSON: {e}") return - # 文件与波段映射（只恢复映射关系，不自动恢复文件对） + # 文件与波段映射 - 从 anc_files 恢复文件对和波段映射 if "anc_files" in config: self.band_mappings = {} + self.file_pairs = [] + self.file_list_widget.clear() + for rfl_path, anc_info in config["anc_files"].items(): + # 恢复波段映射 mapping = {} + obs_path = None for param, val in anc_info.items(): if isinstance(val, list) and len(val) >= 2: mapping[param] = val[1] + # 提取第一个参数的obs路径作为该文件对的obs路径 + if obs_path is None: + obs_path = val[0] self.band_mappings[rfl_path] = mapping - # 刷新当前选中的文件对映射显示 - if self.file_pairs and self.file_list_widget.currentItem(): - self.on_file_selected(self.file_list_widget.currentItem()) + + # 恢复文件对（不检查文件是否存在，支持加载其他机器创建的配置） + if obs_path: + self.file_pairs.append((rfl_path, obs_path)) + rfl_basename = os.path.basename(rfl_path) + obs_basename = os.path.basename(obs_path) + # 标记文件是否存在 + exists_marker = "" if os.path.exists(rfl_path) else " [文件不存在]" + self.file_list_widget.addItem(f"{rfl_basename} -> {obs_basename}{exists_marker}") + + # 自动选择第一个文件并刷新显示 + if self.file_pairs: + self.file_list_widget.setCurrentRow(0) + # 先设置文件夹路径（在on_file_selected之前） + first_rfl_path = self.file_pairs[0][0] + first_obs_path = self.file_pairs[0][1] + input_dir = os.path.dirname(first_rfl_path) + anc_dir = os.path.dirname(first_obs_path) + self.input_folder_edit.setText(input_dir) + self.anc_folder_edit.setText(anc_dir) + # 触发UI更新 + QApplication.processEvents() + # 然后选择文件并刷新波段映射 + self.on_file_selected(self.file_list_widget.item(0)) # corrections corrections = config.get("corrections", []) self.corr_topo_cb.setChecked("topo" in corrections) self.corr_brdf_cb.setChecked("brdf" in corrections) self.corr_glint_cb.setChecked("glint" in corrections) + + # 更新参数组显示状态 + self.on_correction_changed() # topo topo = config.get("topo", {}) @@ -1198,6 +2316,7 @@ class MainWindow(QMainWindow): self.masks_cb.setChecked(export.get("masks", True)) subset = export.get("subset_waves", []) self.subset_waves_edit.setText(",".join(str(w) for w in subset)) + self.num_cpus_spin.setValue(config.get("num_cpus", 2)) self.resample_cb.setChecked(config.get("resample", False)) QMessageBox.information(self, "成功", f"已加载配置: {file_path}") @@ -1211,7 +2330,7 @@ class MainWindow(QMainWindow): config["input_files"] = [rfl for rfl, _ in self.file_pairs] config["file_type"] = "envi" config["bad_bands"] = [] - config["num_cpus"] = 2 + config["num_cpus"] = self.num_cpus_spin.value() anc_files_dict = {} for rfl_path, obs_path in self.file_pairs: @@ -1292,6 +2411,467 @@ class MainWindow(QMainWindow): QMessageBox.critical(self, "错误", f"保存失败: {e}") +# ==================== 水体BRDF模块窗口 ==================== +class WaterBRDFWindow(QMainWindow): + def __init__(self): + super().__init__() + self.setWindowTitle("水体 BRDF 校正 - Ocean BRDF") + self.setMinimumSize(1000, 700) + + # 命令执行相关 + self.command_thread = None + + # 设置样式和界面 + self.setup_water_styles() + self.setup_water_interface() + + def setup_water_styles(self): + """设置水体BRDF模块的样式""" + self.setStyleSheet(""" + QMainWindow { + background-color: #000000; + } + QWidget { + background-color: #000000; + color: #ffffff; + font-family: Arial, Helvetica, sans-serif; + font-size: 14px; + } + QGroupBox { + background-color: #1a1a1a; + border: 2px solid #0080ff; + border-radius: 8px; + margin-top: 15px; + padding-top: 15px; + font-weight: bold; + font-size: 16px; + } + QGroupBox::title { + subcontrol-origin: margin; + left: 15px; + padding: 5px 10px; + color: #ffffff; + background-color: #0080ff; + border-radius: 4px; + font-weight: bold; + } + QPushButton { + background-color: transparent; + color: #ffffff; + border: 2px solid #0080ff; + border-radius: 6px; + padding: 10px 20px; + font-weight: bold; + font-size: 14px; + } + QPushButton:hover { + background-color: #0080ff; + color: #000000; + } + QPushButton:pressed { + background-color: #0066cc; + border-color: #0066cc; + } + QLineEdit { + background-color: #1a1a1a; + color: #ffffff; + border: 1px solid #0080ff; + border-radius: 4px; + padding: 8px; + font-size: 14px; + } + QLineEdit:focus { + border: 2px solid #0080ff; + } + QComboBox { + background-color: #1a1a1a; + color: #ffffff; + border: 1px solid #0080ff; + border-radius: 4px; + padding: 8px; + font-size: 14px; + } + QComboBox::drop-down { + border: none; + width: 25px; + } + QComboBox::down-arrow { + image: none; + border-left: 5px solid transparent; + border-right: 5px solid transparent; + border-top: 6px solid #0080ff; + margin-right: 8px; + } + QComboBox QAbstractItemView { + background-color: #1a1a1a; + color: #ffffff; + border: 1px solid #0080ff; + selection-background-color: #0080ff; + } + QCheckBox { + color: #ffffff; + font-size: 14px; + spacing: 8px; + } + QCheckBox::indicator { + width: 18px; + height: 18px; + border: 2px solid #0080ff; + border-radius: 2px; + background-color: transparent; + } + QCheckBox::indicator:checked { + background-color: #0080ff; + } + QSpinBox, QDoubleSpinBox { + background-color: #1a1a1a; + color: #ffffff; + border: 1px solid #0080ff; + border-radius: 4px; + padding: 5px; + } + QTextEdit { + background-color: #1a1a1a; + color: #ffffff; + border: 1px solid #0080ff; + border-radius: 4px; + font-family: 'Consolas', 'Monaco', monospace; + font-size: 12px; + } + QProgressBar { + background-color: #1a1a1a; + border: 1px solid #0080ff; + border-radius: 4px; + text-align: center; + color: #ffffff; + } + QProgressBar::chunk { + background-color: #0080ff; + border-radius: 4px; + } + """) + + def setup_water_interface(self): + """创建水体BRDF模块界面""" + central_widget = QWidget() + self.setCentralWidget(central_widget) + main_layout = QVBoxLayout(central_widget) + main_layout.setContentsMargins(30, 30, 30, 30) + main_layout.setSpacing(20) + + # 标题 + title_label = QLabel("水体 BRDF 校正") + title_label.setAlignment(Qt.AlignCenter) + title_label.setStyleSheet(""" + font-size: 28px; + font-weight: bold; + color: #0080ff; + margin-bottom: 10px; + """) + main_layout.addWidget(title_label) + + subtitle_label = QLabel("Ocean BRDF Correction - 支持 L11, M02, M02SeaDAS, O25 模型") + subtitle_label.setAlignment(Qt.AlignCenter) + subtitle_label.setStyleSheet(""" + font-size: 14px; + color: #cccccc; + margin-bottom: 20px; + """) + main_layout.addWidget(subtitle_label) + + # 创建分割器：左侧配置，右侧日志 + splitter = QSplitter(Qt.Horizontal) + + # 左侧：参数配置 + left_widget = QWidget() + left_layout = QVBoxLayout(left_widget) + + # 输入文件组 + input_group = QGroupBox("输入文件") + input_layout = QFormLayout(input_group) + input_layout.setFieldGrowthPolicy(QFormLayout.ExpandingFieldsGrow) + input_layout.setLabelAlignment(Qt.AlignLeft) + input_layout.setFormAlignment(Qt.AlignLeft) + + self.hyperspectral_edit = QLineEdit() + self.hyperspectral_edit.setPlaceholderText("选择高光谱数据文件 (.img, .dat, .hdr)") + hyperspectral_btn = QPushButton("浏览...") + hyperspectral_btn.clicked.connect(lambda: self.browse_file(self.hyperspectral_edit, "高光谱数据文件", "ENVI Files (*.img *.dat *.hdr);;All Files (*)")) + input_layout.addRow("高光谱文件:", self.create_file_row(self.hyperspectral_edit, hyperspectral_btn)) + + self.angle_edit = QLineEdit() + self.angle_edit.setPlaceholderText("选择角度文件 (.dat, .bsq, .bip, .bil, .hdr)") + angle_btn = QPushButton("浏览...") + angle_btn.clicked.connect(lambda: self.browse_file(self.angle_edit, "角度文件", "ENVI Files (*.dat *.bsq *.bip *.bil *.hdr);;All Files (*)")) + input_layout.addRow("角度文件:", self.create_file_row(self.angle_edit, angle_btn)) + + self.mask_edit = QLineEdit() + self.mask_edit.setPlaceholderText("选择水体掩膜文件 (.tif, .tiff)") + mask_btn = QPushButton("浏览...") + mask_btn.clicked.connect(lambda: self.browse_file(self.mask_edit, "水体掩膜文件", "GeoTIFF Files (*.tif *.tiff);;All Files (*)")) + input_layout.addRow("水体掩膜:", self.create_file_row(self.mask_edit, mask_btn)) + + left_layout.addWidget(input_group) + + # BRDF模型参数组 + model_group = QGroupBox("BRDF 模型参数") + model_layout = QFormLayout(model_group) + model_layout.setFieldGrowthPolicy(QFormLayout.ExpandingFieldsGrow) + model_layout.setLabelAlignment(Qt.AlignLeft) + model_layout.setFormAlignment(Qt.AlignLeft) + + self.brdf_model_combo = NoWheelComboBox() + self.brdf_model_combo.addItems(["L11", "M02", "O25"]) + self.brdf_model_combo.setCurrentText("L11") + model_layout.addRow("BRDF 模型:", self.brdf_model_combo) + + # 输出变量多选列表 + output_vars_label = QLabel("输出变量:") + output_vars_label.setStyleSheet("color: #ffffff; font-weight: bold;") + model_layout.addRow(output_vars_label) + + self.output_vars_list = QListWidget() + self.output_vars_list.setSelectionMode(QListWidget.MultiSelection) + self.output_vars_list.addItems([ + "Rw_brdf - BRDF校正后的水体反射率", + "rho_ex_w - 水体反射率 (nrrs * π)", + "nrrs - 归一化遥感反射率", + "C_brdf - BRDF校正因子", + "brdf_unc - BRDF不确定度", + "nrrs_unc - nrrs不确定度" + ]) + # 默认选中 Rw_brdf, nrrs, C_brdf + self.output_vars_list.item(0).setSelected(True) + self.output_vars_list.item(2).setSelected(True) + self.output_vars_list.item(3).setSelected(True) + self.output_vars_list.setMaximumHeight(120) + self.output_vars_list.setStyleSheet(""" + QListWidget { + background-color: #1a1a1a; + color: #ffffff; + border: 1px solid #0080ff; + border-radius: 4px; + padding: 5px; + } + QListWidget::item { + padding: 5px; + border-radius: 2px; + } + QListWidget::item:selected { + background-color: #0080ff; + color: #ffffff; + } + QListWidget::item:hover { + background-color: #0046a4; + } + """) + model_layout.addRow(self.output_vars_list) + + # 添加提示标签 + hint_label = QLabel("提示：按住 Ctrl 键可多选；brdf_unc 和 nrrs_unc 需要勾选计算不确定性") + hint_label.setStyleSheet("color: #aaaaaa; font-size: 11px; margin-top: 5px;") + hint_label.setWordWrap(True) + model_layout.addRow(hint_label) + + self.chunk_size_spin = NoWheelSpinBox() + self.chunk_size_spin.setRange(1024, 16384) + self.chunk_size_spin.setValue(4096) + model_layout.addRow("处理块大小:", self.chunk_size_spin) + + self.block_size_spin = NoWheelSpinBox() + self.block_size_spin.setRange(256, 2048) + self.block_size_spin.setValue(512) + model_layout.addRow("空间块大小:", self.block_size_spin) + + self.compute_uncertainty_cb = QCheckBox("计算不确定性") + model_layout.addRow("", self.compute_uncertainty_cb) + + left_layout.addWidget(model_group) + + # 输出设置组 + output_group = QGroupBox("输出设置") + output_layout = QFormLayout(output_group) + output_layout.setFieldGrowthPolicy(QFormLayout.ExpandingFieldsGrow) + output_layout.setLabelAlignment(Qt.AlignLeft) + output_layout.setFormAlignment(Qt.AlignLeft) + + self.output_file_edit = QLineEdit() + self.output_file_edit.setPlaceholderText("输出文件前缀路径") + output_btn = QPushButton("选择...") + output_btn.clicked.connect(self.select_output_file) + output_layout.addRow("输出路径:", self.create_file_row(self.output_file_edit, output_btn)) + + self.output_format_combo = NoWheelComboBox() + self.output_format_combo.addItems(["ENVI"]) + output_layout.addRow("输出格式:", self.output_format_combo) + + left_layout.addWidget(output_group) + + # 运行按钮 + run_btn = QPushButton("运行水体 BRDF 校正") + run_btn.setStyleSheet(""" + QPushButton { + background-color: #0080ff; + color: #ffffff; + border: 2px solid #0080ff; + font-weight: bold; + padding: 12px 24px; + font-size: 16px; + } + QPushButton:hover { + background-color: #0066cc; + border-color: #0066cc; + } + """) + run_btn.clicked.connect(self.run_water_brdf) + left_layout.addWidget(run_btn) + + # 右侧：日志和进度显示 + right_widget = QWidget() + right_layout = QVBoxLayout(right_widget) + + log_label = QLabel("执行日志") + log_label.setStyleSheet("font-weight: bold; color: #0080ff; font-size: 16px;") + right_layout.addWidget(log_label) + + self.log_text = QTextEdit() + self.log_text.setMaximumWidth(400) + self.log_text.setPlainText("等待执行命令...\n") + right_layout.addWidget(self.log_text) + + self.progress_bar = QProgressBar() + self.progress_bar.setVisible(False) + right_layout.addWidget(self.progress_bar) + + # 添加到分割器 + splitter.addWidget(left_widget) + splitter.addWidget(right_widget) + splitter.setStretchFactor(0, 2) + splitter.setStretchFactor(1, 1) + + main_layout.addWidget(splitter) + + def create_file_row(self, line_edit, button): + """创建文件选择行""" + widget = QWidget() + layout = QHBoxLayout(widget) + layout.setContentsMargins(0, 0, 0, 0) + layout.setSpacing(10) + layout.addWidget(line_edit, 1) + layout.addWidget(button) + return widget + + def browse_file(self, line_edit, title, file_filter): + """浏览选择文件""" + file_path, _ = QFileDialog.getOpenFileName(self, f"选择{title}", "", file_filter) + if file_path: + line_edit.setText(file_path) + + def select_output_file(self): + """选择输出文件路径""" + file_path, _ = QFileDialog.getSaveFileName(self, "选择输出路径", "", "All Files (*)") + if file_path: + self.output_file_edit.setText(file_path) + + def run_water_brdf(self): + """运行水体BRDF校正""" + # 检查必需参数 + hyperspectral_file = self.hyperspectral_edit.text().strip() + angle_file = self.angle_edit.text().strip() + mask_file = self.mask_edit.text().strip() + output_file = self.output_file_edit.text().strip() + + if not all([hyperspectral_file, angle_file, mask_file, output_file]): + QMessageBox.warning(self, "警告", "请填写所有必需的文件路径！") + return + + # 检查文件是否存在 + for file_path, name in [(hyperspectral_file, "高光谱文件"), + (angle_file, "角度文件"), + (mask_file, "掩膜文件")]: + if not os.path.exists(file_path): + QMessageBox.warning(self, "警告", f"{name} 不存在：\n{file_path}") + return + + # 构建命令 + ocbrdf_path = os.path.join(os.path.dirname(__file__), "..", "ocbrdf", "ocbrdf_main3.py") + ocbrdf_path = os.path.abspath(ocbrdf_path) + ocbrdf_dir = os.path.dirname(ocbrdf_path) + + if not os.path.exists(ocbrdf_path): + QMessageBox.critical(self, "错误", f"找不到 ocbrdf_main3.py 文件:\n{ocbrdf_path}") + return + + # 使用当前Python解释器（而不是简单的'python'命令） + python_exe = sys.executable + + # 获取选中的输出变量（解析列表项，取空格前的变量名） + selected_items = self.output_vars_list.selectedItems() + output_vars = [] + for item in selected_items: + var_name = item.text().split(' - ')[0].strip() + output_vars.append(var_name) + + if not output_vars: + QMessageBox.warning(self, "警告", "请至少选择一个输出变量！") + return + + # 构建命令参数 - 注意：ocbrdf_main3.py 需要从 ocbrdf 目录运行 + # 以确保能正确导入同级目录的模块 + cmd_parts = [ + f'"{python_exe}"', + f'"{ocbrdf_path}"', + f'"{hyperspectral_file}"', + f'"{angle_file}"', + f'"{mask_file}"', + f'"{output_file}"', + f'--brdf-model {self.brdf_model_combo.currentText()}', + f'--output-var {" ".join(output_vars)}', + f'--output-format {self.output_format_combo.currentText()}', + f'--chunk-size {self.chunk_size_spin.value()}', + f'--block-size {self.block_size_spin.value()}' + ] + + if self.compute_uncertainty_cb.isChecked(): + cmd_parts.append('--compute-uncertainty') + + command = ' '.join(cmd_parts) + + # 显示进度条和日志 + self.progress_bar.setVisible(True) + self.progress_bar.setRange(0, 0) # 不确定进度 + self.log_text.clear() + self.log_text.append("开始执行水体 BRDF 校正...") + self.log_text.append(f"BRDF 模型: {self.brdf_model_combo.currentText()}") + self.log_text.append(f"输出变量: {', '.join(output_vars)}") + self.log_text.append(f"Python: {python_exe}") + self.log_text.append(f"命令: {command}\n") + + # 启动命令执行线程 - 在ocbrdf目录下执行，确保能找到本地模块 + env = {"PYTHONPATH": ocbrdf_dir} # 添加ocbrdf目录到Python路径 + self.command_thread = CommandThread(command, ocbrdf_dir, env=env) + self.command_thread.output_updated.connect(self.append_log) + self.command_thread.finished_signal.connect(self.command_finished) + self.command_thread.start() + + def append_log(self, text): + """添加日志文本""" + self.log_text.append(text) + cursor = self.log_text.textCursor() + cursor.movePosition(cursor.End) + self.log_text.setTextCursor(cursor) + + def command_finished(self, success, message): + """命令执行完成""" + self.progress_bar.setVisible(False) + self.log_text.append(f"\n{message}") + + if success: + QMessageBox.information(self, "成功", "水体 BRDF 校正完成！") + else: + QMessageBox.warning(self, "警告", f"水体 BRDF 校正失败:\n{message}") + + def main(): app = QApplication(sys.argv) window = MainWindow() diff --git a/GUI/temp_flexbrdf_config.json b/GUI/temp_flexbrdf_config.json new file mode 100644 index 0000000..c31a170 --- /dev/null +++ b/GUI/temp_flexbrdf_config.json @@ -0,0 +1,133 @@ +{ + "input_files": [ + "E:/is2/yaopu/output/before.dat" + ], + "file_type": "envi", + "bad_bands": [], + "num_cpus": 10, + "anc_files": { + "E:/is2/yaopu/output/before.dat": { + "path_length": [ + "E:/is2/yaopu/output/angel.dat", + 0 + ], + "sensor_az": [ + "E:/is2/yaopu/output/angel.dat", + 9 + ], + "sensor_zn": [ + "E:/is2/yaopu/output/angel.dat", + 8 + ], + "solar_az": [ + "E:/is2/yaopu/output/angel.dat", + 7 + ], + "solar_zn": [ + "E:/is2/yaopu/output/angel.dat", + 6 + ], + "phase": [ + "E:/is2/yaopu/output/angel.dat", + 0 + ], + "slope": [ + "E:/is2/yaopu/output/angel.dat", + 0 + ], + "aspect": [ + "E:/is2/yaopu/output/angel.dat", + 0 + ], + "cosine_i": [ + "E:/is2/yaopu/output/angel.dat", + 0 + ], + "utc_time": [ + "E:/is2/yaopu/output/angel.dat", + 0 + ] + } + }, + "corrections": [ + "brdf" + ], + "topo": { + "type": "mod_minneart", + "calc_mask": [], + "apply_mask": [], + "c_fit_type": "nnls" + }, + "brdf": { + "type": "flex", + "grouped": true, + "geometric": "li_dense_r", + "volume": "ross_thick", + "b/r": 2.5, + "h/b": 2.0, + "sample_perc": 0.1, + "interp_kind": "linear", + "calc_mask": [ + [ + "water", + { + "band_1": 850, + "band_2": 660, + "threshold": 290.0 + } + ], + [ + "kernel_finite", + {} + ], + [ + "ancillary", + { + "name": "sensor_zn", + "min": 0.0349, + "max": 1000000.0 + } + ], + [ + "ndi", + { + "band_1": 550, + "band_2": 2150, + "min": -1.0, + "max": 0.0 + } + ] + ], + "apply_mask": [ + [ + "water", + { + "band_1": 850, + "band_2": 660, + "threshold": 290.0 + } + ] + ], + "bin_type": "dynamic", + "num_bins": 18, + "ndvi_bin_min": 0.05, + "ndvi_bin_max": 1.0, + "ndvi_perc_min": 10, + "ndvi_perc_max": 95, + "solar_zn_type": "scene" + }, + "glint": { + "type": "hochberg", + "correction_wave": 2150, + "apply_mask": [] + }, + "export": { + "coeffs": false, + "image": true, + "masks": true, + "subset_waves": [], + "output_dir": "E:/is2/yaopu/Flex", + "suffix": "topo_brdf_glint" + }, + "resample": false +} \ No newline at end of file diff --git a/GUI/软件使用说明.md b/GUI/软件使用说明.md new file mode 100644 index 0000000..2a48c89 --- /dev/null +++ b/GUI/软件使用说明.md @@ -0,0 +1,451 @@ +# BRDF校正工具 - 使用说明文档 + +## 一、软件概述 + +BRDF校正工具是一款用于高光谱图像双向反射分布函数（BRDF）校正的统一GUI应用程序，支持两种场景： + +1. **陆地植被BRDF校正**（FlexBRDF算法） + - 适用于陆地植被和地表的BRDF、地形和镜面反射校正 + - 支持AVIRIS、AVIRIS-NG等机载高光谱数据 + +2. **水体BRDF校正**（Ocean BRDF Correction） + - 专用于水体和海洋的BRDF校正 + - 支持多种BRDF模型和不确定性估算 + +--- + +## 二、快速开始 + +### 2.1 启动软件 + +```bash +python brdf_gui.py +``` + +启动后将显示主界面，选择需要使用的模块： +- 点击"进入陆地模块"打开陆地植被BRDF校正 +- 点击"进入水体模块"打开水体BRDF校正 + +--- + +## 三、陆地植被BRDF校正模块（FlexBRDF） + +### 3.1 界面布局 + +模块界面分为左右两部分： +- **左侧**：配置选项卡（文件与波段映射、校正参数、输出设置） +- **右侧**：执行日志和进度显示 + +### 3.2 文件与波段映射页 + +#### 3.2.1 文件夹设置 + +| 参数 | 说明 | 默认值 | +|------|------|--------| +| 输入文件夹 | 反射率文件所在目录（ENVI格式） | 空 | +| 辅助文件夹 | 观测几何数据目录 | 空 | +| CPU核心数 | 并行处理使用的CPU核心数 | 1 | + +**操作步骤**： +1. 点击"浏览..."按钮选择输入文件夹和辅助文件夹 +2. 点击"自动匹配文件"按钮，软件会自动匹配反射率文件（`_rfl`后缀）和观测几何文件（`_obs`后缀） +3. 匹配成功的文件对将显示在列表中 + +#### 3.2.2 文件匹配规则 + +- 反射率文件命名格式：`{prefix}_rfl{ext}` +- 观测几何文件命名格式：`{prefix}_obs{ext}` +- 软件根据`{prefix}`部分自动匹配对应文件对 + +#### 3.2.3 手动添加文件对 + +如需手动指定文件对： +1. 在"手动添加文件对"区域，分别选择反射率文件和观测几何文件 +2. 点击"添加文件对"按钮 + +#### 3.2.4 波段映射设置 + +每个输入文件需要映射以下观测几何参数到辅助文件的波段： + +| 参数名称 | 说明 | 默认波段 | +|----------|------|----------| +| path_length | 路径长度 | 0 | +| sensor_zn | 传感器天顶角 | 1 | +| sensor_az | 传感器方位角 | 2 | +| solar_zn | 太阳天顶角 | 3 | +| solar_az | 太阳方位角 | 4 | +| phase | 相位角 | 5 | +| slope | 坡度 | 6 | +| aspect | 坡向 | 7 | +| cosine_i | 入射角余弦 | 8 | + +**设置方法**： +1. 在文件列表中选择一个文件 +2. 在右侧波段映射表格中设置各参数对应的波段索引（从0开始） +3. 点击"应用当前映射"保存设置 + +### 3.3 校正参数页 + +#### 3.3.1 选择校正类型 + +勾选需要执行的校正类型（可单选或多选）： + +| 校正类型 | 说明 | +|----------|------| +| Topo (地形校正) | 消除地形起伏对反射率的影响 | +| BRDF (双向反射分布函数) | 校正各向异性反射效应 | +| Glint (耀斑校正) | 消除水体或湿润表面的镜面反射 | + +#### 3.3.2 地形校正参数（Topo Correction） + +| 参数 | 选项 | 默认值 | 说明 | +|------|------|--------|------| +| 类型 (type) | mod_minneart, scs+c, cosine, c, scs | mod_minneart | 地形校正算法类型 | +| 拟合类型 (c_fit_type) | nnls, ols, wls | nnls | 系数拟合方法 | + +**掩膜条件设置**： +- **计算掩膜 (calc_mask)**：定义参与计算校正因子的像素条件 +- **应用掩膜 (apply_mask)**：定义应用校正的像素条件 + +掩膜条件类型详见[第5章 掩膜条件说明](#五掩膜条件说明)。 + +#### 3.3.3 BRDF校正参数（BRDF Correction） + +| 参数 | 选项/范围 | 默认值 | 说明 | +|------|-----------|--------|------| +| 算法类型 (type) | flex, universal | flex | BRDF算法类型 | +| 启用分组 (grouped) | 勾选/取消 | 勾选 | 是否按NDVI分组计算 | +| 几何核 (geometric) | li_sparse, li_dense, li_sparse_r, li_dense_r, roujean | li_sparse | 几何光学核函数 | +| 体积核 (volume) | ross_thin, ross_thick, hotspot, roujean | ross_thick | 体积散射核函数 | +| b/r 比率 | 0.1 - 10.0 | 2.5 | 树冠宽度/高度比 | +| h/b 比率 | 0.1 - 5.0 | 2.0 | 树高/树冠宽度比 | +| 采样比例 (sample_perc) | 0.01 - 1.0 | 0.1 | 用于拟合的采样比例 | +| 插值方法 (interp_kind) | linear, nearest, cubic | linear | 插值方法 | +| 分箱类型 (bin_type) | dynamic, fixed | dynamic | NDVI分箱类型 | +| 分箱数量 (num_bins) | 1 - 100 | 18 | NDVI分箱数量 | +| NDVI最小值 (ndvi_bin_min) | 0.0 - 1.0 | 0.05 | 分箱NDVI下限 | +| NDVI最大值 (ndvi_bin_max) | 0.0 - 1.0 | 1.0 | 分箱NDVI上限 | +| NDVI百分比最小值 (ndvi_perc_min) | 0 - 100 | 10 | 有效样本百分比下限 | +| NDVI百分比最大值 (ndvi_perc_max) | 0 - 100 | 95 | 有效样本百分比上限 | +| 太阳天顶角类型 (solar_zn_type) | scene, fixed, mean, custom | scene | 太阳天顶角计算方法 | + +**核函数组合推荐**： +- 植被覆盖区：`li_sparse` + `ross_thick` +- 稀疏植被区：`roujean` + `roujean` +- 浓密植被区：`li_dense` + `ross_thick` + +#### 3.3.4 耀斑校正参数（Glint Correction） + +| 参数 | 选项 | 默认值 | 说明 | +|------|------|--------|------| +| 类型 (type) | mod_minneart, scs+c, cosine, c, scs, nir, mask | nir | 耀斑校正方法 | +| 拟合类型 (c_fit_type) | nnls, ols, wls | nnls | 系数拟合方法 | + +**掩膜条件设置**： +- **计算掩膜 (calc_mask)**：定义参与计算耀斑因子的像素 +- **应用掩膜 (apply_mask)**：定义应用耀斑校正的像素 + +### 3.4 输出设置页 + +| 参数 | 说明 | 默认值 | +|------|------|--------| +| 输出目录 | 校正结果保存路径 | 空 | +| 后缀 (Suffix) | 输出文件名后缀 | topo_brdf_glint | +| 输出系数 (coeffs) | 输出校正系数文件 | 不勾选 | +| 输出图像 (image) | 输出校正后图像 | 勾选 | +| 输出掩膜 (masks) | 输出使用的掩膜 | 不勾选 | + +**输出文件名格式**： +``` +{原文件名}_{后缀}.dat +``` + +### 3.5 操作按钮 + +| 按钮 | 功能 | +|------|------| +| 加载 JSON 配置 | 从JSON文件加载之前保存的配置 | +| 保存 JSON 配置 | 将当前配置保存为JSON文件 | +| 生成配置并运行 | 生成配置并启动BRDF校正处理 | + +--- + +## 四、水体BRDF校正模块（Ocean BRDF） + +### 4.1 输入文件设置 + +| 参数 | 文件格式 | 说明 | +|------|----------|------| +| 高光谱文件 | ENVI (.dat, .bsq, .bip, .bil, .hdr) | 高光谱遥感影像 | +| 角度文件 | ENVI (.dat, .bsq, .bip, .bil, .hdr) | 观测几何角度数据 | +| 水体掩膜 | GeoTIFF (.tif, .tiff) | 水体区域掩膜 | + +### 4.2 BRDF模型参数 + +| 参数 | 选项 | 默认值 | 说明 | +|------|------|--------|------| +| BRDF模型 | L11, M02, M02SeaDAS, O25 | L11 | BRDF校正模型 | +| 处理块大小 | 1024 - 16384 | 4096 | 光谱维分块大小 | +| 空间块大小 | 256 - 2048 | 512 | 空间维分块大小 | +| 计算不确定性 | 勾选/取消 | 取消 | 是否计算BRDF不确定度 | + +**BRDF模型说明**： +- **L11**：Lee et al. (2011) 模型，适用于开阔大洋 +- **M02**：Morel et al. (2002) 模型，经典水体BRDF模型 +- **M02SeaDAS**：SeaDAS实现的M02模型 +- **O25**：最新优化模型 + +### 4.3 输出变量选择 + +支持多选（按住Ctrl键选择多个）： + +| 变量 | 说明 | +|------|------| +| Rw_brdf | BRDF校正后的水体反射率 | +| rho_ex_w | 水体反射率 (nrrs * π) | +| nrrs | 归一化遥感反射率 | +| C_brdf | BRDF校正因子 | +| brdf_unc | BRDF不确定度（需勾选计算不确定性） | +| nrrs_unc | nrrs不确定度（需勾选计算不确定性） | + +**默认选中**：Rw_brdf, nrrs, C_brdf + +### 4.4 输出设置 + +| 参数 | 选项 | 默认值 | +|------|------|--------| +| 输出路径 | 文件前缀路径 | 空 | +| 输出格式 | ENVI | ENVI | + +### 4.5 运行 + +点击"运行水体BRDF校正"按钮开始处理。 + +--- + +## 五、掩膜条件说明 + +掩膜用于定义参与计算或应用校正的像素范围。所有掩膜条件取交集（同时满足）。 + +### 5.1 支持的掩膜类型 + +| 类型 | 说明 | 参数 | +|------|------|------| +| **ndi** | 归一化差值指数 | band_1, band_2, min, max | +| **ancillary** | 辅助数据掩膜 | name, min, max | +| **neon_edge** | NEON数据边缘掩膜 | radius | +| **kernel_finite** | 核函数有限性掩膜 | 无参数 | +| **water** | 水体掩膜 | band_1, band_2, threshold | +| **external** | 外部掩膜文件 | file_path, class | +| **band** | 单波段阈值 | band, min, max | + +### 5.2 各类型详细参数 + +#### ndi（归一化差值指数） + +``` +ndi = (band_1 - band_2) / (band_1 + band_2) +保留范围: [min, max] 内的像素 +``` + +| 参数 | 类型 | 默认值 | 说明 | +|------|------|--------|------| +| band_1 | 整数 | 550 | 第一个波段波长(nm) | +| band_2 | 整数 | 2150 | 第二个波段波长(nm) | +| min | 浮点数 | -1.0 | NDI最小值 | +| max | 浮点数 | 0.0 | NDI最大值 | + +**典型应用**： +- NDVI植被掩膜：band_1=850, band_2=680, min=0.2, max=1.0 +- NDWI水体掩膜：band_1=860, band_2=1240, min=-1.0, max=0.0 + +#### ancillary（辅助数据掩膜） + +基于观测几何参数创建掩膜。 + +| 参数 | 类型 | 默认值 | 可选值 | +|------|------|--------|--------| +| name | 字符串 | slope | slope, aspect, cosine_i, sensor_zn, sensor_az, solar_zn, solar_az, phase, path_length, utc_time | +| min | 浮点数 | 0.0 | - | +| max | 浮点数 | 1.0 | - | + +**典型应用**： +- 坡度掩膜：只处理坡度小于30度的区域 + - name=slope, min=0, max=30 +- 太阳天顶角掩膜：排除大太阳天顶角 + - name=solar_zn, min=0, max=60 + +#### neon_edge（NEON边缘掩膜） + +用于消除NEON高光谱数据图像边缘的噪声区域。 + +| 参数 | 类型 | 默认值 | 说明 | +|------|------|--------|------| +| radius | 整数 | 1 | 边缘裁剪像素数 | + +#### kernel_finite（核函数有限性掩膜） + +自动排除会导致核函数计算溢出或不稳定的像素。无需参数。 + +#### water（水体掩膜） + +使用特定波段组合识别水体。 + +| 参数 | 类型 | 默认值 | 说明 | +|------|------|--------|------| +| band_1 | 整数 | 550 | 第一个波段波长(nm) | +| band_2 | 整数 | 850 | 第二个波段波长(nm) | +| threshold | 浮点数 | 0.0 | 水体识别阈值 | + +#### external（外部掩膜文件） + +使用外部分类或掩膜文件。 + +| 参数 | 类型 | 默认值 | 说明 | +|------|------|--------|------| +| file_path | 字符串 | 空 | 外部掩膜文件路径 | +| class | 整数 | 1 | 要保留的类别值 | + +#### band（单波段阈值） + +基于单个波段的值创建掩膜。 + +| 参数 | 类型 | 默认值 | 说明 | +|------|------|--------|------| +| band | 整数 | 650 | 波段波长(nm) | +| min | 浮点数 | 0.0 | 最小值 | +| max | 浮点数 | 1.0 | 最大值 | + +### 5.3 掩膜条件设置步骤 + +1. 在校正参数页找到掩膜设置区域 +2. 点击"+ 添加条件"按钮 +3. 选择掩膜类型 +4. 设置对应参数值 +5. 可添加多个条件（取交集） +6. 点击"删除"按钮移除不需要的条件 + +--- + +## 六、配置文件的保存与加载 + +### 6.1 保存配置 + +1. 完成所有参数设置 +2. 点击"保存 JSON 配置"按钮 +3. 选择保存路径和文件名 +4. 配置将以JSON格式保存，包含： + - 输入文件列表 + - 波段映射关系 + - 校正参数设置 + - 掩膜条件 + - 输出设置 + +### 6.2 加载配置 + +1. 点击"加载 JSON 配置"按钮 +2. 选择之前保存的JSON配置文件 +3. 软件将自动恢复所有设置 + +**注意**： +- 加载配置后会自动匹配文件对，请确保输入文件和辅助文件路径仍然有效 +- 如果文件路径已改变，需要重新手动设置文件夹路径 + +--- + +## 七、典型工作流程 + +### 7.1 陆地植被BRDF校正流程 + +1. **启动软件** → 选择"陆地植被BRDF校正" +2. **设置文件夹** → 选择输入文件夹和辅助文件夹 +3. **匹配文件** → 点击"自动匹配文件" +4. **检查映射** → 确认波段映射设置正确 +5. **选择校正** → 勾选需要的校正类型（Topo/BRDF/Glint） +6. **设置参数** → 根据数据特点调整校正参数 +7. **设置掩膜** → 添加合适的掩膜条件（如NDVI>0.2） +8. **设置输出** → 指定输出目录和后缀 +9. **保存配置** → （可选）保存当前配置 +10. **运行** → 点击"生成配置并运行" + +### 7.2 水体BRDF校正流程 + +1. **启动软件** → 选择"水体BRDF校正" +2. **选择文件** → 依次选择高光谱文件、角度文件、水体掩膜 +3. **选择模型** → 选择适合的BRDF模型（通常用L11） +4. **选择变量** → 勾选需要的输出变量 +5. **调整参数** → （可选）调整块大小以优化性能 +6. **设置输出** → 指定输出文件前缀路径 +7. **运行** → 点击"运行水体BRDF校正" + +--- + +## 八、常见问题与解决方案 + +### 8.1 文件匹配问题 + +**问题**：自动匹配找不到文件对 + +**解决**： +- 检查文件命名是否符合 `{prefix}_rfl` 和 `{prefix}_obs` 格式 +- 确保输入文件夹和辅助文件夹路径正确 +- 使用手动添加功能指定文件对 + +### 8.2 波段映射问题 + +**问题**：提示波段索引超出范围 + +**解决**： +- 检查辅助文件的实际波段数 +- 确认波段索引从0开始计数 +- 验证辅助文件格式正确 + +### 8.3 内存不足问题 + +**问题**：处理大文件时内存溢出 + +**解决**： +- 减少CPU核心数设置 +- 对于水体模块，减小处理块大小（chunk_size） +- 分块处理（FlexBRDF会自动分块） + +### 8.4 掩膜问题 + +**问题**：校正后结果为空或过少 + +**解决**： +- 检查掩膜条件是否过于严格 +- 验证NDVI或波段阈值设置是否合理 +- 查看掩膜类型是否与数据匹配 + +--- + +## 九、技术参数参考 + +### 9.1 支持的文件格式 + +| 类型 | 格式 | +|------|------| +| 输入反射率 | ENVI标准格式 (.dat + .hdr) | +| 观测几何 | ENVI标准格式 (.dat + .hdr) | +| 水体掩膜 | GeoTIFF (.tif) | +| 输出结果 | ENVI标准格式 | +| 配置文件 | JSON (.json) | + +### 9.2 系统要求 + +- **操作系统**：Windows 10/11, Linux, macOS +- **Python版本**：3.8+ +- **内存**：建议8GB以上（处理大图像需要更多） +- **依赖库**：PyQt5, numpy, scipy等（详见requirements.txt） + +--- + +## 十、联系与支持 + +如有问题或建议，请查看项目文档或联系开发团队。 + +--- + +*文档版本：v1.0* +*最后更新：2026年4月* diff --git a/angle/angle_compute.py b/angle/angle_compute.py new file mode 100644 index 0000000..201826a --- /dev/null +++ b/angle/angle_compute.py @@ -0,0 +1,1417 @@ +""" +无人机高光谱线阵相机角度计算程序 + +计算每个像素的太阳天顶角、太阳方位角、传感器天顶角、传感器方位角、相对方位角 +以南向北为正向 +""" + +import os +import sys +import numpy as np +import pandas as pd +from datetime import datetime, timedelta +import time +import re +from typing import Dict, List, Tuple, Optional +from pyproj import Transformer +import pvlib.solarposition as solarposition +from pathlib import Path +import spectral +import warnings + +warnings.filterwarnings('ignore') + +try: + from tqdm import tqdm + HAS_TQDM = True +except ImportError: + HAS_TQDM = False + print("警告: 未安装tqdm库，无法显示进度条。建议安装: pip install tqdm") + + +class UAVAngleCalculator: + """无人机高光谱线阵相机角度计算器""" + + def __init__(self, config: Dict): + """ + 初始化计算器 + + Args: + config: 配置字典 + """ + self.config = config + self.data_dir = Path(config.get('data_dir', './Data')) + # times 根目录：内部包含多个“航带名”子目录（或兼容旧结构：直接放 *.times） + self.times_root_dir = Path(config.get('times_root_dir', config.get('times_dir', config.get('data_dir', './Data')))) + self.base_height = config.get('base_height', 254) # 基准高度 + self.output_dir = Path(config.get('output_dir', './Output')) + + # 创建输出目录 + self.output_dir.mkdir(exist_ok=True) + + # 初始化坐标变换器 + self.utm_crs = "EPSG:32651" # UTM Zone 51N + self.wgs84_crs = "EPSG:4326" + self.transformer = Transformer.from_crs(self.wgs84_crs, self.utm_crs, always_xy=True) + + # 缓存 + self.gps_data = None + self.times_data = {} + self.image_info = None + self.image_data = None + + def read_bip_image(self, bip_file: Path) -> Tuple[np.ndarray, Dict]: + """ + 读取BIP格式图像文件 + + Args: + bip_file: BIP文件路径 + + Returns: + 图像数据数组和图像信息字典 + """ + print(f"读取BIP图像文件: {bip_file}") + + # 使用spectral库读取ENVI文件 + print("使用spectral库读取图像数据...") + start_time = time.time() + + try: + # spectral会自动读取HDR文件并解析元数据 + image_data = spectral.open_image(str(bip_file)) + + # 获取图像数组 + image_array = image_data.load() + + # 解析图像信息 + image_info = self._extract_spectral_info(image_data) + + read_time = time.time() - start_time + print(f"图像数据读取完成，耗时: {read_time:.2f} 秒") + print(f"图像形状: {image_array.shape}") + + return image_array, image_info + + except Exception as e: + print(f"使用spectral读取失败: {e}") + print("回退到手动解析方法...") + + # 回退到原来的方法 + return self._read_bip_manual(bip_file) + + def _extract_band_names_from_hdr(self, hdr_file: Path, band_count: int) -> List[str]: + """ + 从HDR文件中提取band names，若不存在则生成默认名称 + """ + if not hdr_file.exists(): + return [f'Band {i + 1}' for i in range(band_count)] + + try: + with open(hdr_file, 'r', encoding='utf-8', errors='ignore') as f: + lines = f.read().splitlines() + except Exception: + return [f'Band {i + 1}' for i in range(band_count)] + + names = [] + in_block = False + for line in lines: + line_stripped = line.strip() + if not in_block: + if line_stripped.lower().startswith('band names'): + in_block = True + if '{' in line_stripped: + line_stripped = line_stripped.split('{', 1)[1].strip() + else: + line_stripped = '' + if in_block: + if '}' in line_stripped: + line_stripped = line_stripped.split('}', 1)[0].strip() + in_block = False + if line_stripped: + parts = [p.strip().strip("'\"") for p in line_stripped.split(',') if p.strip()] + names.extend(parts) + + if not in_block and names: + break + + if len(names) < band_count: + names = [f'Band {i + 1}' for i in range(band_count)] + else: + names = names[:band_count] + + return names + + def _extract_spectral_info(self, spectral_img) -> Dict: + """ + 从spectral图像对象中提取图像信息 + + Args: + spectral_img: spectral图像对象 + + Returns: + 图像信息字典 + """ + info = {} + + # 基本信息 + info['samples'] = spectral_img.ncols + info['lines'] = spectral_img.nrows + info['bands'] = spectral_img.nbands + info['data_type'] = spectral_img.dtype + + # 数据类型字节数 + dtype_map = { + np.uint8: 1, + np.int16: 2, + np.float32: 4, + np.float64: 8, + np.uint16: 2, + np.uint32: 4, + } + info['data_type_bytes'] = dtype_map.get(spectral_img.dtype, 4) + + # 交织方式 + info['interleave'] = spectral_img.interleave + + # 元数据 + metadata = spectral_img.metadata + print(f"元数据键: {list(metadata.keys())}") + if 'map info' in metadata: + # 解析地图信息 + map_info = metadata['map info'] + print(f"map info 类型: {type(map_info)}") + print(f"map info 内容: {map_info}") + + # 如果是字典格式（spectral库已解析） + if isinstance(map_info, dict): + info['projection'] = map_info.get('projection', 'UTM') + info['x_ul'] = map_info.get('x_ul', map_info.get('ul_x', 0)) + info['y_ul'] = map_info.get('y_ul', map_info.get('ul_y', 0)) + info['pixel_x'] = map_info.get('pixel_x', map_info.get('pixel_size_x', 1)) + info['pixel_y'] = map_info.get('pixel_y', map_info.get('pixel_size_y', 1)) + info['utm_zone'] = str(map_info.get('zone', 51)) + info['hemisphere'] = map_info.get('hemisphere', 'North') + info['datum'] = map_info.get('datum', 'WGS-84') + # 为兼容性设置x_start和y_start + info['x_start'] = info['x_ul'] + info['y_start'] = info['y_ul'] + else: + # 如果是字符串格式，手动解析 + map_info_str = str(map_info) + # 移除花括号 + map_info_str = map_info_str.strip('{}') + map_parts = map_info_str.split(',') + if len(map_parts) >= 9: + info['projection'] = map_parts[0].strip() + + # ENVI map info 格式： + # [投影, ref_col, ref_row, ref_x, ref_y, pixel_x, pixel_y, zone, hemisphere, ...] + ref_col = float(map_parts[1].strip().strip("'\"")) + ref_row = float(map_parts[2].strip().strip("'\"")) + ref_x = float(map_parts[3].strip().strip("'\"")) + ref_y = float(map_parts[4].strip().strip("'\"")) + pixel_x = float(map_parts[5].strip().strip("'\"")) + pixel_y = float(map_parts[6].strip().strip("'\"")) + info['utm_zone'] = map_parts[7].strip().strip("'\"") + info['hemisphere'] = map_parts[8].strip().strip("'\"") + info['datum'] = map_parts[9].strip().strip("'\"") if len(map_parts) > 9 else 'WGS-84' + + # 从参考像元回推出左上角像元左上角的坐标 + x_ul = ref_x - (ref_col - 0.5) * pixel_x + y_ul = ref_y + (ref_row - 0.5) * pixel_y + + # 保存解析结果 + info['ref_col'] = ref_col + info['ref_row'] = ref_row + info['ref_x'] = ref_x + info['ref_y'] = ref_y + info['pixel_x'] = pixel_x + info['pixel_y'] = pixel_y + info['x_ul'] = x_ul + info['y_ul'] = y_ul + + # 构建GDAL风格的仿射变换参数 + # [x_ul, pixel_x, 0, y_ul, 0, -pixel_y] + info['geo_transform'] = [ + info['x_ul'], # X坐标偏移 + info['pixel_x'], # X像素大小 + 0.0, # X旋转（通常为0） + info['y_ul'], # Y坐标偏移 + 0.0, # Y旋转（通常为0） + -info['pixel_y'] # Y像素大小（负号表示Y方向向下） + ] + + print(f"图像信息: {info['samples']}x{info['lines']}x{info['bands']}像素") + print(f"数据类型: {info['data_type']} ({info['data_type_bytes']}字节)") + if 'utm_zone' in info: + print(f"地理参考: UTM Zone {info['utm_zone']}{info.get('hemisphere', 'N')}") + print(f"左上角坐标: ({info['x_ul']:.2f}, {info['y_ul']:.2f})") + print(f"像素大小: {info['pixel_x']}x{info['pixel_y']}米") + + return info + + def _read_bip_manual(self, bip_file: Path) -> Tuple[np.ndarray, Dict]: + """ + 手动读取BIP格式图像文件（回退方法） + + Args: + bip_file: BIP文件路径 + + Returns: + 图像数据数组和图像信息字典 + """ + print(f"手动读取BIP图像文件: {bip_file}") + + # 读取HDR文件获取图像信息 + hdr_file = bip_file.with_suffix('.hdr') + image_info = self._parse_envi_header(hdr_file) + + # 计算文件大小验证 + expected_size = (image_info['samples'] * image_info['lines'] * + image_info['bands'] * image_info['data_type_bytes']) + actual_size = bip_file.stat().st_size + + if actual_size != expected_size: + raise ValueError(f"文件大小不匹配。期望: {expected_size} bytes, 实际: {actual_size} bytes") + + # 读取BIP数据 + print("读取图像数据...") + start_time = time.time() + + # 根据数据类型设置numpy dtype + dtype_map = { + 1: np.uint8, # Byte + 2: np.int16, # Integer (2 bytes) + 4: np.float32, # Float (4 bytes) + 5: np.float64, # Double (8 bytes) + 12: np.uint16, # Unsigned integer (2 bytes) + 13: np.uint32, # Unsigned long (4 bytes) + } + + dtype = dtype_map.get(image_info['data_type'], np.float32) + + # 使用内存映射读取大文件 + image_data = np.memmap(str(bip_file), dtype=dtype, mode='r', + shape=(image_info['lines'], image_info['samples'], image_info['bands']), + order='C') # BIP格式是按像素交织的 + + read_time = time.time() - start_time + print(f"图像数据读取完成，耗时: {read_time:.2f} 秒") + return image_data, image_info + + def _parse_envi_header(self, hdr_file: Path) -> Dict: + """ + 解析ENVI格式的HDR头文件 + + Args: + hdr_file: HDR文件路径 + + Returns: + 图像信息字典 + """ + info = {} + + with open(hdr_file, 'r', encoding='utf-8', errors='ignore') as f: + content = f.read() + + # 解析基本参数 + lines = content.split('\n') + for line in lines: + line = line.strip() + if '=' in line: + key, value = line.split('=', 1) + key = key.strip().lower() + value = value.strip() + + # 移除花括号和等号 + value = value.strip('{}=') + + if key == 'samples': + info['samples'] = int(value) + elif key == 'lines': + info['lines'] = int(value) + elif key == 'bands': + info['bands'] = int(value) + elif key == 'data type': + info['data_type'] = int(value) + elif key == 'byte order': + info['byte_order'] = int(value) + elif key == 'interleave': + info['interleave'] = value.lower() + elif key == 'data ignore value': + info['ignore_value'] = float(value) if '.' in value else int(value) + elif key == 'map info': + # 解析地图信息 + map_parts = value.split(',') + if len(map_parts) >= 9: + info['projection'] = map_parts[0].strip() + + # ENVI map info 格式： + # [投影, ref_col, ref_row, ref_x, ref_y, pixel_x, pixel_y, zone, hemisphere, ...] + # ref_col/row: 参考像元列/行号 (ENVI 1基, 通常为像元中心) + # ref_x/y: 该参考像元的地图坐标 + ref_col = float(map_parts[1].strip()) + ref_row = float(map_parts[2].strip()) + ref_x = float(map_parts[3].strip()) + ref_y = float(map_parts[4].strip()) + pixel_x = float(map_parts[5].strip()) + pixel_y = float(map_parts[6].strip()) + info['utm_zone'] = map_parts[7].strip() + info['hemisphere'] = map_parts[8].strip() + info['datum'] = map_parts[9].strip() if len(map_parts) > 9 else 'WGS-84' + + # 从参考像元回推出左上角像元左上角的坐标 + # ENVI参考像元为1基且为中心坐标，GDAL原点为(0,0)像元左上角 + x_ul = ref_x - (ref_col - 0.5) * pixel_x # 向西偏移得到左上角X + y_ul = ref_y + (ref_row - 0.5) * pixel_y # 向北偏移得到左上角Y + + # 保存解析结果 + info['ref_col'] = ref_col + info['ref_row'] = ref_row + info['ref_x'] = ref_x + info['ref_y'] = ref_y + info['pixel_x'] = pixel_x + info['pixel_y'] = pixel_y + info['x_ul'] = x_ul + info['y_ul'] = y_ul + + # 构建GDAL风格的仿射变换参数 + # [x_ul, pixel_x, 0, y_ul, 0, -pixel_y] + info['geo_transform'] = [ + x_ul, # X坐标偏移 (左上角像元左上角X) + pixel_x, # X像素大小 + 0.0, # X旋转（通常为0） + y_ul, # Y坐标偏移 (左上角像元左上角Y) + 0.0, # Y旋转（通常为0） + -pixel_y # Y像素大小（负号表示Y方向向下） + ] + + # 根据data type设置字节数 + data_type_bytes = { + 1: 1, # Byte + 2: 2, # Integer + 4: 4, # Float + 5: 8, # Double + 12: 2, # Unsigned integer + 13: 4, # Unsigned long + } + info['data_type_bytes'] = data_type_bytes.get(info['data_type'], 4) + + print(f"图像信息: {info['samples']}x{info['lines']}x{info['bands']}像素") + print(f"数据类型: {info['data_type']} ({info['data_type_bytes']}字节)") + print(f"地理参考: UTM Zone {info['utm_zone']}{info['hemisphere']}") + print(f"左上角坐标: ({info['x_ul']:.2f}, {info['y_ul']:.2f})") + print(f"像素大小: {info['pixel_x']}x{info['pixel_y']}米") + + return info + + def extract_pixel_info(self, image_data: np.ndarray, image_info: Dict) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + 从图像数据中提取像素信息 + + Args: + image_data: BIP图像数据 + image_info: 图像信息 + + Returns: + 列号、行号、航带号数组 + """ + print("提取像素几何信息...") + + # BIP格式：bands顺序为RGB + 列号 + 行号 + 航带号 + # 波段3：列号，波段4：行号，波段5：航带号 + col_indices = image_data[:, :, 3].astype(np.int32) # 列号 + row_indices = image_data[:, :, 4].astype(np.int32) # 行号 + swath_ids = image_data[:, :, 5].astype(np.int32) # 航带号 + + # 排除无效像素（所有波段值均为0的像素为填充像素） + valid_mask = np.any(image_data != 0, axis=2) + + print(f"总像素数: {col_indices.size}") + print(f"有效像素数: {np.sum(valid_mask)}") + print(f"航带数量: {len(np.unique(swath_ids[valid_mask]))}") + + return col_indices, row_indices, swath_ids + + def load_gps_data(self, gps_file: Path) -> pd.DataFrame: + """ + 加载GPS/IMU数据 + + Args: + gps_file: GPS数据文件路径 + + Returns: + GPS数据DataFrame + """ + print(f"加载GPS数据: {gps_file}") + + # 读取GPS数据（空格分隔的文本文件） + gps_data = pd.read_csv(gps_file, sep=r'\s+', header=None, engine='python', + names=['Date', 'Time', 'roll', 'pitch', 'yaw', + 'Lon', 'Lat', 'Height']) + + # 合并日期时间并转换为datetime + gps_data['datetime'] = pd.to_datetime(gps_data['Date'] + ' ' + gps_data['Time']) + gps_data = gps_data.drop(columns=['Date', 'Time']) + + # GPS数据已经是UTC时间，保持为普通datetime格式（无时区信息） + # 以便与times时间戳格式匹配 + + # 添加UTM坐标 + utm_x, utm_y = self.transformer.transform(gps_data['Lon'].values, gps_data['Lat'].values) + gps_data['X_utm'] = utm_x + gps_data['Y_utm'] = utm_y + + print(f"GPS数据加载完成，共 {len(gps_data)} 个记录") + print(f"时间范围: {gps_data['datetime'].min()} 到 {gps_data['datetime'].max()}") + print(f"位置范围: 经度 {gps_data['Lon'].min():.6f}°-{gps_data['Lon'].max():.6f}°, " + f"纬度 {gps_data['Lat'].min():.6f}°-{gps_data['Lat'].max():.6f}°") + + return gps_data + + def load_times_data(self, times_file: Path) -> np.ndarray: + """ + 加载时间戳数据 + + Args: + times_file: times文件路径 + + Returns: + 时间戳数组（周内秒） + """ + print(f"加载时间戳数据: {times_file}") + + # 读取时间戳数据 + times_data = np.loadtxt(times_file) + + print(f"时间戳数据加载完成，共 {len(times_data)} 个记录") + print(f"时间范围: {times_data[0]:.6f} - {times_data[-1]:.6f} 秒") + + return times_data + + def convert_week_seconds_to_datetime(self, week_seconds: float, base_date: datetime) -> datetime: + """ + 将GPS周内秒转换为绝对时间 + + Args: + week_seconds: GPS周内秒 + base_date: 基准日期 + + Returns: + 绝对时间 + """ + # 找到基准日期所在周的周日0点 + days_since_sunday = base_date.weekday() + 1 if base_date.weekday() != 6 else 0 + week_start = base_date - timedelta(days=days_since_sunday) + week_start = datetime(week_start.year, week_start.month, week_start.day, 0, 0, 0) + + # 周日起点 + 周内秒 = 绝对时间 + absolute_time = week_start + timedelta(seconds=week_seconds) + + return absolute_time + + + + + + def get_pixel_coordinates(self, col_idx: int, row_idx: int, swath_id: int, + image_info: Dict) -> Tuple[float, float]: + """ + 根据像素行列号计算像素的UTM坐标 + + 使用从ENVI头文件中读取的GDAL仿射变换参数进行坐标变换计算 + + Args: + col_idx: 当前图像中的列号（像素列索引） + row_idx: 当前图像中的行号（像素行索引） + swath_id: 航带号 + image_info: 图像信息（包含geo_transform） + + Returns: + UTM坐标 (X, Y) + """ + # 从头文件中获取仿射变换参数 + geo_transform = image_info['geo_transform'] + + # 使用GDAL标准的像素到地理坐标变换公式 + pixel_x = geo_transform[0] + col_idx * geo_transform[1] + row_idx * geo_transform[2] + pixel_y = geo_transform[3] + col_idx * geo_transform[4] + row_idx * geo_transform[5] + + return pixel_x, pixel_y + + def _precompute_pixel_coordinates(self, image_info: Dict, valid_pixels: np.ndarray) -> np.ndarray: + """ + 预计算所有有效像素的UTM坐标（向量化）- 参照raster_to_points.py的实现 + + Args: + image_info: 图像信息，包含geo_transform + valid_pixels: 有效像素位置数组 (N, 2) - [row, col] + + Returns: + UTM坐标数组 (N, 2) - [x, y] + """ + if len(valid_pixels) == 0: + return np.empty((0, 2), dtype=np.float64) + + geo_transform = image_info['geo_transform'] + rows, cols = valid_pixels[:, 0], valid_pixels[:, 1] + + # 使用GDAL标准的像素到地理坐标变换公式 + # pixel_x = geo_transform[0] + col * geo_transform[1] + row * geo_transform[2] + # pixel_y = geo_transform[3] + col * geo_transform[4] + row * geo_transform[5] + pixel_x = geo_transform[0] + cols * geo_transform[1] + rows * geo_transform[2] + pixel_y = geo_transform[3] + cols * geo_transform[4] + rows * geo_transform[5] + + return np.column_stack([pixel_x, pixel_y]) + + def _process_pixels_chunked_optimized(self, valid_pixels: np.ndarray, col_indices: np.ndarray, + row_indices: np.ndarray, swath_ids: np.ndarray, + times_data_cache: Dict, base_date: datetime, + image_info: Dict) -> np.ndarray: + """ + 优化的处理版本 - 使用主数据表和最近邻GPS匹配，实现高性能向量化计算 + """ + total_start_time = time.time() + + angle_results = np.full((image_info['lines'], image_info['samples'], 5), np.nan, dtype=np.float32) + total_pixels = len(valid_pixels) + + print(f"开始构建主数据表，共 {total_pixels} 个有效像素") + build_start = time.time() + + # 构建主数据表 (Master Table) + master_table = self._build_master_table(valid_pixels, col_indices, row_indices, swath_ids, + times_data_cache, base_date, image_info) + + build_time = time.time() - build_start + print(".2f") + if len(master_table) == 0: + print("警告: 没有有效的像素数据") + return angle_results + + print(f"主数据表构建完成，包含 {len(master_table)} 条记录") + + # 一次性向量化角度计算 + print("开始向量化角度计算...") + angle_calc_start = time.time() + + if HAS_TQDM: + # 使用tqdm显示进度，但由于是向量化计算，我们用一个简单的进度指示 + print(f"正在计算 {len(master_table)} 个像素的角度...") + angles = self._calculate_angles_vectorized_from_master_table(master_table, base_date) + else: + angles = self._calculate_angles_vectorized_from_master_table(master_table, base_date) + + angle_calc_time = time.time() - angle_calc_start + print(".2f") + # 将结果写回图像矩阵 + print("将结果写回图像矩阵...") + write_back_start = time.time() + + if HAS_TQDM: + from tqdm import tqdm + for i, (idx, row) in enumerate(tqdm(master_table.iterrows(), desc="写回结果", unit="pixel")): + img_row, img_col = int(row['img_row']), int(row['img_col']) + angle_results[img_row, img_col] = angles[i] + else: + for i, (idx, row) in enumerate(master_table.iterrows()): + img_row, img_col = int(row['img_row']), int(row['img_col']) + angle_results[img_row, img_col] = angles[i] + + write_back_time = time.time() - write_back_start + print(".2f") + total_time = time.time() - total_start_time + print(".2f") + print("角度计算完成") + + return angle_results + + def _build_master_table(self, valid_pixels: np.ndarray, col_indices: np.ndarray, + row_indices: np.ndarray, swath_ids: np.ndarray, + times_data_cache: Dict, base_date: datetime, + image_info: Dict) -> pd.DataFrame: + """ + 构建主数据表，包含像素UTM坐标、行列号、航带号和匹配的GPS数据 + + Args: + valid_pixels: 有效像素位置数组 (N, 2) - [row, col] + col_indices: 列号数组 + row_indices: 行号数组 + swath_ids: 航带号数组 + times_data_cache: times数据缓存 + base_date: 基准日期 + image_info: 图像信息 + + Returns: + 主数据表DataFrame + """ + # 预计算所有像素的UTM坐标 + pixel_coords = self._precompute_pixel_coordinates(image_info, valid_pixels) + + # 构建基础数据 + master_data = [] + + for idx, (img_row, img_col) in enumerate(valid_pixels): + img_row = int(img_row) # 转换为 Python int + img_col = int(img_col) # 转换为 Python int + # 使用 .item() 安全地提取标量值 + orig_col = int(col_indices[img_row, img_col].item()) + orig_row = int(row_indices[img_row, img_col].item()) + swath_id = int(swath_ids[img_row, img_col].item()) + + # 检查航带有对应的times数据 + if swath_id not in times_data_cache: + continue + + times_data = times_data_cache[swath_id] + + # 检查行号是否有效 + if orig_row >= len(times_data) - 1: + continue + + # 计算像素UTM坐标 + pixel_x, pixel_y = pixel_coords[idx] + + # 计算该行的采集时间（周内秒）- 使用行开始时间 + t_center = times_data[orig_row] + + # 转换为绝对时间 + absolute_time = self.convert_week_seconds_to_datetime(t_center, base_date) + + master_data.append({ + 'img_row': img_row, + 'img_col': img_col, + 'pixel_x': pixel_x, + 'pixel_y': pixel_y, + 'col_idx': orig_col, + 'row_idx': orig_row, + 'swath_id': swath_id, + 'timestamp': absolute_time, + 'week_seconds': t_center + }) + + # 转换为DataFrame + master_df = pd.DataFrame(master_data) + + if len(master_df) == 0: + return master_df + + # 使用最近邻匹配找到对应的GPS数据 + print("进行GPS数据最近邻匹配...") + + # 为了提高效率，我们对每个唯一的 (swath_id, row_idx) 组合进行一次GPS匹配 + # 因为同一行的所有像素具有相同的时间戳 + unique_time_rows = master_df[['swath_id', 'row_idx', 'timestamp']].drop_duplicates() + + # GPS最近邻匹配 + gps_matched = pd.merge_asof( + unique_time_rows.sort_values('timestamp'), + self.gps_data.sort_values('datetime')[['datetime', 'X_utm', 'Y_utm', 'Height', 'Lon', 'Lat']], + left_on='timestamp', + right_on='datetime', + direction='nearest' + ) + + # 重命名GPS列 + gps_matched = gps_matched.rename(columns={ + 'X_utm': 'sensor_x', + 'Y_utm': 'sensor_y', + 'Height': 'sensor_height', + 'Lon': 'sensor_lon', + 'Lat': 'sensor_lat' + }) + + # 计算相对高程 + gps_matched['relative_height'] = gps_matched['sensor_height'] - self.base_height + + # 合并回主数据表 + master_df = master_df.merge( + gps_matched[['swath_id', 'row_idx', 'sensor_x', 'sensor_y', 'sensor_height', + 'relative_height', 'sensor_lon', 'sensor_lat']], + on=['swath_id', 'row_idx'], + how='left' + ) + + # 添加本地时间（从UTC转换为北京时间） + master_df['local_time'] = master_df['timestamp'].dt.tz_localize('UTC').dt.tz_convert('Asia/Shanghai') + + print(f"GPS匹配完成，共匹配 {len(master_df)} 条记录") + + return master_df + + def _calculate_angles_vectorized_from_master_table(self, master_chunk: pd.DataFrame, + base_date: datetime) -> np.ndarray: + """ + 从主数据表向量化计算角度 + + Args: + master_chunk: 主数据表的一个数据块 + base_date: 基准日期 + + Returns: + 角度数组 (N, 5) - [SZA, SAA, VZA, VAA, RAA] + """ + n_pixels = len(master_chunk) + angles = np.empty((n_pixels, 3), dtype=np.float32) + + # 批量计算太阳角度 + solar_angles = self._calculate_solar_angles_batch( + master_chunk['local_time'].values, + master_chunk['sensor_lon'].values, + master_chunk['sensor_lat'].values + ) + + # 向量化计算传感器角度 + # 提取所有必要的坐标数据 + sensor_x = master_chunk['sensor_x'].values + sensor_y = master_chunk['sensor_y'].values + sensor_z = master_chunk['relative_height'].values + pixel_x = master_chunk['pixel_x'].values + pixel_y = master_chunk['pixel_y'].values + + # 批量计算传感器天顶角 (VZA) + horizontal_dist = np.sqrt((sensor_x - pixel_x) ** 2 + (sensor_y - pixel_y) ** 2) + with np.errstate(divide='ignore', invalid='ignore'): + vza = 90.0 - np.degrees(np.arctan(sensor_z / horizontal_dist)) + # 约束到 0~90（从天顶的夹角） + vza = np.clip(vza, 0.0, 90.0) + + # 批量计算观测方位角 (VAA) + dx = sensor_x - pixel_x # 东方向差 + dy = sensor_y - pixel_y # 北方向差 + angle_rad = np.arctan2(dx, dy) + vaa = np.degrees(angle_rad) % 360.0 + + # 批量计算相对方位角 (RAA) + saa_array = solar_angles[:, 1] # 太阳方位角数组 + diff = np.abs(saa_array - vaa) + raa = np.minimum(diff, 360 - diff) + + # 组合所有角度 + angles = np.column_stack([solar_angles[:, 0], solar_angles[:, 1], vza, vaa, raa]) # [SZA, SAA, VZA, VAA, RAA] + + return angles + + def _calculate_solar_angles_batch(self, datetime_array, lon_array, lat_array) -> np.ndarray: + """ + 批量计算太阳角度（完全向量化版本） + + Args: + datetime_array: 时间戳数组 + lon_array: 经度数组 + lat_array: 纬度数组 + + Returns: + 太阳角度数组 (N, 2) - [SZA, SAA] + """ + n_points = len(datetime_array) + + # 创建DatetimeIndex数组，转换为北京时间 + datetime_idx = pd.DatetimeIndex(datetime_array).tz_localize('UTC').tz_convert('Asia/Shanghai') + + # 批量计算时角方程 + equation_of_time = solarposition.equation_of_time_spencer71(datetime_idx.dayofyear) + + # 批量计算时角 + hour_angle = solarposition.hour_angle(datetime_idx, lon_array, equation_of_time) + + # 批量计算太阳赤纬 + declination = solarposition.declination_cooper69(datetime_idx.dayofyear) + + # 批量计算太阳天顶角 + lat_rad = np.radians(lat_array) + hour_angle_rad = np.radians(hour_angle) + declination_rad = np.radians(declination) + sza_rad = solarposition.solar_zenith_analytical(lat_rad, hour_angle_rad, declination_rad) + + # 批量计算太阳方位角 + saa_rad = solarposition.solar_azimuth_analytical(lat_rad, hour_angle_rad, declination_rad, sza_rad) + + # 转换为度 + sza_deg = np.degrees(sza_rad) + # 约束太阳天顶角到 0~90（太阳在地平线以下时会 >90，这里按需求裁剪） + sza_deg = np.clip(sza_deg, 0.0, 90.0) + + saa_deg = np.degrees(saa_rad) % 360.0 + + # 返回角度数组 + return np.column_stack([sza_deg, saa_deg]) + + + def _extract_swath_name_from_bip_filename(self, bip_file: Path) -> str: + """ + 从输入的 BIP/HDR 文件名中提取“航带名称”（仅包含数字与下划线的前缀）。 + + 示例： + 2025_9_2_3_53_45_202592_35252_0_rad_rgbxyz_geo_registered.bip + -> 2025_9_2_3_53_45_202592_35252_0 + """ + name = bip_file.name + base = name + lower = base.lower() + if lower.endswith('.bip.hdr'): + base = base[:-len('.bip.hdr')] + elif lower.endswith('.hdr'): + base = base[:-len('.hdr')] + elif lower.endswith('.bip'): + base = base[:-len('.bip')] + + match = re.match(r'^(\d+(?:_\d+)*)', base) + if match: + return match.group(1) + + # 兜底：常见格式会在 _rad... 之前出现航带名 + if '_rad' in base: + return base.split('_rad', 1)[0] + return Path(base).stem + + def _find_times_file_for_bip(self, bip_file: Path, times_root_dir: Optional[Path] = None) -> Optional[Path]: + """ + 根据输入的 BIP/HDR 文件名解析“航带名称”，并在 times 根目录下查找对应的 times 文件。 + + 期望目录结构： + times_root_dir/ + 航带名称/ + 航带名称.bil.times + + Args: + bip_file: 输入的 BIP/HDR 文件 + times_root_dir: times 根目录（可选，优先于配置） + + Returns: + times文件路径，如果找不到则返回None + """ + root = Path(times_root_dir) if times_root_dir is not None else self.times_root_dir + + if not root.exists(): + print(f"警告: times根目录不存在 {root}") + return None + + swath_name = self._extract_swath_name_from_bip_filename(bip_file) + expected = root / swath_name / f"{swath_name}.bil.times" + if expected.exists(): + return expected + + # 兼容：旧结构可能直接在根目录下 + flat = root / f"{swath_name}.bil.times" + if flat.exists(): + return flat + + # 再兜底：在子目录中搜索同名 times 文件 + try: + for p in root.glob(f"*/{swath_name}.bil.times"): + return p + for p in root.rglob(f"{swath_name}.bil.times"): + return p + except Exception as e: + print(f"警告: 搜索times文件失败: {e}") + + return None + + def _parse_base_date_from_times_filename(self, times_filename: str) -> datetime: + """ + 从times文件名中解析基准日期 + + Args: + times_filename: times文件名（不含扩展名） + + Returns: + 基准日期时间 + """ + # 示例文件名: 2025_9_2_3_11_52_202592_31037_2 + # 格式: YYYY_M_D_H_M_S_... + + # 去掉扩展名并分割 + parts = times_filename.replace('.bil.times', '').split('_') + + if len(parts) >= 6: + try: + year, month, day, hour, minute, second = map(int, parts[:6]) + return datetime(year, month, day, hour, minute, second) + except (ValueError, IndexError) as e: + print(f"警告: 无法从文件名解析日期 {times_filename}: {e}") + # 返回默认日期 + return datetime(2025, 9, 2, 3, 11, 52) + + print(f"警告: 文件名格式不正确 {times_filename}") + return datetime(2025, 9, 2, 3, 11, 52) + + def process_image(self, bip_file: Path, gps_file: Path, times_root_dir: Optional[Path] = None) -> np.ndarray: + """ + 处理单个BIP图像，计算所有像素的角度 + + 对图像的每一个有效像素（非全零像素）进行角度计算： + - 使用像素在图像中的行列号计算UTM坐标 + - 使用波段4、5、6的信息确定时间和传感器状态 + - 计算太阳天顶角、太阳方位角、传感器天顶角、传感器方位角、相对方位角 + + Args: + bip_file: BIP图像文件路径 + gps_file: GPS数据文件路径 + times_root_dir: times根目录（可选，优先于配置） + + Returns: + 角度结果数组 (lines, samples, 5) - [SZA, SAA, VZA, VAA, RAA] + """ + print(f"\n开始处理图像: {bip_file.name}") + start_time = time.time() + + # 1. 读取图像数据 + image_data, image_info = self.read_bip_image(bip_file) + self.image_info = image_info + self.image_data = image_data + + # 1.5 预处理优化：提前过滤无效像素 + print("预处理：识别有效像素...") + valid_mask = np.any(image_data != 0, axis=2) # 所有波段都不为0的像素 + valid_pixels = np.column_stack(np.where(valid_mask)) # (N, 2) 数组 + print( + f"发现 {len(valid_pixels)} 个有效像素，约占总像素的 {len(valid_pixels) / (image_info['lines'] * image_info['samples']) * 100:.1f}%") + + # 2. 提取像素信息 + col_indices, row_indices, swath_ids = self.extract_pixel_info(image_data, image_info) + + # 3. 加载GPS数据 + if self.gps_data is None: + self.gps_data = self.load_gps_data(gps_file) + + # 4. 获取唯一的航带ID + unique_swath_ids = np.unique(swath_ids[swath_ids > 0]) + print(f"发现航带: {unique_swath_ids}") + + # 5. 从该 BIP/HDR 对应的 times 文件解析基准日期 + base_date = None + if len(unique_swath_ids) > 0: + times_file = self._find_times_file_for_bip(bip_file, times_root_dir=times_root_dir) + if not times_file: + raise ValueError("无法找到该BIP对应的times文件") + base_date = self._parse_base_date_from_times_filename(times_file.stem) + print(f"从times文件解析的基准日期: {base_date}") + else: + raise ValueError("未识别到有效的航带ID，无法关联times文件") + + # 6. 初始化结果数组 + result_shape = (image_info['lines'], image_info['samples'], 5) # SZA, SAA, VZA, VAA, RAA + angle_results = np.full(result_shape, np.nan, dtype=np.float32) + + # 7. 创建航带到times数据的映射 + # 说明：按需求，times 文件由“输入 BIP/HDR 的航带名”唯一确定；若 swath_id 出现多个值，统一使用同一份 times 数据 + times_data_cache = {} + times_file = self._find_times_file_for_bip(bip_file, times_root_dir=times_root_dir) + if not times_file: + raise ValueError("无法找到该BIP对应的times文件") + if str(times_file) not in self.times_data: + self.times_data[str(times_file)] = self.load_times_data(times_file) + shared_times_data = self.times_data[str(times_file)] + for swath_id in unique_swath_ids: + times_data_cache[swath_id] = shared_times_data + + # 8. 对图像的每一个像素进行处理（支持并行处理） + total_pixels = image_info['lines'] * image_info['samples'] + + + # 使用优化的处理（向量化和纯NumPy实现） + print(f"开始计算每个像素的角度... 处理 {len(valid_pixels)} 个有效像素") + print("使用优化的向量化处理模式（纯NumPy实现）") + + angle_results = self._process_pixels_chunked_optimized( + valid_pixels, col_indices, row_indices, swath_ids, + times_data_cache, base_date, image_info + ) + + total_time = time.time() - start_time + processed_pixels = np.sum(~np.isnan(angle_results[:, :, 0])) + print(f"\n图像处理完成，共处理 {processed_pixels} 个有效像素，总耗时: {total_time:.2f} 秒") + print(f"平均处理速度: {processed_pixels / total_time:.0f} 像素/秒") + return angle_results + + def save_bil_results(self, angle_results: np.ndarray, output_file: Path, image_info: Dict, original_bip_file: Path = None): + """ + 将角度结果保存为ENVI BSQ格式（使用GDAL），HDR文件基于原始图像 + + Args: + angle_results: 角度结果数组 (lines, samples, 3) + output_file: 输出文件路径 + image_info: 原始图像信息 + original_bip_file: 原始BIP文件路径，用于复制HDR文件 + """ + print(f"使用GDAL保存结果为ENVI BSQ格式: {output_file}") + + from osgeo import gdal, osr + + # 获取图像尺寸 + lines, samples, bands = angle_results.shape + assert bands == 5, "角度结果必须有5个波段" + + # 创建GDAL驱动 + driver = gdal.GetDriverByName('ENVI') + + # 创建输出数据集 + dataset = driver.Create(str(output_file), samples, lines, bands, gdal.GDT_Float32, + options=['INTERLEAVE=BSQ']) + + if dataset is None: + raise RuntimeError(f"无法创建输出文件: {output_file}") + + # 设置地理变换参数（如果有的话） + if 'geo_transform' in image_info: + dataset.SetGeoTransform(image_info['geo_transform']) + + # 设置投影信息 + if 'projection' in image_info: + srs = osr.SpatialReference() + # 如果是WKT格式 + if isinstance(image_info['projection'], str) and 'PROJCS' in image_info['projection']: + srs.ImportFromWkt(image_info['projection']) + else: + # 尝试设置为UTM + utm_zone = image_info.get('utm_zone', 51) + if isinstance(utm_zone, str): + utm_zone = int(''.join(filter(str.isdigit, utm_zone))) + hemisphere = image_info.get('hemisphere', 'North') + + if hemisphere.upper().startswith('N'): + epsg_code = 32600 + utm_zone + else: + epsg_code = 32700 + utm_zone + + srs.ImportFromEPSG(epsg_code) + + dataset.SetProjection(srs.ExportToWkt()) + + # 设置波段数据和元数据 + band_names = ['Solar Zenith Angle (SZA)', 'Solar Azimuth Angle (SAA)', 'Sensor Zenith Angle (VZA)', 'Sensor Azimuth Angle (VAA)', 'Relative Azimuth Angle (RAA)'] + + for band_idx in range(bands): + band = dataset.GetRasterBand(band_idx + 1) + band_data = angle_results[:, :, band_idx] + + # 处理NaN值 + band_data = np.nan_to_num(band_data, nan=0) + + # 写入数据 + band.WriteArray(band_data) + band.SetNoDataValue(0) + + # 设置波段描述 + band.SetDescription(band_names[band_idx]) + + # 设置数据集级别的元数据 + dataset.SetMetadataItem('description', 'Angle calculation results for UAV hyperspectral imagery') + dataset.SetMetadataItem('band_names', '{' + ', '.join(band_names) + '}') + + # 清理资源 + dataset.FlushCache() + dataset = None + + # GDAL会自动生成HDR文件，现在我们用基于原始HDR的自定义HDR文件覆盖它 + self._create_hdr_from_original(output_file, original_bip_file, bands, band_names) + + print("结果保存完成") + + def save_bip_with_angles(self, angle_results: np.ndarray, output_file: Path, + image_info: Dict, original_bip_file: Path): + """ + 将角度结果追加到原始BIP数据后面，保存为新的BIP文件 + """ + print(f"保存带角度的BIP文件: {output_file}") + + from osgeo import gdal, osr + + if self.image_data is None or self.image_info is None: + self.image_data, self.image_info = self.read_bip_image(original_bip_file) + + lines, samples, orig_bands = self.image_data.shape + angle_bands = angle_results.shape[2] + total_bands = orig_bands + angle_bands + + driver = gdal.GetDriverByName('ENVI') + dataset = driver.Create(str(output_file), samples, lines, total_bands, gdal.GDT_Float32, + options=['INTERLEAVE=BIP']) + if dataset is None: + raise RuntimeError(f"无法创建输出文件: {output_file}") + + if 'geo_transform' in image_info: + dataset.SetGeoTransform(image_info['geo_transform']) + + if 'projection' in image_info: + srs = osr.SpatialReference() + if isinstance(image_info['projection'], str) and 'PROJCS' in image_info['projection']: + srs.ImportFromWkt(image_info['projection']) + else: + utm_zone = image_info.get('utm_zone', 51) + if isinstance(utm_zone, str): + utm_zone = int(''.join(filter(str.isdigit, utm_zone))) + hemisphere = image_info.get('hemisphere', 'North') + epsg_code = 32600 + utm_zone if hemisphere.upper().startswith('N') else 32700 + utm_zone + srs.ImportFromEPSG(epsg_code) + dataset.SetProjection(srs.ExportToWkt()) + + # 写入原始波段 + original_data = self.image_data.astype(np.float32, copy=False) + for band_idx in range(orig_bands): + band = dataset.GetRasterBand(band_idx + 1) + band.WriteArray(original_data[:, :, band_idx]) + + # 写入角度波段 + angle_band_names = [ + 'Solar Zenith Angle (SZA)', + 'Solar Azimuth Angle (SAA)', + 'Sensor Zenith Angle (VZA)', + 'Sensor Azimuth Angle (VAA)', + 'Relative Azimuth Angle (RAA)' + ] + for idx in range(angle_bands): + band = dataset.GetRasterBand(orig_bands + idx + 1) + band_data = np.nan_to_num(angle_results[:, :, idx], nan=0) + band.WriteArray(band_data) + band.SetNoDataValue(0) + if idx < len(angle_band_names): + band.SetDescription(angle_band_names[idx]) + + dataset.FlushCache() + dataset = None + + # 生成HDR，保留原始band names并追加角度band names + if original_bip_file is not None: + hdr_file = original_bip_file if original_bip_file.suffix.lower() == '.hdr' else original_bip_file.with_suffix('.hdr') + orig_names = self._extract_band_names_from_hdr(hdr_file, orig_bands) + else: + orig_names = [f'Band {i + 1}' for i in range(orig_bands)] + + appended_names = orig_names + angle_band_names[:angle_bands] + self._create_hdr_from_original(output_file, original_bip_file, total_bands, appended_names) + + def _create_hdr_from_original(self, output_file: Path, original_bip_file: Path, bands: int, band_names: list): + """ + 基于原始HDR文件创建新的HDR文件 + + Args: + output_file: 输出文件路径 + original_bip_file: 原始BIP文件路径 + bands: 输出文件的波段数 + band_names: 波段名称列表 + """ + # 查找原始HDR文件 + if original_bip_file is None: + print(f"警告: 未提供原始文件路径，将创建基本的HDR文件") + return + + # 如果输入已经是hdr文件 + if original_bip_file.suffix.lower() == '.hdr': + original_hdr_file = original_bip_file + else: + # 尝试不同的HDR文件命名方式 + original_hdr_file = original_bip_file.with_suffix('.hdr') + if not original_hdr_file.exists(): + original_hdr_file = original_bip_file.with_suffix('.bip.hdr') + if not original_hdr_file.exists(): + original_hdr_file = Path(str(original_bip_file) + '.hdr') + + if not original_hdr_file.exists(): + print(f"警告: 找不到原始HDR文件 {original_hdr_file}，将创建基本的HDR文件") + self._create_basic_hdr_file(output_file, bands, band_names) + return + + output_hdr_file = output_file.with_suffix('.hdr') + + print(f"基于原始HDR文件 {original_hdr_file} 创建新的HDR文件") + + try: + with open(original_hdr_file, 'r', encoding='utf-8', errors='ignore') as f: + hdr_content = f.read() + + # 分割为行 + lines = hdr_content.split('\n') + + # 修改相关参数 + modified_lines = [] + skip_wavelength = False + + for line in lines: + line_lower = line.lower().strip() + + # 修改波段数 + if line_lower.startswith('bands'): + modified_lines.append(f'bands = {bands}') + continue + + # 修改数据类型为Float32 (4) + if line_lower.startswith('data type'): + modified_lines.append('data type = 4') + continue + + # 添加或修改data ignore value + if line_lower.startswith('data ignore value'): + modified_lines.append('data ignore value = 0') + continue + + # 修改交织方式为BSQ + if line_lower.startswith('interleave'): + modified_lines.append('interleave = bip') + continue + + # 跳过wavelength相关行 + if line_lower.startswith('wavelength') or line_lower.startswith('wavelength units'): + skip_wavelength = True + continue + elif skip_wavelength and (line.strip() == '' or line.strip().startswith('}')): + skip_wavelength = False + continue + elif skip_wavelength: + continue + + # 跳过band names（后面会重新添加） + if line_lower.startswith('band names'): + continue + + modified_lines.append(line) + + # 添加新的band names + modified_lines.insert(-1, f'band names = {{') + for i, name in enumerate(band_names): + if i < len(band_names) - 1: + modified_lines.insert(-1, f' {name},') + else: + modified_lines.insert(-1, f' {name}') + modified_lines.insert(-1, '}') + + # 更新description + for i, line in enumerate(modified_lines): + if line.lower().strip().startswith('description'): + modified_lines[i] = 'description = {Angle calculation results for UAV hyperspectral imagery}' + break + + # 写回HDR文件 + with open(output_hdr_file, 'w', encoding='utf-8') as f: + f.write('\n'.join(modified_lines)) + + print(f"HDR文件创建完成: {output_hdr_file}") + + except Exception as e: + print(f"复制HDR文件失败: {e}，将创建基本的HDR文件") + self._create_basic_hdr_file(output_file, bands, band_names) + + def _create_basic_hdr_file(self, output_file: Path, bands: int, band_names: list): + """ + 创建基本的HDR文件 + + Args: + output_file: 输出文件路径 + bands: 波段数 + band_names: 波段名称列表 + """ + output_hdr_file = output_file.with_suffix('.hdr') + + # 创建基本的HDR文件 + basic_hdr = f"""ENVI +description = {{ +Angle calculation results for UAV hyperspectral imagery +}} +samples = {self.image_info['samples']} +lines = {self.image_info['lines']} +bands = {bands} +header offset = 0 +file type = ENVI Standard +data type = 4 +interleave = bip +byte order = 0 +band names = {{ +{chr(10).join(f' {name}' for name in band_names)} +}} +""" + + # 添加地图信息（如果有） + if 'geo_transform' in self.image_info: + geo_transform = self.image_info['geo_transform'] + map_info_parts = [ + self.image_info.get('projection', 'UTM'), + '1.00000', # dummy x start + '1.00000', # dummy y start + str(geo_transform[0]), # x_ul + str(geo_transform[3]), # y_ul + str(geo_transform[1]), # pixel_x + str(abs(geo_transform[5])), # pixel_y + str(self.image_info.get('utm_zone', 51)), + self.image_info.get('hemisphere', 'North'), + 'WGS-84', + 'units=m', + 'rotation=0.000' + ] + basic_hdr += f"map info = {{{', '.join(map_info_parts)}}}\n" + + with open(output_hdr_file, 'w', encoding='utf-8') as f: + f.write(basic_hdr) + + print(f"基本HDR文件创建完成: {output_hdr_file}") + + +def main(): + """主函数""" + # 配置参数 + config = { + 'data_dir': r'D:\BaiduNetdiskDownload\20250902\_3_52_52\Geoout\Geoout', + # times根目录：内部包含多个“航带名”子文件夹，每个子文件夹内包含 航带名.bil.times + 'times_root_dir': r"D:\BaiduNetdiskDownload\20250902\_3_52_52\202592_35252", + # GPS文件：所有HDR文件共用这一个GPS文件 + 'gps_file': r"D:\BaiduNetdiskDownload\20250902\_3_52_52\2025_9_2_3_53_45.gps", + 'output_dir': r'D:\BaiduNetdiskDownload\20250902\_3_52_52\316\agnle', + 'base_height': 254, # 基准高度 + } + + # 创建计算器 + calculator = UAVAngleCalculator(config) + + data_dir = Path(config['data_dir']) + times_root_dir = Path(config['times_root_dir']) + gps_file = Path(config['gps_file']) + + if not gps_file.exists(): + raise FileNotFoundError(f"GPS文件不存在: {gps_file}") + + # 批量处理：自动查找 data_dir 下所有 .hdr 文件 + hdr_files = sorted(data_dir.glob("*.hdr")) + if not hdr_files: + raise FileNotFoundError(f"在目录中未找到任何 .hdr 文件: {data_dir}") + + # 使用进度条处理文件 + if HAS_TQDM: + hdr_iterator = tqdm(hdr_files, desc="处理HDR文件", unit="文件") + else: + hdr_iterator = hdr_files + print(f"开始处理 {len(hdr_files)} 个HDR文件...") + + for hdr_file in hdr_iterator: + try: + if HAS_TQDM: + hdr_iterator.set_description(f"处理: {hdr_file.name}") + else: + print(f"\n处理文件: {hdr_file.name}") + + angle_results = calculator.process_image(hdr_file, gps_file, times_root_dir=times_root_dir) + + base_name = hdr_file.name + if base_name.lower().endswith('.bip.hdr'): + base_name = base_name[:-len('.bip.hdr')] + elif base_name.lower().endswith('.hdr'): + base_name = base_name[:-len('.hdr')] + + output_file = Path(config['output_dir']) / f"{base_name}_angles.bip" + calculator.save_bip_with_angles(angle_results, output_file, calculator.image_info, hdr_file) + + if not HAS_TQDM: + print(f"完成: {hdr_file.name} -> {output_file.name}") + except Exception as e: + if HAS_TQDM: + hdr_iterator.set_description(f"失败: {hdr_file.name}") + print(f"处理失败: {hdr_file.name}: {e}") + import traceback + traceback.print_exc() + continue + + print("\n处理完成!") + print(f"输出目录: {config['output_dir']}") + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/ocbrdf/ocbrdf_main3.py b/ocbrdf/ocbrdf_main3.py index 6407ea1..8836f55 100644 --- a/ocbrdf/ocbrdf_main3.py +++ b/ocbrdf/ocbrdf_main3.py @@ -177,11 +177,11 @@ def _block_windows(src, block_size=512): def _resolve_output_envi_path(output_file): """Return the ENVI image path used for rasterio writes.""" - if output_file.lower().endswith(".img"): + if output_file.lower().endswith(".dat"): return output_file if output_file.lower().endswith(".hdr"): - return output_file[:-4] + ".img" - return output_file + ".img" + return output_file[:-4] + ".dat" + return output_file + ".dat" def _read_wavelengths_from_header(hyperspectral_file): @@ -259,7 +259,7 @@ def _get_output_filename(base_output_file: str, var_name: str) -> str: input.hdr + brdf_unc -> input_brdf_unc.hdr """ base, ext = os.path.splitext(base_output_file) - if ext.lower() in ('.hdr', '.img'): + if ext.lower() in ('.hdr', '.dat'): base = base # keep base without extension suffix = f"_{var_name}" if var_name.lower() != "rw_brdf" else "" return f"{base}{suffix}{ext}" @@ -486,6 +486,29 @@ def run_brdf_correction_block( for dst in output_dsts.values(): dst.close() + # Write ENVI header files with wavelength information from input + for var_name, out_path in output_paths.items(): + hdr_path = out_path[:-4] + ".hdr" if out_path.lower().endswith(".dat") else out_path + ".hdr" + # Prepare ENVI header content with wavelengths from input file + header_lines = [ + "ENVI", + f"description = {{BRDF corrected result: {var_name}}}", + f"samples = {n_cols}", + f"lines = {n_rows}", + f"bands = {n_bands}", + "header offset = 0", + "file type = ENVI Standard", + "data type = 4", # float32 + "interleave = bsq", + "byte order = 0", # little-endian + f"wavelength = {{{', '.join(str(v) for v in wavelengths.tolist())}}}", + "wavelength units = nm", + ] + # Write header file + with open(hdr_path, 'w') as f: + f.write('\n'.join(header_lines)) + print(f" Written header: {hdr_path}") + # Print performance statistics total_blocks = blocks_processed + blocks_skipped if total_blocks > 0: @@ -864,7 +887,7 @@ def main(): parser.add_argument("hyperspectral_file", help="Input ENVI BSQ hyperspectral file or .hdr path.") parser.add_argument("angle_file", help="Input ENVI BIP angle file or .hdr path.") parser.add_argument("mask_file", help="Input water mask GeoTIFF; pixels == 1 are corrected.") - parser.add_argument("output_file", help="Output path prefix. ENVI writes .hdr/.img, NetCDF writes .nc.") + parser.add_argument("output_file", help="Output path prefix. ENVI writes .hdr/.dat, NetCDF writes .nc.") parser.add_argument("--brdf-model", default="L11", choices=["L11", "M02", "M02SeaDAS", "O25"]) parser.add_argument("--output-var", nargs='+', default=["Rw_brdf"], help="Output variables to save (one or more). " diff --git a/requirements.txt b/requirements.txt index 7fcdae0..716de42 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,14 +1,42 @@ +# BRDF_GUI 项目依赖 +# 安装命令: pip install -r requirements.txt +# +# 本项目包含两个主要模块: +# 1. 陆地植被BRDF校正 (FlexBRDF) - 使用 HyTools 库 +# 2. 水体BRDF校正 (Ocean BRDF) - 使用 ocbrdf 模块 + +# ========== 核心科学计算 ========== numpy>=1.21.0 -ray>=2.0.0 -h5py>=3.0.0 -h5netcdf>=0.13.0 -matplotlib>=3.5.0 -pandas>=1.3.0 -scikit-learn>=1.0.0 scipy>=1.7.0 +pandas>=1.3.0 xarray>=0.20.0 -rasterio>=1.2.0 -spectral>=0.22.0 + +# ========== GUI 界面 ========== +PyQt5>=5.15.0 + +# ========== 高光谱数据处理 ========== +spectral>=0.22.0 # ENVI文件读写 +rasterio>=1.3.0 # 栅格数据处理 (GeoTIFF/ENVI) +geopandas>=0.10.0 # 地理空间数据处理 +shapely>=2.0.0 # 几何形状处理 + +# ========== HDF5/NetCDF 支持 ========== +h5py>=3.0.0 # HDF5文件支持 +h5netcdf>=0.13.0 # NetCDF via HDF5 +netCDF4>=1.6.0 # NetCDF文件支持 (可选) + +# ========== 并行处理 ========== +ray>=2.0.0 + +# ========== 机器学习 ========== +scikit-learn>=1.0.0 + +# ========== 可视化 ========== +matplotlib>=3.5.0 + +# ========== 观测几何计算 ========== +pvlib>=0.10.0 # 太阳位置计算 +pyproj>=3.4.0 # 坐标投影转换 + +# ========== 进度条/工具 ========== tqdm>=4.60.0 -geopandas>=0.10.0 -shapely>=2.0.0 diff --git a/slope/caijian.py b/slope/caijian.py new file mode 100644 index 0000000..e1a7ef6 --- /dev/null +++ b/slope/caijian.py @@ -0,0 +1,103 @@ +import rasterio +from rasterio.warp import reproject, Resampling +import numpy as np + +def merge_angles_with_slope_aspect(angle_file, slope_aspect_file, output_file, + resampling=Resampling.bilinear): + """ + 将坡度坡向波段添加到角度文件波段之后，输出合并的ENVI DAT文件。 + 输出文件的空间范围、分辨率、投影与角度文件完全一致。 + 掩膜基于角度文件的忽略值生成：只有所有角度波段均为忽略值(0.0)的像素才被设为忽略值。 + 输出文件的忽略值统一设为0.0。 + + Parameters + ---------- + angle_file : str + 包含角度信息的ENVI DAT文件（如传感器角度、太阳角度等，至少1个波段）。 + 该文件定义了输出的目标空间参数和掩膜基准。 + slope_aspect_file : str + 包含坡度坡向的ENVI DAT文件（通常2个波段）。 + output_file : str + 输出合并后的ENVI DAT文件路径（将自动生成.hdr头文件）。 + resampling : rasterio.warp.Resampling, optional + 重采样方法，默认为双线性（适用于连续参数）。 + """ + # 1. 从角度文件获取目标空间参数、数据、掩膜 + with rasterio.open(angle_file) as src_angle: + # 目标空间参数 + dst_transform = src_angle.transform + dst_width = src_angle.width + dst_height = src_angle.height + dst_crs = src_angle.crs + + # 读取角度文件所有波段 + angle_data = src_angle.read() # shape: (bands, rows, cols) + angle_count = src_angle.count + + # 获取角度文件的忽略值（如果未定义，默认为0.0，但应给出警告） + angle_nodata = src_angle.nodata + if angle_nodata is None: + print("警告：角度文件未定义忽略值，将使用0.0作为默认忽略值。") + angle_nodata = 0.0 + + # 构建掩膜：所有波段均不等于忽略值视为有效 + mask = np.all(angle_data != angle_nodata, axis=0) + + # 2. 处理坡度坡向文件：重采样到目标空间 + with rasterio.open(slope_aspect_file) as src_sa: + src_sa_data = src_sa.read() # shape: (bands, src_sa_height, src_sa_width) + src_sa_count = src_sa.count + src_sa_crs = src_sa.crs + src_sa_transform = src_sa.transform + + sa_resampled = np.empty((src_sa_count, dst_height, dst_width), + dtype=src_sa_data.dtype) + + reproject( + source=src_sa_data, + destination=sa_resampled, + src_transform=src_sa_transform, + src_crs=src_sa_crs, + dst_transform=dst_transform, + dst_crs=dst_crs, + resampling=resampling + ) + + # 3. 合并波段：角度在前，坡度坡向在后 + merged_data = np.concatenate([angle_data, sa_resampled], axis=0) + merged_count = merged_data.shape[0] + + # 4. 准备输出元数据（基于角度文件的profile） + with rasterio.open(angle_file) as src_angle: + out_profile = src_angle.profile.copy() + output_nodata = 0.0 # 强制输出忽略值为0.0（用户要求） + out_profile.update({ + 'driver': 'ENVI', + 'height': dst_height, + 'width': dst_width, + 'transform': dst_transform, + 'crs': dst_crs, + 'count': merged_count, + 'dtype': merged_data.dtype, + 'nodata': output_nodata + }) + + # 5. 应用掩膜：将无效区域的像素设为忽略值0.0 + for b in range(merged_count): + merged_data[b, ~mask] = output_nodata + + # 6. 写入输出文件 + with rasterio.open(output_file, 'w', **out_profile) as dst: + dst.write(merged_data) + + print(f"处理完成！合并后的文件已保存至：{output_file}") + print(f"总波段数：{merged_count}（角度 {angle_count} 波段 + 坡度坡向 {src_sa_count} 波段）") + print(f"有效像素区域基于角度文件所有波段判断，忽略值统一设为 0.0。") + +# 示例用法 +if __name__ == "__main__": + angle_file = r"E:\is2\yaopu\output\angel.dat" # 角度DAT文件 + slope_aspect_file = r"E:\is2\yaopu\DEM\slope_aspect_cosine.dat" # 原始坡度坡向DAT + out_file = r"E:\is2\yaopu\DEM\angles_with_slope_aspect.dat" # 合并输出 + + merge_angles_with_slope_aspect(angle_file, slope_aspect_file, out_file) \ No newline at end of file diff --git a/slope/readme.txt b/slope/readme.txt new file mode 100644 index 0000000..69d0e42 --- /dev/null +++ b/slope/readme.txt @@ -0,0 +1,197 @@ +================================================================================ + 地形校正模块 - 已知问题与改进方向 +================================================================================ + +【问题描述】 +----------- +当前版本的 slope_aspectV2.py 在处理推扫式无人机高光谱图像时，存在一个关键假设： +整幅图像使用单一的太阳天顶角和方位角来计算 cosine_i（入射角余弦）。 + +然而，实际数据采集情况是： +- 高光谱图像由3个架次的无人机飞行拼接而成 +- 每个架次的飞行时间不同（通常间隔数分钟至数十分钟） +- 因此，不同扫描线（或不同图像区域）对应的太阳角度是不同的 +- 但目前代码仅基于DEM中心点坐标计算一组太阳角度，并应用于整幅图像 + +【技术影响】 +----------- +这个简化会导致以下问题： + +1. 时间跨度误差 + - 假设3个架次间隔30分钟，太阳方位角可能变化7-8度（取决于纬度） + - 在日出/日落时段，太阳高度角变化可达10-15度 + - 这会导致不同架次交界处的 cosine_i 出现系统性偏差 + +2. 地形校正不准确 + - 东向坡面：过早架次会低估照明，过晚架次会高估照明 + - 西向坡面：与东向坡面相反 + - 交界区域会出现明显的"拼接缝"现象 + +3. BRDF校正累积误差 + - cosine_i 是许多BRDF模型的关键输入参数 + - 角度误差会传播到后续的BRDF地形校正结果 + +【当前代码实现】 +--------------- +文件: slope/slope_aspectV2.py + +当前流程： + 1. 读取DEM -> 2. 计算坡度/坡向 -> 3. 获取DEM中心点坐标(lat, lon) + 4. 根据单一时间计算一组太阳角度(solar_zn, solar_az) + 5. 将这组角度应用于整幅图像的所有像素 + +关键代码段（约110-130行）： + ```python + # 获取DEM中心点坐标 + latitude, longitude = get_dem_center_coords(dem) # 仅一个中心点！ + + # 计算太阳角度（仅一组） + solar_zn_deg, solar_az_deg = calc_solar_angles(latitude, longitude, observation_time) + + # 应用于整幅图像 + cosine_i = calc_cosine_i(solar_zn_rad, solar_az_rad, aspect_rad, slope_rad) + ``` + +【需要的改进方向】 +----------------- +接手人员需要进行以下修改： + +方案A: 分架次处理（推荐） +---------------------------- +如果已知3个架次的地理分界线或时间范围： + + 1. 输入格式修改 + - 接收3组时间参数（或从飞行日志自动提取） + - 例如: -t1 "2024-06-15 08:30:00" -t2 "2024-06-15 09:00:00" -t3 "2024-06-15 09:30:00" + + 2. 空间分区处理 + - 根据每架次的覆盖范围（可从高光谱图像的GLT或地理坐标确定） + - 将DEM分割为3个区域 + - 每个区域使用对应架次的时间计算太阳角度 + + 3. 修改后的伪代码示例： + ```python + # 定义3个架次的时间和地理范围 + flight_segments = [ + {'time': '2024-06-15 08:30:00', 'row_range': (0, 1200)}, # 架次1: 行0-1199 + {'time': '2024-06-15 09:00:00', 'row_range': (1200, 2800)}, # 架次2: 行1200-2799 + {'time': '2024-06-15 09:30:00', 'row_range': (2800, 4000)}, # 架次3: 行2800-3999 + ] + + # 初始化cosine_i数组 + cosine_i = np.zeros_like(slope_data) + + # 分段计算 + for segment in flight_segments: + dt = parse_datetime(segment['time']) + solar_zn, solar_az = calc_solar_angles(lat, lon, dt) + + y1, y2 = segment['row_range'] + cosine_i[y1:y2, :] = calc_cosine_i( + np.deg2rad(solar_zn), np.deg2rad(solar_az), + aspect_rad[y1:y2, :], slope_rad[y1:y2, :] + ) + ``` + +方案B: 逐行/逐像素太阳角度（更精确） +------------------------------------- +如果需要更高精度（考虑太阳角度在单架次内的缓慢变化）： + + 1. 输入要求 + - 高光谱图像的每一行对应的时间戳（可从IMU/GPS数据提取） + - 或者：每个像素的精确UTC时间 + + 2. 逐行计算 + ```python + # 假设times是一个与图像行数相同长度的一维数组 + # times[i] = 第i行的观测时间（datetime对象） + + for row in range(n_rows): + dt = observation_times[row] # 该行的具体时间 + solar_zn, solar_az = calc_solar_angles(lat, lon, dt) + + cosine_i[row, :] = calc_cosine_i( + np.deg2rad(solar_zn), np.deg2rad(solar_az), + aspect_rad[row, :], slope_rad[row, :] + ) + ``` + + 3. 性能优化 + - 由于pvlib计算较慢，可对时间相近的行批量计算 + - 例如：每10行计算一次太阳角度，中间行插值 + +方案C: 从飞行日志自动提取（最自动化） +-------------------------------------- +如果飞行平台记录了GPS/IMU数据： + + 1. 输入数据 + - 读取POS数据（位置、姿态、时间） + - 或从高光谱图像的GLT（地理位置查找表）提取 + + 2. 建立时间-位置映射 + - 根据无人机位置和时间建立每行图像的观测时间 + + 3. 分区域计算cosine_i + +【数据结构参考】 +--------------- +典型的推扫式高光谱数据结构： + +高光谱图像 (lines, samples, bands) + - lines: 扫描线数（飞行方向的像素数） + - samples: 每行的像素数（垂直飞行方向） + - 第i行对应无人机的第i个位置 + +关键信息来源： + 1. 图像头文件 (.hdr) - 通常包含起始时间和积分时间 + 2. GLT文件 (Geographic Lookup Table) - 每个像素的地理坐标 + 3. 飞行日志/POS数据 - 每行对应的时间戳 + 4. 帧辅助数据 (Frame Auxiliary Data) - 部分传感器自带每帧时间 + +【建议的输入接口修改】 +---------------------- +当前命令行接口： + python slope_aspectV2.py -i dem.tif -o output.dat -t "2024-06-15 10:00:00" + +建议扩展为： + # 方式1: 手动指定3个架次的时间和行范围 + python slope_aspectV2.py -i dem.tif -o output.dat \ + --segments "0,1200,2024-06-15T08:30:00" \ + --segments "1200,2800,2024-06-15T09:00:00" \ + --segments "2800,4000,2024-06-15T09:30:00" + + # 方式2: 从飞行日志文件自动读取 + python slope_aspectV2.py -i dem.tif -o output.dat --time-log flight_log.csv + + # 方式3: 从高光谱图像头文件提取时间信息 + python slope_aspectV2.py -i dem.tif -o output.dat --hyperspectral-ref image.hdr + +【验证方法】 +----------- +修改后应验证以下指标： + +1. 连续性检查 + - 目视检查3个架次交界处的cosine_i是否连续 + - 统计交界行前后10行的cosine_i差异 + +2. 物理合理性 + - 东向坡的cosine_i应随时间增加而增加（上午飞行） + - 西向坡的cosine_i应随时间增加而减小（上午飞行） + +3. 与GLT对比 + - 如果高光谱有GLT，可根据每像素的太阳角度公式验证 + +【相关文件】 +----------- +- slope_aspectV2.py - 当前主程序（需要修改） +- angle_compute.py - 角度计算工具（可能包含相关逻辑） +- 高光谱图像的.hdr文件 - 可能包含时间元数据 +- 飞行日志/POS数据 - 外部时间参考 + +【优先级建议】 +-------------- +1. 【高】首先实现方案A（分架次处理）- 工作量适中，效果显著 +2. 【中】评估是否需要方案B（逐行）- 取决于太阳角度变化幅度 +3. 【低】长期可考虑方案C（全自动）- 需要整合飞行数据解析 + +简单来说在计算的时候，需要加入角度文件，进行逐像素的cosi计算 diff --git a/slope/slope_aspect.py b/slope/slope_aspect.py new file mode 100644 index 0000000..d96cfce --- /dev/null +++ b/slope/slope_aspect.py @@ -0,0 +1,83 @@ +import xdem +import rasterio +import numpy as np + +def calc_cosine_i(solar_zn, solar_az, aspect, slope): + """Generate cosine i image. The cosine of the incidence angle (i) is + defined as the angle between the normal to the pixel surface + and the solar zenith direction. + All input geometry units must be in radians. + + Args: + solar_az (numpy.ndarray): Solar azimuth angle. + solar_zn (numpy.ndarray): Solar zenith angle. + aspect (numpy.ndarray): Ground aspect. + slope (numpy.ndarray): Ground slope. + + Returns: + numpy.ndarray: Cosine i image. + """ + relative_az = aspect - solar_az + cosine_i = (np.cos(solar_zn) * np.cos(slope) + + np.sin(solar_zn) * np.sin(slope) * np.cos(relative_az)) + return cosine_i + +# ========== 1. 读取 DEM ========== +dem = xdem.DEM("E:/is2/yaopu/DEM/dsm.tif") + +# ========== 2. 计算坡度、坡向 ========== +slope = xdem.terrain.slope(dem, method='ZevenbergThorne') +aspect = xdem.terrain.aspect(dem) + +# ========== 3. 获取空间参考信息 ========== +transform = dem.transform +crs = dem.crs + +# ========== 4. 转换为 numpy 数组 ========== +slope_data = slope.data +aspect_data = aspect.data + +# ========== 5. 太阳角度（单位：度，根据实际飞行时间和地理位置设置） ========== +solar_zn_deg = 30.0 # 太阳天顶角（从天顶到太阳的夹角），例如 30° +solar_az_deg = 180.0 # 太阳方位角（从北顺时针），例如 180°（正南） + +# 转换为弧度（因为三角函数需要弧度） +solar_zn_rad = np.deg2rad(solar_zn_deg) +solar_az_rad = np.deg2rad(solar_az_deg) + +# 坡度和坡向也转换为弧度 +slope_rad = np.deg2rad(slope_data) +aspect_rad = np.deg2rad(aspect_data) + +# ========== 6. 计算 cosine_i ========== +cosine_i = calc_cosine_i(solar_zn_rad, solar_az_rad, aspect_rad, slope_rad) + +# ========== 7. 处理无效值（NaN） ========== +# 注意：将 NaN 替换为 0 可能会与真实值为 0 的像元混淆（如水平面坡度为0，阴影边缘cosine_i=0）。 +# 如需严格区分，建议使用特殊值（如 -9999）作为 NoData。 +slope_data = np.nan_to_num(slope_data, nan=0.0) +aspect_data = np.nan_to_num(aspect_data, nan=0.0) +cosine_i = np.nan_to_num(cosine_i, nan=0.0) + +# ========== 8. 堆叠为多波段数组（坡度、坡向、cosine_i） ========== +stacked = np.stack([slope_data, aspect_data, cosine_i], axis=0) + +# ========== 9. 输出 ENVI 格式文件 ========== +output_path = "E:/is2/yaopu/DEM/slope_aspect_cosine.dat" + +with rasterio.open( + output_path, + 'w', + driver='ENVI', + height=stacked.shape[1], + width=stacked.shape[2], + count=stacked.shape[0], + dtype=stacked.dtype, + crs=crs, + transform=transform, + nodata=0.0, # 与上面替换的 NaN 值一致 +) as dst: + dst.write(stacked) + +print(f"已保存多波段 ENVI 文件：{output_path}") +print(f"波段顺序：1-坡度(°), 2-坡向(°), 3-cosine_i") \ No newline at end of file diff --git a/slope/slope_aspectV2.py b/slope/slope_aspectV2.py new file mode 100644 index 0000000..2fad3ba --- /dev/null +++ b/slope/slope_aspectV2.py @@ -0,0 +1,381 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +DEM地形处理工具 - 计算坡度、坡向和太阳入射角余弦 + +该工具从DEM文件中提取坐标信息，根据给定时间计算太阳位置， +并生成坡度、坡向和cosine_i的多波段ENVI文件。 + +作者: BRDF_GUI Team +版本: 2.0.0 +""" + +import xdem +import rasterio +import numpy as np +import pandas as pd +import argparse +import sys +from datetime import datetime +from pathlib import Path +import pvlib.solarposition as solarposition +from rasterio.crs import CRS +import pyproj + + +def calc_cosine_i(solar_zn, solar_az, aspect, slope): + """ + 计算入射角余弦（cosine i）。 + + 入射角余弦定义为地表法线与太阳天顶方向之间夹角的余弦值。 + 所有角度输入必须为弧度。 + + Args: + solar_az (numpy.ndarray): 太阳方位角（弧度） + solar_zn (numpy.ndarray): 太阳天顶角（弧度） + aspect (numpy.ndarray): 地面坡向（弧度） + slope (numpy.ndarray): 地面坡度（弧度） + + Returns: + numpy.ndarray: 入射角余弦图像 + """ + relative_az = aspect - solar_az + cosine_i = (np.cos(solar_zn) * np.cos(slope) + + np.sin(solar_zn) * np.sin(slope) * np.cos(relative_az)) + return cosine_i + + +def get_dem_center_coords(dem): + """ + 获取DEM图像中心点的经纬度坐标。 + + Args: + dem: xdem.DEM对象 + + Returns: + tuple: (纬度, 经度) - 单位：度 + """ + # 获取DEM的范围边界 + bounds = dem.bounds # (left, bottom, right, top) + center_x = (bounds.left + bounds.right) / 2 + center_y = (bounds.bottom + bounds.top) / 2 + + # 如果DEM是投影坐标系，需要转换为WGS84经纬度 + if dem.crs is not None and not dem.crs.is_geographic: + # 创建从投影坐标系到WGS84的转换 + transformer = pyproj.Transformer.from_crs(dem.crs, CRS.from_epsg(4326), always_xy=True) + lon, lat = transformer.transform(center_x, center_y) + else: + # 已经是地理坐标系（经纬度） + lon, lat = center_x, center_y + + return lat, lon + + +def parse_datetime(datetime_str): + """ + 解析日期时间字符串。 + + 支持的格式： + - YYYY-MM-DD HH:MM:SS + - YYYY-MM-DDTHH:MM:SS + - YYYYMMDD_HHMMSS + + Args: + datetime_str (str): 日期时间字符串 + + Returns: + datetime: datetime对象（UTC时间） + """ + formats = [ + "%Y-%m-%d %H:%M:%S", + "%Y-%m-%dT%H:%M:%S", + "%Y%m%d_%H%M%S", + "%Y-%m-%d", + ] + + for fmt in formats: + try: + dt = datetime.strptime(datetime_str, fmt) + # 如果只输入了日期，默认时间为正午12:00 + if fmt == "%Y-%m-%d": + dt = dt.replace(hour=12, minute=0, second=0) + return dt + except ValueError: + continue + + raise ValueError(f"无法解析日期时间: {datetime_str}. 支持的格式: YYYY-MM-DD HH:MM:SS, YYYY-MM-DDTHH:MM:SS, YYYYMMDD_HHMMSS") + + +def get_user_input_datetime(): + """ + 交互式获取用户输入的日期和时间。 + + Returns: + datetime: datetime对象（UTC时间） + """ + print("请输入观测日期和时间（UTC时间）：") + + try: + year = int(input(" 年份 (YYYY): ")) + month = int(input(" 月份 (MM): ")) + day = int(input(" 日期 (DD): ")) + hour = int(input(" 小时 (HH, 24小时制): ")) + minute = int(input(" 分钟 (MM): ")) + second = int(input(" 秒 (SS): ")) + + dt = datetime(year, month, day, hour, minute, second) + print(f"\n使用的时间: {dt.strftime('%Y-%m-%d %H:%M:%S')} UTC") + return dt + except ValueError as e: + print(f"输入错误: {e}") + print("将使用默认时间 2024-06-15 10:00:00 UTC") + return datetime(2024, 6, 15, 10, 0, 0) + + +def calc_solar_angles(lat, lon, dt): + """ + 使用pvlib计算太阳天顶角和方位角。 + + Args: + lat (float): 纬度（度） + lon (float): 经度（度） + dt (datetime): datetime对象（UTC时间） + + Returns: + tuple: (太阳天顶角, 太阳方位角) - 单位：度 + 天顶角: 0=天顶, 90=地平线 + 方位角: 从北顺时针，0=北，90=东，180=南，270=西 + """ + # 创建时间索引 + times = pd.DatetimeIndex([dt]) + + # 使用pvlib计算太阳位置 + solpos = solarposition.get_solarposition(times, lat, lon) + + solar_zn_deg = solpos['zenith'].values[0] # 天顶角 + solar_az_deg = solpos['azimuth'].values[0] # 方位角 + + return solar_zn_deg, solar_az_deg + + +def write_envi_output(output_path, data, transform, crs, metadata=None): + """ + 将数据写入ENVI格式文件。 + + Args: + output_path (str): 输出文件路径 + data (numpy.ndarray): 数据数组 (bands, rows, cols) + transform: 地理变换参数 + crs: 坐标参考系 + metadata (dict, optional): 额外的元数据 + """ + with rasterio.open( + output_path, + 'w', + driver='ENVI', + height=data.shape[1], + width=data.shape[2], + count=data.shape[0], + dtype=data.dtype, + crs=crs, + transform=transform, + nodata=0.0, + ) as dst: + dst.write(data) + + # 写入ENVI头文件元数据 + if metadata: + dst.update_tags(**metadata) + + +def process_dem(dem_path, output_path, datetime_input=None, interactive=False): + """ + 主处理函数：读取DEM，计算坡度坡向，计算太阳角度，输出结果。 + + Args: + dem_path (str): DEM输入文件路径 + output_path (str): 输出文件路径 + datetime_input (str, optional): 日期时间字符串，如果为None则使用交互模式 + interactive (bool): 是否使用交互式时间输入 + + Returns: + dict: 处理结果信息 + """ + # ========== 1. 读取 DEM ========== + print(f"读取DEM文件: {dem_path}") + dem = xdem.DEM(dem_path) + + # ========== 2. 计算坡度、坡向 ========== + print("计算坡度和坡向...") + slope = xdem.terrain.slope(dem, method='ZevenbergThorne') + aspect = xdem.terrain.aspect(dem) + + # ========== 3. 获取空间参考信息 ========== + transform = dem.transform + crs = dem.crs + + # ========== 4. 转换为 numpy 数组 ========== + slope_data = slope.data + aspect_data = aspect.data + + # ========== 5. 获取DEM中心点坐标并计算太阳角度 ========== + print("\n获取DEM坐标信息...") + latitude, longitude = get_dem_center_coords(dem) + print(f" 中心点纬度: {latitude:.6f}°") + print(f" 中心点经度: {longitude:.6f}°") + print(f" 坐标系: WGS84") + + # 获取观测时间 + if datetime_input: + observation_time = parse_datetime(datetime_input) + print(f"\n使用命令行提供的时间: {observation_time.strftime('%Y-%m-%d %H:%M:%S')} UTC") + elif interactive or not sys.stdin.isatty(): + observation_time = get_user_input_datetime() + else: + # 默认时间（如果非交互式且未提供时间） + observation_time = datetime(2024, 6, 15, 10, 0, 0) + print(f"\n使用默认时间: {observation_time.strftime('%Y-%m-%d %H:%M:%S')} UTC") + + # 计算太阳角度 + print("\n计算太阳位置...") + solar_zn_deg, solar_az_deg = calc_solar_angles(latitude, longitude, observation_time) + + print(f" 太阳天顶角: {solar_zn_deg:.2f}°") + print(f" 太阳高度角: {90 - solar_zn_deg:.2f}°") + print(f" 太阳方位角: {solar_az_deg:.2f}°") + + # 转换为弧度 + solar_zn_rad = np.deg2rad(solar_zn_deg) + solar_az_rad = np.deg2rad(solar_az_deg) + + # 坡度和坡向也转换为弧度 + slope_rad = np.deg2rad(slope_data) + aspect_rad = np.deg2rad(aspect_data) + + # ========== 6. 计算 cosine_i ========== + print("\n计算入射角余弦 (cosine_i)...") + cosine_i = calc_cosine_i(solar_zn_rad, solar_az_rad, aspect_rad, slope_rad) + + # ========== 7. 处理无效值（NaN） ========== + slope_data = np.nan_to_num(slope_data, nan=0.0) + aspect_data = np.nan_to_num(aspect_data, nan=0.0) + cosine_i = np.nan_to_num(cosine_i, nan=0.0) + + # ========== 8. 堆叠为多波段数组 ========== + stacked = np.stack([slope_data, aspect_data, cosine_i], axis=0) + + # ========== 9. 准备元数据 ========== + metadata = { + 'description': 'Slope, Aspect, and Cosine_i derived from DEM', + 'band_names': 'slope(deg),aspect(deg),cosine_i', + 'solar_zenith_angle': str(solar_zn_deg), + 'solar_azimuth_angle': str(solar_az_deg), + 'solar_elevation_angle': str(90 - solar_zn_deg), + 'observation_time_utc': observation_time.strftime('%Y-%m-%d %H:%M:%S'), + 'dem_center_lat': str(latitude), + 'dem_center_lon': str(longitude), + } + + # ========== 10. 输出 ENVI 格式文件 ========== + print(f"\n保存结果到: {output_path}") + write_envi_output(output_path, stacked, transform, crs, metadata) + + print("\n处理完成!") + print(f" 输出文件: {output_path}") + print(f" 波段顺序: 1-坡度(°), 2-坡向(°), 3-cosine_i") + + return { + 'output_path': output_path, + 'solar_zenith': solar_zn_deg, + 'solar_azimuth': solar_az_deg, + 'center_lat': latitude, + 'center_lon': longitude, + 'observation_time': observation_time, + } + + +def main(): + """主函数 - 命令行入口""" + parser = argparse.ArgumentParser( + description='DEM地形处理工具 - 计算坡度、坡向和太阳入射角余弦', + formatter_class=argparse.RawDescriptionHelpFormatter, + epilog=''' +使用示例: + # 基本用法（交互式输入时间） + python slope_aspectV2.py -i dem.tif -o output.dat + + # 通过命令行指定时间 + python slope_aspectV2.py -i dem.tif -o output.dat -t "2024-06-15 10:30:00" + + # 仅指定日期（默认使用正午12:00） + python slope_aspectV2.py -i dem.tif -o output.dat -t "2024-06-15" + + # 显示详细处理信息 + python slope_aspectV2.py -i dem.tif -o output.dat -t "2024-06-15 10:30:00" -v + +支持的日期时间格式: + - "YYYY-MM-DD HH:MM:SS" (推荐) + - "YYYY-MM-DDTHH:MM:SS" (ISO格式) + - "YYYYMMDD_HHMMSS" + - "YYYY-MM-DD" (默认使用12:00:00) + ''' + ) + + parser.add_argument('-i', '--input', required=True, + help='输入DEM文件路径 (支持GeoTIFF等格式)') + + parser.add_argument('-o', '--output', required=True, + help='输出ENVI文件路径 (.dat)') + + parser.add_argument('-t', '--time', + help='观测日期时间 (UTC)，格式: "YYYY-MM-DD HH:MM:SS"。' + '如果不提供，将进入交互式输入模式') + + parser.add_argument('--interactive', action='store_true', + help='强制使用交互式时间输入模式') + + parser.add_argument('-v', '--verbose', action='store_true', + help='显示详细处理信息') + + parser.add_argument('--version', action='version', version='%(prog)s 2.0.0') + + args = parser.parse_args() + + # 验证输入文件存在 + if not Path(args.input).exists(): + print(f"错误: 输入文件不存在: {args.input}", file=sys.stderr) + sys.exit(1) + + # 创建输出目录（如果不存在） + output_dir = Path(args.output).parent + if not output_dir.exists(): + output_dir.mkdir(parents=True, exist_ok=True) + print(f"创建输出目录: {output_dir}") + + try: + # 执行处理 + result = process_dem( + dem_path=args.input, + output_path=args.output, + datetime_input=args.time, + interactive=args.interactive + ) + + if args.verbose: + print("\n详细结果:") + for key, value in result.items(): + print(f" {key}: {value}") + + sys.exit(0) + + except Exception as e: + print(f"\n处理失败: {e}", file=sys.stderr) + if args.verbose: + import traceback + traceback.print_exc() + sys.exit(1) + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/说明文档.docx b/说明文档.docx new file mode 100644 index 0000000..f261b82 Binary files /dev/null and b/说明文档.docx differ