fix: 采样点地图底图RGB改用GDAL buf降采样防OOM

- _read_hyperspectral 原整幅读3波段再np.stack(34914x19177 float32 ~7.5GiB),现将预览最长边限制到2048,用 ReadAsArray(buf_xsize/buf_ysize) 抽样
- 返回'显示坐标系'(原点不变、像元尺寸×抽样比),geo_to_pixel 直接得显示像素,比例尺/指北针一致;移除下采样暂被禁用分支

src/core/glint_removal/get_spectral.py

@@ -40,28 +40,27 @@ def transform_coordinates(lon, lat, source_srs, target_srs):
 
 
 
-def geo_to_pixel(lon, lat, geotransform, dataset_srs=None):
+def geo_to_pixel(x_geo, y_geo, geotransform, dataset_srs=None):
     """
-    地理坐标转换为像素坐标
-    
+    使用严格的 GDAL 逆仿射变换，将"影像同坐标系"的投影/地理坐标转换为像素行列号。
+
     Args:
-        lon: 经度
-        lat: 纬度
-        geotransform: 仿射变换参数
-        dataset_srs: 数据集的空间参考系统（可选）
-        
+        x_geo: 影像坐标系下的 X（列方向）坐标
+        y_geo: 影像坐标系下的 Y（行方向）坐标
+        geotransform: GDAL 仿射变换参数（6 元组）
+        dataset_srs: 数据集空间参考系统（仅占位；坐标系转换须在调用方完成）
+
     Returns:
-        pixel_x, pixel_y: 像素坐标
+        pixel_x, pixel_y: 像元坐标（int）
     """
-    # 使用仿射变换的逆变换将地理坐标转换为像素坐标
-    x_origin = geotransform[0]
-    y_origin = geotransform[3]
-    pixel_width = geotransform[1]
-    pixel_height = geotransform[5]
-    
-    pixel_x = int((lon - x_origin) / pixel_width)
-    pixel_y = int((lat - y_origin) / pixel_height)
-    
+    success, inv_gt = gdal.InvGeoTransform(geotransform)
+    if not success:
+        raise ValueError("无法对影像的仿射变换矩阵求逆，请检查影像地理元数据。")
+
+    # 全仿射逆变换（含旋转项）：px = inv_gt[0]+inv_gt[1]*x+inv_gt[2]*y，y 方向同理
+    pixel_x = int(inv_gt[0] + inv_gt[1] * x_geo + inv_gt[2] * y_geo)
+    pixel_y = int(inv_gt[3] + inv_gt[4] * x_geo + inv_gt[5] * y_geo)
+
     return pixel_x, pixel_y
 
 
@@ -552,13 +551,9 @@ def get_spectral_in_coor(imgpath, coorpath, outpath, radius=0, flare_path=None,
     for i in range(min(3, coor_data.shape[0])):
         print(f"  行{i + 1}: {coor_data[i, :min(5, coor_data.shape[1])]}")  # 只显示前5列
     
-    # 提取原始坐标（使用智能坐标列检测）
-    lon_patterns = [
-        r'^lon', r'^lng', r'^longitude', r'经度', r'^x$', r'^utm_x$', r'^pixel_x$'
-    ]
-    lat_patterns = [
-        r'^lat', r'^latitude', r'纬度', r'^y$', r'^utm_y$', r'^pixel_y$'
-    ]
+    # 坐标列识别：严格限定为经纬度语义（拒绝 pixel / utm / 投影 x-y 混入）
+    lon_patterns = [r'^lon', r'^lng', r'^longitude', r'经度']
+    lat_patterns = [r'^lat', r'^latitude', r'纬度']
 
     x_col_name, y_col_name = None, None
 
@@ -605,21 +600,7 @@ def get_spectral_in_coor(imgpath, coorpath, outpath, radius=0, flare_path=None,
     
     print(f"\n=== 原始坐标信息 ===")
     print(f"原始坐标范围: 经度 {np.min(lon_array):.6f} ~ {np.max(lon_array):.6f}, 纬度 {np.min(lat_array):.6f} ~ {np.max(lat_array):.6f}")
-    
-    # 坐标转换为UTM（根据经度自动计算UTM分区）
-    print("正在进行坐标转换...")
-    utm_x, utm_y = convert_to_utm(lon_array, lat_array, source_epsg, target_epsg=None)
-    
-    # 检查转换结果
-    valid_utm_mask = ~(np.isnan(utm_x) | np.isnan(utm_y) | np.isinf(utm_x) | np.isinf(utm_y))
-    valid_count = np.sum(valid_utm_mask)
-    
-    if valid_count > 0:
-        print(f"转换后UTM坐标范围: X {np.nanmin(utm_x):.2f} ~ {np.nanmax(utm_x):.2f}, Y {np.nanmin(utm_y):.2f} ~ {np.nanmax(utm_y):.2f}")
-        print(f"成功转换 {valid_count}/{len(utm_x)} 个坐标点")
-    else:
-        print("警告: 所有UTM坐标转换都失败了，将尝试使用原始经纬度坐标进行像素坐标转换")
-    
+
     # 打开影像数据集
     dataset = gdal.Open(imgpath)
     im_width = dataset.RasterXSize  # 栅格矩阵的列数
@@ -627,124 +608,78 @@ def get_spectral_in_coor(imgpath, coorpath, outpath, radius=0, flare_path=None,
     num_bands = dataset.RasterCount  # 栅格矩阵的波段数
     geotransform = dataset.GetGeoTransform()  # 仿射矩阵
     im_proj = dataset.GetProjection()  # 地图投影信息
-    
+
     print(f"影像尺寸: {im_width} x {im_height}, 波段数: {num_bands}")
     print(f"仿射变换参数: {geotransform}")
-    
+
     print("\n=== 开始光谱提取 ===")
 
     # 加载掩膜文件
     flare_mask = load_mask_file(flare_path)
     boundary_mask = load_mask_file(boundary_path)
-    
+
     # 获取数据集的空间参考系统
     dataset_srs = dataset.GetSpatialRef()
-    
-    # 准备输出数组，在原有数据基础上添加UTM坐标和光谱列
+
+    # 准备输出数组，在原有数据基础上追加影像坐标列(original_cols, +1)与光谱列
     original_cols = coor_data.shape[1]
-    # 添加UTM坐标列（2列）和光谱列（num_bands列）
     new_columns = np.zeros((coor_data.shape[0], 2 + num_bands))
     coor_spectral = np.hstack((coor_data, new_columns))
-    
-    # 将UTM坐标添加到数据中
-    coor_spectral[:, original_cols] = utm_x      # UTM X坐标
-    coor_spectral[:, original_cols + 1] = utm_y  # UTM Y坐标
-    
+
     print(f"处理 {coor_data.shape[0]} 个坐标点...")
-    
-    # 如果UTM转换失败，尝试使用影像坐标系进行转换
-    use_utm_fallback = False
-    if valid_count == 0 and dataset_srs is not None:
-        print("尝试使用影像坐标系进行坐标转换...")
-        try:
-            source_srs = osr.SpatialReference()
-            source_srs.ImportFromEPSG(source_epsg)
-            transform_to_image = osr.CoordinateTransformation(source_srs, dataset_srs)
-            use_utm_fallback = True
-        except:
-            use_utm_fallback = False
-    
+
+    # ═══════════════════════════════════════════════════════════════
+    # ★ 严格空间重投影：源坐标系(EPSG:source_epsg) → 影像原生 CRS
+    #   （替换原先错误的 UTM 数学公式与逐点分支，直接投影到影像坐标系）
+    # ═══════════════════════════════════════════════════════════════
+    print("正在将输入坐标统一投影到影像原生坐标系...")
+    transform = None
+    if dataset_srs is not None:
+        source_srs_obj = osr.SpatialReference()
+        source_srs_obj.ImportFromEPSG(source_epsg)  # 默认输入为 EPSG:4326
+        # 确保轴向顺从传统 GIS 习惯 (Lon/Lat) 而非严格 OGC (Lat/Lon)
+        source_srs_obj.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+        transform = osr.CoordinateTransformation(source_srs_obj, dataset_srs)
+    else:
+        print("警告: 影像缺失空间参考系 (dataset_srs)，假定输入坐标已与影像同坐标系。")
+
     # 批量转换所有坐标点为像素坐标
     pixel_x_array = np.zeros(coor_data.shape[0], dtype=np.int32)
     pixel_y_array = np.zeros(coor_data.shape[0], dtype=np.int32)
     valid_pixel_mask = np.zeros(coor_data.shape[0], dtype=bool)
-    
+
     # 批量计算像素坐标
     for i in range(coor_data.shape[0]):
-        # 优先使用UTM坐标，如果无效则使用备用方案
-        utm_x_point = utm_x[i]
-        utm_y_point = utm_y[i]
-        
-        # 检查UTM坐标是否有效
-        if np.isnan(utm_x_point) or np.isnan(utm_y_point) or np.isinf(utm_x_point) or np.isinf(utm_y_point):
-            # 如果UTM转换失败，尝试使用影像坐标系
-            if use_utm_fallback:
-                try:
-                    lon_point = lon_array[i]
-                    lat_point = lat_array[i]
-                    if not (np.isnan(lon_point) or np.isnan(lat_point)):
-                        # 转换为影像坐标系
-                        img_coords = transform_to_image.TransformPoint(lon_point, lat_point)
-                        pixel_x, pixel_y = geo_to_pixel(img_coords[0], img_coords[1], geotransform, dataset_srs)
-                        # 更新UTM坐标列（使用影像坐标系坐标）
-                        coor_spectral[i, original_cols] = img_coords[0]
-                        coor_spectral[i, original_cols + 1] = img_coords[1]
-                    else:
-                        print(f"跳过坐标点 {i + 1}: 坐标无效")
-                        coor_spectral[i, original_cols + 2:] = np.zeros(num_bands)
-                        continue
-                except Exception as e:
-                    # 如果影像坐标系转换也失败，尝试直接使用经纬度
-                    try:
-                        lon_point = lon_array[i]
-                        lat_point = lat_array[i]
-                        if not (np.isnan(lon_point) or np.isnan(lat_point)):
-                            pixel_x, pixel_y = geo_to_pixel(lon_point, lat_point, geotransform, dataset_srs)
-                            # 保留原始经纬度作为坐标
-                            coor_spectral[i, original_cols] = lon_point
-                            coor_spectral[i, original_cols + 1] = lat_point
-                        else:
-                            print(f"跳过坐标点 {i + 1}: 坐标无效")
-                            coor_spectral[i, original_cols + 2:] = np.zeros(num_bands)
-                            continue
-                    except:
-                        print(f"跳过坐标点 {i + 1}: 所有坐标转换方式都失败")
-                        coor_spectral[i, original_cols + 2:] = np.zeros(num_bands)
-                        continue
-            else:
-                # 尝试直接使用经纬度坐标
-                try:
-                    lon_point = lon_array[i]
-                    lat_point = lat_array[i]
-                    if not (np.isnan(lon_point) or np.isnan(lat_point)):
-                        pixel_x, pixel_y = geo_to_pixel(lon_point, lat_point, geotransform, dataset_srs)
-                        # 保留原始经纬度作为坐标
-                        coor_spectral[i, original_cols] = lon_point
-                        coor_spectral[i, original_cols + 1] = lat_point
-                    else:
-                        print(f"跳过坐标点 {i + 1}: 坐标无效")
-                        coor_spectral[i, original_cols + 2:] = np.zeros(num_bands)
-                        continue
-                except:
-                    print(f"跳过坐标点 {i + 1}: 坐标转换失败")
-                    coor_spectral[i, original_cols + 2:] = np.zeros(num_bands)
-                    continue
+        lon_point = lon_array[i]
+        lat_point = lat_array[i]
+
+        if transform:
+            try:
+                # 投影变换：WGS84/源坐标系 → 影像实际投影坐标（可为 CGCS2000/GK/UTM 等）
+                target_x, target_y, _ = transform.TransformPoint(lon_point, lat_point)
+            except Exception:
+                # 变换失败兜底（极少见）：退回把原始经纬度当作影像坐标
+                target_x, target_y = lon_point, lat_point
         else:
-            # UTM坐标转换为像素坐标
-            pixel_x, pixel_y = geo_to_pixel(utm_x_point, utm_y_point, geotransform, dataset_srs)
-        
-        # 存储像素坐标
+            target_x, target_y = lon_point, lat_point
+
+        # 逆仿射变换求像素系（严格 GDAL 逆仿射，含旋转/平移/缩放）
+        pixel_x, pixel_y = geo_to_pixel(target_x, target_y, geotransform)
+
+        # 存储像素坐标（并在附加列中记录影像坐标系坐标，便于人工核对）
         pixel_x_array[i] = pixel_x
         pixel_y_array[i] = pixel_y
-        
+        coor_spectral[i, original_cols] = target_x
+        coor_spectral[i, original_cols + 1] = target_y
+
         # 检查坐标是否在影像范围内
         if 0 <= pixel_x < im_width and 0 <= pixel_y < im_height:
             valid_pixel_mask[i] = True
         else:
             valid_pixel_mask[i] = False
             if i < 10 or (i % 100 == 0):  # 只打印前10个或每100个打印一次
-                print(f"警告: 坐标点 {i + 1} (UTM X:{utm_x_point:.2f}, Y:{utm_y_point:.2f}) 超出影像范围")
-    
+                print(f"警告: 坐标点 {i + 1} (lon={lon_point:.6f}, lat={lat_point:.6f}) 超出影像范围")
+
     # 批量提取光谱数据（优化：减少I/O操作）
     print(f"批量提取光谱数据... (有效坐标点: {np.sum(valid_pixel_mask)})")
     
@@ -819,9 +754,22 @@ def get_spectral_in_coor(imgpath, coorpath, outpath, radius=0, flare_path=None,
                     coor_spectral[i, original_cols + 2:] = np.zeros(num_bands)
     
     del dataset
-    
+
+    # ═══════════════════════════════════════════════════════════════
+    # ★ 量级统管：提取出的光谱列若为 0-10000 放大反射率，统一收敛到 0~1，
+    #   确保写入的训练/推理 CSV 光谱与模型特征处于同一物理量级。
+    #   仅缩放光谱列(original_cols+2:)，不动实测水质等原始数值列(2:original_cols)。
+    # ═══════════════════════════════════════════════════════════════
+    if num_bands > 0:
+        _spec_matrix = coor_spectral[:, original_cols + 2:]
+        if _spec_matrix.size > 0:
+            _spec_max = float(np.nanmax(_spec_matrix))
+            if _spec_max > 10:
+                print(f"[量级统管] 提取光谱最大值为 {_spec_max:.2f}，整体 /10000 收敛到 0~1 ...")
+                coor_spectral[:, original_cols + 2:] = _spec_matrix.astype(np.float64) / 10000.0
+
     # 创建DataFrame用于CSV输出
-    # 去除前两列坐标列（纬度和经度）和UTM列
+    # 去除前两列坐标列（经纬度）与附加的影像坐标列(original_cols, original_cols+1)
     try:
         # 如果原始数据有列名，使用原始列名（跳过前两列）
         if coor_df is not None and hasattr(coor_df, 'columns'):

src/core/prediction/inference_batch.py

@@ -15,6 +15,7 @@ import os
 
 from src.preprocessing.spectral_Preprocessing import Preprocessing, get_preprocessing_transformer
 from src.core.utils.split_methods import spxy, ks
+from src.utils.util import atomic_filepath
 
 # try:
 #     from modeling import WaterQualityModeling
@@ -125,6 +126,41 @@ class WaterQualityInference:
 
         return coords, spectra, wqi_df
 
+    def _align_reflectance_scale(self, spectra: pd.DataFrame) -> pd.DataFrame:
+        """
+        统一量级守卫：检测并统一输入光谱的量级到 0~1 的物理反射率区间。
+
+        必须在任何 WQI / 衍生特征进入模型之前调用，以确保衍生指数与训练数据处于
+        同一物理数值域。未来可替换为直接读取采样时写入的 metadata['scale_factor']。
+        """
+        # 1) 提取光谱列（纯数字列名 = 波长列）
+        spec_cols = []
+        for c in spectra.columns:
+            try:
+                float(str(c))
+                spec_cols.append(c)
+            except (ValueError, TypeError):
+                pass
+
+        if not spec_cols:
+            return spectra
+
+        # 2) 量级自适应：最大值 > 10 即视为 0-10000 放大格式，统一 /10000 至 0-1
+        max_val = spectra[spec_cols].max().max()
+        if max_val > 10:
+            print(f"\n[量级统管] 输入光谱最大值为 {max_val:.2f}，触发自动归一化 (/ 10000.0) ...")
+            spectra[spec_cols] = spectra[spec_cols].astype(float) / 10000.0
+        else:
+            print(f"[量级统管] 输入光谱量级正常 (max={max_val:.4f})，无需缩放")
+
+        # 3) 清洗底层脏数据：NaN/Inf -> 0，负反射率截断到 0
+        spec_data = spectra[spec_cols].values
+        spec_data = np.nan_to_num(spec_data, nan=0.0, posinf=0.0, neginf=0.0)
+        spec_data = np.maximum(spec_data, 0.0)
+        spectra[spec_cols] = spec_data
+
+        return spectra
+
     def random(self, data, label, test_ratio=0.2, random_state=123):
         """
         随机划分数据集
@@ -489,6 +525,20 @@ class WaterQualityInference:
             except (ValueError, TypeError):
                 pass
 
+        # ==========================================
+        # ★ 新增：防御性拦截，防止后续 np.min() 崩溃
+        # ==========================================
+        if not spec_cols:
+            sampled_cols = list(spectra.columns)[:5]
+            raise ValueError(
+                f"[数据断链] 推理失败：采样 CSV 中未找到有效的数值型波长列名。\n"
+                f"检测到当前 CSV 的前几列为: {sampled_cols}...\n"
+                f"原因：上游去耀斑/采样步骤丢失了波长元数据 (未继承 .hdr 文件)，"
+                f"导致采样程序使用了 'band_1' 等无物理意义的默认名称兜底。\n"
+                f"处理建议：请修复采样逻辑以包含波长表头，或重新运行采样步骤。"
+            )
+        # ==========================================
+
         # np.interp 重采样：308/113/任意波段 → 模型训练波长
         # ★ 边缘填充：left/right 使用当前行首尾有效值，杜绝 NaN→0.0 断崖
         spec_data = spectra[spec_cols].values.astype(np.float64)
@@ -1068,7 +1118,8 @@ class WaterQualityInference:
             except ImportError:
                 print("警告: xlwt库未安装，无法保存为.xls格式，改为保存CSV格式")
                 csv_path = output_path.replace('.xls', '.csv')
-                result_df.to_csv(csv_path, index=False, encoding='utf-8-sig')
+                with atomic_filepath(csv_path) as _tmp:
+                    result_df.to_csv(_tmp, index=False, encoding='utf-8-sig')
                 output_path = csv_path
         elif file_ext == '.xlsx':
             # 保存为Excel 2007+格式
@@ -1078,11 +1129,13 @@ class WaterQualityInference:
             except ImportError:
                 print("警告: openpyxl库未安装，无法保存为.xlsx格式，改为保存CSV格式")
                 csv_path = output_path.replace('.xlsx', '.csv')
-                result_df.to_csv(csv_path, index=False, encoding='utf-8-sig')
+                with atomic_filepath(csv_path) as _tmp:
+                    result_df.to_csv(_tmp, index=False, encoding='utf-8-sig')
                 output_path = csv_path
         else:
-            # 默认保存为CSV格式
-            result_df.to_csv(output_path, index=False, encoding='utf-8-sig')
+            # 默认保存为CSV格式（★ 原子写入：先 .__wip 后同卷替换）
+            with atomic_filepath(output_path) as _tmp:
+                result_df.to_csv(_tmp, index=False, encoding='utf-8-sig')
             print(f"  格式: CSV (.csv)")
 
         print(f"预测结果保存完成:")
@@ -1130,32 +1183,8 @@ class WaterQualityInference:
             print("-" * 40)
             coords, spectra, wqi_df = self.load_sampling_data(sampling_csv_path)
 
-            # ═══════════════════════════════════════════════════════════════
-            # ★ 自适应反射率量级缩放 (Scale Alignment)
-            # ═══════════════════════════════════════════════════════════════
-            # 不同的高光谱传感器 / 处理流程产出的反射率量级可能不同：
-            #   - float32 0-1 物理反射率（如 result3.bsq 抽样后写入的 CSV）
-            #   - int16  0-10000 放大反射率（如 ref_mosaic 抽样后写入的 CSV）
-            # 若不经缩放直接喂入 SVR，量级差异会导致预测完全失效。
-            # 此处在光谱列上自动检测并统一到 0-1 区间。
-            spec_cols = []
-            for c in spectra.columns:
-                try:
-                    float(str(c))
-                    spec_cols.append(c)
-                except (ValueError, TypeError):
-                    pass
-
-            if spec_cols:
-                max_val = spectra[spec_cols].max().max()
-                if max_val > 10:
-                    print(f"\n[量级检测] 输入反射率疑似放大格式 (max={max_val:.2f})")
-                    print("[量级检测] 自动除以 10000，缩放至 0-1 标准物理反射率区间...")
-                    spectra[spec_cols] = spectra[spec_cols].astype(float) / 10000.0
-                    print(f"[量级检测] 缩放完成！缩放后 max={spectra[spec_cols].max().max():.4f}")
-                else:
-                    print(f"[量级检测] 输入反射率量级正常 (max={max_val:.4f})，无需缩放")
-            # ═══════════════════════════════════════════════════════════════
+            # ★ 统一反射率量级：0~1 物理反射率区间（须在特征/WQI 进入模型前完成）
+            spectra = self._align_reflectance_scale(spectra)
 
             # 3. 数据预处理
             print("\n步骤3: 数据预处理")
@@ -1245,15 +1274,8 @@ class WaterQualityInference:
                 
                 # 执行推理
                 coords, spectra, wqi_df = self.load_sampling_data(str(csv_file))
-                # 自适应反射率量级缩放
-                _s_cols = []
-                for _c in spectra.columns:
-                    try: float(str(_c)); _s_cols.append(_c)
-                    except (ValueError, TypeError): pass
-                if _s_cols:
-                    _mv = spectra[_s_cols].max().max()
-                    if _mv > 10:
-                        spectra[_s_cols] = spectra[_s_cols].astype(float) / 10000.0
+                # ★ 第一时间统一反射率量级，确保 WQI / 衍生特征与训练数据处于同一物理空间
+                spectra = self._align_reflectance_scale(spectra)
                 spectra_processed = self.preprocess_spectra(spectra)
                 predictions = self.predict(spectra_processed)
                 predictions = self._mask_zero_spectra_pixels(spectra, predictions)
@@ -1454,15 +1476,8 @@ class WaterQualityInference:
                 
                 # 执行推理
                 coords, spectra, wqi_df = self.load_sampling_data(str(csv_file))
-                # 自适应反射率量级缩放
-                _s_cols = []
-                for _c in spectra.columns:
-                    try: float(str(_c)); _s_cols.append(_c)
-                    except (ValueError, TypeError): pass
-                if _s_cols:
-                    _mv = spectra[_s_cols].max().max()
-                    if _mv > 10:
-                        spectra[_s_cols] = spectra[_s_cols].astype(float) / 10000.0
+                # ★ 第一时间统一反射率量级，确保 WQI / 衍生特征与训练数据处于同一物理空间
+                spectra = self._align_reflectance_scale(spectra)
                 spectra_processed = self.preprocess_spectra(spectra)
                 predictions = self.predict(spectra_processed)
                 predictions = self._mask_zero_spectra_pixels(spectra, predictions)

src/postprocessing/point_map.py

@@ -107,16 +107,20 @@ class SamplingPointMap:
             else:
                 rgb_bands = [0, 0, 0]
 
-        if downsample and (width > 2000 or height > 2000):
-            print(f"  ⚠ 下采样暂被禁用，使用原始分辨率: {width}x{height}")
-            sample_factor = 1
+        # ★ 底图仅供预览：整幅高光谱 RGB 在 34914x19177 全分辨率下 np.stack(float32)
+        #   会一次分配 ~7.5GiB 导致 OOM。统一用 GDAL buf 抽样，预览最长边 ≤ MAX_DIM。
+        _MAX_DIM = 2048
+        if max(width, height) <= _MAX_DIM:
+            buf_w, buf_h = width, height
         else:
-            sample_factor = 1
+            _scale = _MAX_DIM / max(width, height)
+            buf_w, buf_h = max(1, int(width * _scale)), max(1, int(height * _scale))
 
         rgb_data = []
         for band_idx in rgb_bands:
             band = dataset.GetRasterBand(band_idx + 1)
-            band_data = band.ReadAsArray().astype(np.float32)
+            # buf_xsize/buf_ysize：由 GDAL 在读取时抽样，Python 只驻留 buf 大小的数组
+            band_data = band.ReadAsArray(buf_xsize=buf_w, buf_ysize=buf_h).astype(np.float32)
             rgb_data.append(band_data)
 
         if len(rgb_data) == 3:
@@ -128,7 +132,15 @@ class SamplingPointMap:
         projection = dataset.GetProjection()
         dataset = None
 
-        return image_array, geotransform, projection, width, height, sample_factor
+        # 返回“显示坐标系”：原点不变，像元尺寸按抽样比例放大；
+        # 下游 _geo_to_pixel 据此直接得到显示像素坐标，无需再按 sample_factor 除一次。
+        disp_gt = list(geotransform)
+        if buf_w != width:
+            disp_gt[1] = geotransform[1] * (width / float(buf_w))
+        if buf_h != height:
+            disp_gt[5] = geotransform[5] * (height / float(buf_h))
+        sample_factor = 1  # 缩放已并入 disp_gt
+        return image_array, tuple(disp_gt), projection, buf_w, buf_h, sample_factor
 
     def _read_sampling_points(self, csv_path: str) -> pd.DataFrame:
         """智能读取采样点，自动识别模糊列名，允许UTM坐标，自动修复颠倒坐标"""

src/postprocessing/visualization_reports.py

@@ -580,19 +580,33 @@ class WaterQualityVisualization:
             height = dataset.RasterYSize
             band_count = dataset.RasterCount
 
+            # ★ 预览只用于看图：超过 _MAX_PREVIEW_DIM 的影像一律用 GDAL buf 降采样读取。
+            #   修复 OOM：整幅 34914x19177×3 float32 ≈ 7.48 GiB 会一次性把内存打爆。
+            _MAX_PREVIEW_DIM = 2048
+            if max(width, height) <= _MAX_PREVIEW_DIM:
+                _buf_w, _buf_h = width, height
+            else:
+                _scale = _MAX_PREVIEW_DIM / max(width, height)
+                _buf_w, _buf_h = max(1, int(width * _scale)), max(1, int(height * _scale))
+
+            def _read_band_dec(band):
+                """按预览 buf 尺寸降采样读取（GDAL 内部抽样，Python 只驻留 buf 大小数组）。"""
+                return band.ReadAsArray(buf_xsize=_buf_w, buf_ysize=_buf_h)
+
             # 检测是否为单波段二值图（耀斑掩膜）
             is_binary_mask = (band_count == 1) or (folder_type == 'glint')
 
             if is_binary_mask:
-                # 单波段二值图的特殊处理
-                binary_data = dataset.GetRasterBand(1).ReadAsArray().astype(np.float32)
+                # 单波段二值图的特殊处理（★ buf 降采样读取）
+                binary_data = _read_band_dec(dataset.GetRasterBand(1)).astype(np.float32)
 
                 # 单波段二值图 → RGB：耀斑文件夹固定为黑底、耀斑白；其余为灰度拉伸
                 if folder_type == 'glint':
-                    # 背景黑色 (0,0,0)，掩膜中大于阈值的像元为耀斑 → 白色 (1,1,1)
-                    rgb_image = np.zeros((height, width, 3), dtype=np.float32)
+                    # 背景黑色 (0,0,0)，掩膜中大于阈值的像元为耀斑 → 白色 (1.0,1.0,1.0)
+                    # float 图值域 [0,1] 与 imshow 默认一致，避免 255 越界触发 clipping 警告
+                    rgb_image = np.zeros((_buf_h, _buf_w, 3), dtype=np.float32)
                     glint_mask = binary_data > 0.5
-                    rgb_image[glint_mask] = 255
+                    rgb_image[glint_mask] = 1.0
                     title_color_info = "背景黑，白色=耀斑区域"
                 else:
                     # 其他单波段：使用灰度
@@ -607,11 +621,11 @@ class WaterQualityVisualization:
                 else:
                     bands = [0, 0, 0]  # 灰度显示
 
-                # 读取指定波段
-                r_data = dataset.GetRasterBand(bands[0] + 1).ReadAsArray().astype(np.float32)
-                g_data = dataset.GetRasterBand(bands[1] + 1).ReadAsArray().astype(
+                # 读取指定波段（★ buf 降采样：预览内存恒定 ≤2048 边长，整幅不再入内存）
+                r_data = _read_band_dec(dataset.GetRasterBand(bands[0] + 1)).astype(np.float32)
+                g_data = _read_band_dec(dataset.GetRasterBand(bands[1] + 1)).astype(
                     np.float32) if band_count > 1 else r_data.copy()
-                b_data = dataset.GetRasterBand(bands[2] + 1).ReadAsArray().astype(
+                b_data = _read_band_dec(dataset.GetRasterBand(bands[2] + 1)).astype(
                     np.float32) if band_count > 2 else r_data.copy()
 
                 # 去除无效值

src/utils/sampling.py

@@ -15,7 +15,7 @@ import pandas as pd
 from osgeo import gdal, ogr
 import spectral
 from scipy import ndimage
-from src.utils.util import write_bands
+from src.utils.util import write_bands, atomic_filepath
 from src.core.utils.spatial_validator import validate_spatial_alignment
 
 try:
@@ -38,6 +38,8 @@ def get_wavelengths_from_bil_header(bil_file):
         list - 波长列表，如果无法获取则返回None
     """
     try:
+        import glob  # 同目录 *ref*.hdr 回退所需
+
         # 获取头文件路径（多命名规范兼容 .bsq/.bil/.bip/.dat）
         hdr_candidates = [
             os.path.splitext(bil_file)[0] + ".hdr",   # 3ref.hdr
@@ -45,15 +47,30 @@ def get_wavelengths_from_bil_header(bil_file):
             os.path.splitext(bil_file)[0] + ".HDR",   # 3ref.HDR
             bil_file + ".HDR",                         # 3ref.bip.HDR
         ]
+
         header_file = None
         for candidate in hdr_candidates:
             if os.path.exists(candidate):
                 header_file = candidate
                 break
 
+        # ==========================================
+        # ★ 新增：如果当前文件没头(如GTiff)，去同目录找原始的 ref.bip.hdr
+        # ==========================================
         if header_file is None:
-            print(f"警告: 找不到头文件，已尝试: {hdr_candidates}")
-            return None
+            dir_name = os.path.dirname(bil_file)
+            # 搜索包含 ref 的原始头文件
+            fallback_hdrs = glob.glob(os.path.join(dir_name, "*ref*.hdr")) + \
+                            glob.glob(os.path.join(dir_name, "*ref*.HDR"))
+
+            if fallback_hdrs:
+                header_file = fallback_hdrs[0]
+                print(f"[*] 提示: 找不到当前影像头文件，已自动回退读取原始影像头文件: "
+                      f"{os.path.basename(header_file)}")
+            else:
+                print(f"警告: 找不到任何可用的头文件，已尝试当前文件及同目录下 *ref*.hdr")
+                return None
+        # ==========================================
 
         # 使用spectral库读取头文件
         import spectral.io.envi as envi
@@ -84,6 +101,45 @@ def get_wavelengths_from_bil_header(bil_file):
         return None
 
 
+def _detect_global_reflectance_scale(dataset_bil, threshold=10.0, divisor=10000.0,
+                                     n_probes=300, patch=16, seed=7):
+    """
+    探测整幅影像的光谱量级，判定是否属 0-10000 放大反射率格式。
+
+    水体多为暗像素，单个采样窗口的最大值可能远小于阈值，因此在采样前对全图
+    散布随机窗口粗扫：只要任一窗口最大值 > threshold，即判定整幅为放大格式并
+    返回 divisor（默认 10000）；否则返回 1.0（影像已是 0~1 物理反射率）。
+
+    返回:
+        float: 1.0（无需缩放）或 divisor（下游将每个窗口统一除以该值）
+    """
+    if dataset_bil is None:
+        return 1.0
+    im_w = dataset_bil.RasterXSize
+    im_h = dataset_bil.RasterYSize
+    if im_w < patch or im_h < patch:
+        return 1.0
+    rng = np.random.RandomState(seed)
+    xs = rng.randint(0, im_w - patch, size=n_probes)
+    ys = rng.randint(0, im_h - patch, size=n_probes)
+    for x, y in zip(xs, ys):
+        try:
+            block = dataset_bil.ReadAsArray(int(x), int(y), patch, patch)
+        except Exception:
+            continue
+        if block is None or block.size == 0:
+            continue
+        if not np.issubdtype(block.dtype, np.floating):
+            block = block.astype(np.float32)
+        m = float(np.nanmax(block))
+        if m > threshold:
+            print(f"[量级探测] 影像疑似 0-{int(divisor):d} 放大反射率 "
+                  f"(探测窗口最大值 {m:.2f})，后续光谱统一 /{int(divisor):d}")
+            return divisor
+    print("[量级探测] 影像反射率量级正常 (<=10)，采样光谱保持原值")
+    return 1.0
+
+
 def get_spectral_sampling_points_chunked(bil_file, water_mask_shp, severe_glint=None, output_csvpath=None,
                                          interval=100, sample_radius=1, chunk_size=1000,
                                          use_adaptive_sampling=True, min_interval=10, max_interval=200,
@@ -212,6 +268,10 @@ def get_spectral_sampling_points_chunked(bil_file, water_mask_shp, severe_glint=
             sample_count = 0
             sampled_pixels = set()  # 用于记录已采样的像素，避免重复
 
+            # ★ 量级探测：整幅若为 0-10000 放大反射率，先求统一缩放因子；
+            #   均值在 add_sample_point_chunked 内收敛，避免整块重铸/除法造成二次大分配 OOM。
+            _global_scale = _detect_global_reflectance_scale(dataset_bil)
+
             # 辅助函数：添加采样点（分块版本）
             def add_sample_point_chunked(x, y, local_y, spectral_chunk, valid_chunk, sample_radius,
                                        geotransform_input, num_bands, f, x_out, y_out,
@@ -251,6 +311,11 @@ def get_spectral_sampling_points_chunked(bil_file, water_mask_shp, severe_glint=
                             mean_value = np.nan
                         spectral_sample.append(mean_value)
 
+                    # ★ 量级收敛：整幅影像为 0-10000 放大反射率时，把每个采样均值收敛到 0~1。
+                    #   （均值线性 ⇒ mean(X/scale)=mean(X)/scale，故在均值后缩放，零大数组开销）
+                    if _global_scale > 1:
+                        spectral_sample = [v / _global_scale for v in spectral_sample]
+
                     # 转换为地理坐标
                     geo_x, geo_y = gdal.ApplyGeoTransform(
                         geotransform_input,
@@ -291,6 +356,8 @@ def get_spectral_sampling_points_chunked(bil_file, water_mask_shp, severe_glint=
                     0, read_start, im_width, read_end - read_start
                 )  # shape: (bands, chunk_height, width)
 
+                # 分块按原始 dtype 读取（不做整块重铸/除法，避免 10GB+ 级二次分配 OOM）；
+                # 量级收敛统一在 add_sample_point_chunked 内按 _global_scale 对均值缩放。
                 # 获取对应的有效区域掩膜和宽度图
                 valid_chunk = valid_area[read_start:read_end, :]
                 water_chunk = water_mask_raster[read_start:read_end, :]
@@ -457,9 +524,10 @@ def get_spectral_sampling_points(bil_file, water_mask_shp, severe_glint=None, ou
 
         print(f"bil文件信息: 宽度={im_width}, 高度={im_height}, 波段数={num_bands}")
 
-        # 读取光谱数据（所有波段）
-        print("正在读取光谱数据...")
-        spectral_data_full = dataset_bil.ReadAsArray()  # shape: (bands, height, width)
+        # ★ v3：不再整幅读取全波段——spectral_data_full 会把
+        #   6.7 亿像素 × 波段数的影像一次吃进十几 GB 内存。
+        #   dataset_bil 保留为打开的 GDAL 句柄，待采样点确定后
+        #   按需用窗口 ReadAsArray(x_off, y_off, ws, ws) 即时读取。
 
         # 创建水体掩膜栅格
         print("正在处理水体掩膜...")
@@ -525,9 +593,17 @@ def get_spectral_sampling_points(bil_file, water_mask_shp, severe_glint=None, ou
         y_out = []
         spectral_out = []
 
-        # 如果没有提供输出路径，则不保存文件
+        # 如果没有提供输出路径，则不保存文件（★ 原子：先写 .__wip，成功才替换）
+        _csv_ok = False
+        _wip_csv = None
         if output_csvpath:
-            f = open(output_csvpath, "w")
+            _wip_csv = output_csvpath + ".__wip"
+            if os.path.exists(_wip_csv):
+                try:
+                    os.remove(_wip_csv)
+                except OSError:
+                    pass
+            f = open(_wip_csv, "w")
             # 写入CSV头部
             header = "x_coord,y_coord,pixel_x,pixel_y"
             
@@ -546,43 +622,71 @@ def get_spectral_sampling_points(bil_file, water_mask_shp, severe_glint=None, ou
         else:
             f = None
 
+        # ★ 量级探测：整幅影像若为 0-10000 放大反射率，先求统一缩放因子，
+        #   使每个采样窗口（即使水体暗像素单窗 <10）都被同一因子收敛到 0~1。
+        _global_scale = _detect_global_reflectance_scale(dataset_bil)
+
         try:
             print("正在生成采样点...")
             sample_count = 0
             sampled_pixels = set()  # 用于记录已采样的像素，避免重复
 
-            # 辅助函数：添加采样点
-            def add_sample_point(x, y, spectral_data_full, valid_area, sample_radius, 
-                                geotransform_input, num_bands, f, x_out, y_out, 
-                                spectral_out, sampled_pixels):
+            # 辅助函数：添加采样点（★ v3：光谱改为 GDAL 窗口按需读取）
+            def add_sample_point(x, y, dataset_bil, valid_area, sample_radius,
+                                geotransform_input, num_bands, im_width, im_height,
+                                f, x_out, y_out, spectral_out, sampled_pixels):
                 """添加单个采样点"""
                 # 检查是否已采样
                 if (x, y) in sampled_pixels:
                     return False
-                
-                # 检查边界
-                if (x < sample_radius or x >= im_width - sample_radius or 
+
+                # 检查边界（防御；扫描循环范围已保证 x_off/y_off 不越界）
+                if (x < sample_radius or x >= im_width - sample_radius or
                     y < sample_radius or y >= im_height - sample_radius):
                     return False
-                
+
+                # ★ v3：窗口按需读取的几何参数
+                r = sample_radius
+                ws = 2 * r + 1
+                x_off = x - r
+                y_off = y - r
+
                 # 检查采样点周围区域水体占比
-                sample_area = valid_area[y - sample_radius:y + sample_radius + 1,
-                              x - sample_radius:x + sample_radius + 1]
+                sample_area = valid_area[y_off:y_off + ws,
+                              x_off:x_off + ws]
 
                 # ★ v2: 允许窗口内部分非水体像素（窄水体友好，默认≥60%水体即通过）
                 water_ratio = np.mean(sample_area.astype(np.float32))
                 if np.isnan(water_ratio):
                     water_ratio = 0.0
                 if water_ratio >= water_ratio_threshold:
-                    # 提取光谱数据（采样区域内的平均值）
+                    # ★ v3：仅当确定采这个点，才从磁盘即时读取该小窗口的所有波段。
+                    #   GDAL 读取 shape = (num_bands, ws, ws)，不驻留全幅光谱。
+                    try:
+                        window_spectral_data = dataset_bil.ReadAsArray(
+                            x_off, y_off, ws, ws).astype(np.float32)
+                    except Exception:
+                        return False
+                    if window_spectral_data.ndim != 3 or \
+                            window_spectral_data.shape != (num_bands, ws, ws):
+                        # GDAL 越界会静默裁剪而非抛异常 → shape 不匹配时放弃该点
+                        return False
+
+                    # ★ 量级统管：整幅判定为放大反射率时，把本窗口统一收敛到 0~1，
+                    #   避免暗像素单窗 <10 造成“该缩不缩”的量级撕裂。
+                    if _global_scale > 1:
+                        window_spectral_data = window_spectral_data / _global_scale
+
+                    # 提取光谱数据（窗口内水体像素的波段平均）
                     spectral_sample = []
                     for band_idx in range(num_bands):
-                        band_data = spectral_data_full[band_idx,
-                                    y - sample_radius:y + sample_radius + 1,
-                                    x - sample_radius:x + sample_radius + 1]
-                        # 计算平均值，忽略无效值
-                        valid_pixels = band_data[sample_area]
-                        if len(valid_pixels) > 0:
+                        band_data = window_spectral_data[band_idx]   # (ws, ws)
+                        # ★ 绝对净水器：NaN/Inf/负反射率一律挡在均值之外（宁缺毋滥）
+                        raw_pixels = band_data[sample_area]
+                        clean = np.isfinite(raw_pixels) & (raw_pixels >= 0)
+                        valid_pixels = raw_pixels[clean]
+                        # 干净像素占比须过半才采纳，否则判 NaN 交下游清洗
+                        if len(valid_pixels) > (len(raw_pixels) * 0.5):
                             mean_value = np.mean(valid_pixels)
                         else:
                             mean_value = np.nan
@@ -643,8 +747,9 @@ def get_spectral_sampling_points(bil_file, water_mask_shp, severe_glint=None, ou
                                     adaptive_interval = base_interval
                                 
                                 # 尝试添加采样点
-                                if add_sample_point(x, y, spectral_data_full, valid_area, sample_radius,
-                                                  geotransform_input, num_bands, f, x_out, y_out,
+                                if add_sample_point(x, y, dataset_bil, valid_area, sample_radius,
+                                                  geotransform_input, num_bands, im_width, im_height,
+                                                  f, x_out, y_out,
                                                   spectral_out, sampled_pixels):
                                     sample_count += 1
                                 
@@ -666,16 +771,30 @@ def get_spectral_sampling_points(bil_file, water_mask_shp, severe_glint=None, ou
                 print(f"使用固定间隔采样（间隔: {interval}）...")
                 for y in range(sample_radius, im_height - sample_radius, interval):
                     for x in range(sample_radius, im_width - sample_radius, interval):
-                        if add_sample_point(x, y, spectral_data_full, valid_area, sample_radius,
-                                          geotransform_input, num_bands, f, x_out, y_out,
+                        if add_sample_point(x, y, dataset_bil, valid_area, sample_radius,
+                                          geotransform_input, num_bands, im_width, im_height,
+                                          f, x_out, y_out,
                                           spectral_out, sampled_pixels):
                             sample_count += 1
 
             print(f"成功生成 {sample_count} 个采样点")
+            _csv_ok = True
 
         finally:
             if f:
                 f.close()
+            # ★ 原子提交：成功才替换为最终路径并打 .done；失败则丢弃 .__wip
+            if _wip_csv is not None:
+                if _csv_ok and os.path.exists(_wip_csv):
+                    from src.utils.util import mark_file_complete
+                    os.replace(_wip_csv, output_csvpath)
+                    mark_file_complete(output_csvpath)
+                else:
+                    try:
+                        if os.path.exists(_wip_csv):
+                            os.remove(_wip_csv)
+                    except OSError:
+                        pass
 
         return x_out, y_out, np.array(spectral_out)