From 4d49753cca6d9ba1017c5e5f803039867888915d Mon Sep 17 00:00:00 2001 From: duxin Date: Tue, 8 Sep 2026 17:26:18 +0800 Subject: [PATCH] =?UTF-8?q?fix:=20=E5=AE=9E=E6=B5=8B=E7=AB=99=E5=9D=90?= =?UTF-8?q?=E6=A0=87=E6=94=B9osr=E9=87=8D=E6=8A=95=E5=BD=B1+GDAL=E9=80=86?= =?UTF-8?q?=E4=BB=BF=E5=B0=84,=E6=8F=90=E5=8F=96=E5=85=89=E8=B0=B1?= =?UTF-8?q?=E5=86=99CSV=E5=89=8D=E7=BB=9F=E4=B8=800~1=E9=87=8F=E7=BA=B2?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - geo_to_pixel 改用 gdal.InvGeoTransform 全仿射(含旋转项);坐标列正则收紧,只认经纬度,剔除 pixel/utm 误判 - get_spectral_in_coor 主流程由错误 UTM 数学公式改为 osr.CoordinateTransformation(EPSG源->影像原生CRS) 重投影 - 提取光谱列写 CSV 前做 >10 则 /10000 量级收敛(不动实测值列) --- src/core/glint_removal/get_spectral.py | 212 ++++++++++--------------- 1 file changed, 80 insertions(+), 132 deletions(-) diff --git a/src/core/glint_removal/get_spectral.py b/src/core/glint_removal/get_spectral.py index 02c514e..5963b02 100644 --- a/src/core/glint_removal/get_spectral.py +++ b/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'):