获取遥感影像的真正的顶点
获取遥感影像的真正的顶点
ytkz遥感影像的真正的顶点
如上图所示，红色的四个点是遥感影像中的真正的顶点，指的是影像有效区域的边界点，而非无效值（如NoData或背景值）的顶点。
而蓝色的四个点，影像的地理范围（边界）的角点。这些顶点是影像的四个角点（左上、右上、右下、左下），表示影像的地理范围。与影像的有效区域无关，即使影像中有大量无效值（如NoData），GDAL的顶点仍然是影像的四个角点。
GDAL只保存左上角信息，其余三个角通过基于影像的地理坐标（经纬度或投影坐标）计算的。
代码如下：
from osgeo import gdal
# 打开影像
dataset = gdal.Open('image.tif')
width = dataset.RasterXSize
height = dataset.RasterYSize
geotransform = dataset.GetGeoTransform()
# 计算四个角点的坐标
corners = [
(geotransform[0], geotransform[3]),
(geotransform[0] + width * geotransform[1], geotransform[3]),
(geotransform[0] + width * geotransform[1], geotransform[3] + height * geotransform[5]),
(geotransform[0], geotransform[3] + height * geotransform[5])
]
print("GDAL顶点坐标：", corners)区别总结
| 特性 | 真正的顶点 | GDAL中的顶点 |
|---|---|---|
| 定义 | 有效区域的边界顶点 | 影像地理范围的四个角点 |
| 与无效值的关系 | 仅包含有效区域的顶点 | 与无效值无关，始终是四个角点 |
| 数量 | 取决于有效区域的形状，可能不止四个点 | 固定为四个点 |
| 用途 | 裁剪、几何校正、有效区域分析 | 地理范围计算、影像配准 |
| 获取方法 | 图像处理（如OpenCV） | GDAL库的地理变换信息 |
示例对比
- 真正的顶点：
- 如果影像中有一个不规则的有效区域（如一个湖泊），真正的顶点是这个湖泊边界的角点。
- 这些顶点是通过图像处理技术提取的。
- GDAL中的顶点：
- 无论影像中有效区域是什么形状，GDAL的顶点始终是影像的四个角点。
- 这些顶点是基于影像的地理坐标计算的。
真正顶点的获取
代码如下：
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time : 2025/2/18 11:32
# @File : get_vertex_point.py
'''
'''
import numpy as np
from osgeo import gdal
# DatasetInfo：遥感影像的属性，增加“获取源影像的非无效值的顶点坐标”
class DatasetInfo():
def __init__(self, path, type=True):
'''
@param path: 影像路径
@param num: 是否获取顶点坐标
'''
# 调用父类的构函
self.dataset = gdal.Open(path)
self.rows = self.dataset.RasterYSize # todo 图像宽度
self.cols = self.dataset.RasterXSize # todo 图像长度
self.bands = self.dataset.RasterCount # TODO 图像波段数量
self.proj = self.dataset.GetProjection() # todo 地图投影信息
self.geotrans = self.dataset.GetGeoTransform() # todo 仿射矩阵
self.longitude = self.geotrans[0] # todo 经度
self.latitude = self.geotrans[3] # todo 纬度
self.pixelWidth = self.geotrans[1] # todo x轴空间间隔
self.pixelHeight = self.geotrans[5]
if type==True:
try:
self.vertex_coordinates_info = self.vertex_coordinates()
except:
print('vertex_coordinates of DatasetInfo have error.')
def vertex_coordinates(self):
'''
获取源影像的非无效值的顶点坐标
输入影像的格式为 必定存在无效值的影像
@return:
'''
# ds = gdal.Open(path)
data = self.dataset.GetRasterBand(1).ReadAsArray()
nandata = data[0, 0]
result = []
# for i in range(self.rows):
conditions = 0
i = 0
while conditions == 0:
data1 = data[i, :]
if np.mean(data1) != nandata:
for j in range(self.cols):
if data1[j] != nandata:
print(data1[j])
result.append([i, j])
conditions = 1
break
i += 1
conditions = 0
j = 0
while conditions == 0:
data1 = data[:, j]
if np.mean(data1) != nandata:
for i in range(self.rows):
if data1[i] != nandata:
print(data1[i])
result.append([i, j])
conditions = 1
break
j += 1
conditions = 0
i = self.rows
while conditions == 0:
data1 = data[i - 1, :]
if np.mean(data1) != nandata:
for j in range(self.cols):
if data1[j] != nandata:
print(data1[j])
result.append([i, j])
conditions = 1
break
i -= 1
conditions = 0
j = self.cols
while conditions == 0:
data1 = data[:, j - 1]
if np.mean(data1) != nandata:
for i in range(self.rows - 1, 0, -1):
if data1[i] != nandata:
print(data1[i])
result.append([i, j])
conditions = 1
break
j -= 1
return result
if __name__ == "__main__":
file = r'D:\orth.tif'
result = DatasetInfo(file)
print(result.vertex_coordinates_info)
改进
上面的代码，得到的是影像有效范围的四个图像坐标点，一般来说，我们需要的影像的地理坐标点。所以我们再进一步把图像坐标转换为地理坐标。
整体代码如下：
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time : 2025/2/27 09:25
# @File : get_vertex_point.py
"""
获取遥感影像非无效值的顶点坐标
依赖库：numpy, gdal (osgeo)
"""
import numpy as np
from osgeo import gdal
class DatasetInfo:
"""遥感影像属性类，用于获取影像元数据及非无效值顶点坐标"""
def __init__(self, path, get_vertices=True):
"""
初始化影像数据集
Args:
path (str): 影像文件路径
get_vertices (bool): 是否计算顶点坐标，默认为True
"""
# 打开数据集并获取基本属性
self.dataset = gdal.Open(path)
self.rows = self.dataset.RasterYSize # 图像高度（行数）
self.cols = self.dataset.RasterXSize # 图像宽度（列数）
self.bands = self.dataset.RasterCount # 波段数
self.proj = self.dataset.GetProjection() # 投影信息
self.geotrans = self.dataset.GetGeoTransform() # 仿射变换参数
# 提取地理坐标信息
self.longitude = self.geotrans[0] # 左上角经度
self.latitude = self.geotrans[3] # 左上角纬度
self.pixelWidth = self.geotrans[1] # 像素宽度（x方向）
self.pixelHeight = self.geotrans[5] # 像素高度（y方向，通常为负值）
# 可选计算顶点坐标
if get_vertices:
try:
self.vertex_coordinates_info = self.vertex_coordinates()
self.vertex_geo_coordinates()
except Exception as e:
print(f'Error in vertex_coordinates calculation: {e}')
def vertex_coordinates(self):
"""
获取影像非无效值的顶点坐标（从四条边向内搜索第一个非无效值点）
Returns:
list: 包含四个顶点坐标的列表，每个坐标为[i, j]形式
"""
# 读取第一个波段的数据
data = self.dataset.GetRasterBand(1).ReadAsArray()
nodata_value = data[0, 0] # 假设左上角值为无效值
result = []
# 1. 从上向下扫描顶部边
i, found = 0, False
while not found and i < self.rows:
row_data = data[i, :]
if np.mean(row_data) != nodata_value:
for j in range(self.cols):
if row_data[j] != nodata_value:
print(f"Top edge value: {row_data[j]}")
result.append([i, j])
found = True
break
i += 1
# 2. 从左向右扫描左侧边
j, found = 0, False
while not found and j < self.cols:
col_data = data[:, j]
if np.mean(col_data) != nodata_value:
for i in range(self.rows):
if col_data[i] != nodata_value:
print(f"Left edge value: {col_data[i]}")
result.append([i, j])
found = True
break
j += 1
# 3. 从下向上扫描底部边
i, found = self.rows, False
while not found and i > 0:
row_data = data[i - 1, :]
if np.mean(row_data) != nodata_value:
for j in range(self.cols):
if row_data[j] != nodata_value:
print(f"Bottom edge value: {row_data[j]}")
result.append([i, j])
found = True
break
i -= 1
# 4. 从右向左扫描右侧边
j, found = self.cols, False
while not found and j > 0:
col_data = data[:, j - 1]
if np.mean(col_data) != nodata_value:
for i in range(self.rows - 1, -1, -1):
if col_data[i] != nodata_value:
print(f"Right edge value: {col_data[i]}")
result.append([i, j])
found = True
break
j -= 1
return result
def geo_coordinates(self, i, j):
"""
将图像坐标转换为地理坐标
Args:
i (int): 图像行索引
j (int): 图像列索引
Returns:
tuple: 包含地理坐标的元组 (longitude, latitude)
"""
longitude = self.longitude + j * self.pixelWidth
latitude = self.latitude + i * self.pixelHeight
return longitude, latitude
def vertex_geo_coordinates(self):
"""
获取影像非无效值的顶点地理坐标
Returns:
list: 包含四个顶点地理坐标的列表，每个坐标为[longitude, latitude]形式
"""
self.vertex_geo_coordinates_info = []
for i, j in self.vertex_coordinates_info:
self.vertex_geo_coordinates_info.append(self.geo_coordinates(i, j))
if __name__ == "__main__":
file_path = r'D:\orth.tif'
dataset = DatasetInfo(file_path)
print("Vertex coordinates:", dataset.vertex_coordinates_info)可视化
为了将四个地理坐标点可视化，可以用Python的matplotlib库绘制出来。以下是实现代码：
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
绘制由经纬度坐标定义的多边形
依赖库：matplotlib
"""
import matplotlib.pyplot as plt
# 定义四个顶点坐标
coordinates = [
(98.845352, 40.200816), # 点1
(98.787688, 40.005784000000006), # 点2
(99.051168, 39.959368000000005), # 点3
(99.10960800000001, 40.154368000000005) # 点4
]
# 解包经度和纬度
longitudes, latitudes = zip(*coordinates)
# 为了闭合多边形，将第一个点添加到末尾
longitudes = list(longitudes) + [longitudes[0]]
latitudes = list(latitudes) + [latitudes[0]]
# 创建图形
plt.figure(figsize=(8, 6))
plt.plot(longitudes, latitudes, 'b-', label='Polygon') # 蓝色线条绘制多边形
plt.scatter(longitudes[:-1], latitudes[:-1], c='red', label='Vertices') # 红色点标记顶点
plt.grid(True)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('Polygon Defined by Given Coordinates')
plt.legend()
# 显示图形
plt.show()输出结果
运行代码后，你会看到一个窗口显示由这四个点构成的多边形。横轴为经度（Longitude），纵轴为纬度（Latitude），顶点用红点标记，多边形边框用蓝色线条绘制。
