duxingchen/WQ_GUI
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

fix: 实测站坐标改osr重投影+GDAL逆仿射,提取光谱写CSV前统一0~1量纲

Browse Source 
- geo_to_pixel 改用 gdal.InvGeoTransform 全仿射(含旋转项);坐标列正则收紧,只认经纬度,剔除 pixel/utm 误判
- get_spectral_in_coor 主流程由错误 UTM 数学公式改为 osr.CoordinateTransformation(EPSG源->影像原生CRS) 重投影
- 提取光谱列写 CSV 前做 >10 则 /10000 量级收敛(不动实测值列)
This commit is contained in:
duxin
2026-09-08 17:26:18 +08:00
parent c6523cebd8
commit 4d49753cca
1 changed files with 80 additions and 132 deletions
Download Patch File Download Diff File Expand all files Collapse all files

212
src/core/glint_removal/get_spectral.py
View File

@ -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'):