Python+GDAL 计算图像四个角的经纬度

本文参考了文末的三个博客并进行修改和总结，在ArcGIS中验证了结果

本文使用的图像是基于WGS UTM投影坐标系，计算出的结果为投影坐标（M），需要将结果转为WGS1984坐标系才可得到经纬度坐标

首先输出图像的投影信息，查看图像的投影格式

# 输出图像的投影信息
print(ds.GetProjectionRef())

输出结果如下所示，是UTM投影格式，如果是WGS1984地理坐标系算出来应该直接是经纬度坐标，不用转换

PROJCS["WGS 84 / UTM zone 44N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",81],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","32644"]]

1.首先计算图像四个角的投影坐标，代码如下

GetGeoTransform方法可以读取地理信息，返回六个参数

GT(0) 左上像素左上角的x坐标。
GT(1) w-e像素分辨率/像素宽度。
GT(2) 行旋转（通常为零）。
GT(3) 左上像素左上角的y坐标。
GT(4) 列旋转（通常为零）。
GT(5) n-s像素分辨率/像素高度（北上图像为负值）

from osgeo import gdal

def project_xy(tif_path):
    dataset = gdal.Open(tif_path)
    geo_information = dataset.GetGeoTransform()
    col = dataset.RasterXSize  # 行数
    row = dataset.RasterYSize  # 列数
    # band = dataset.RasterCount  # 波段

    # 左上角经纬度方向投影坐标
    top_left_corner_lon = geo_information[0]
    top_left_corner_lat = geo_information[3]

    # 左下角经纬度方向投影坐标
    bottom_left_corner_lon = geo_information[0] + row * geo_information[2]
    bottom_left_corner_lat = geo_information[3] + row * geo_information[5]

    # 右上角经纬度方向投影坐标
    top_right_corner_lon = geo_information[0] + col * geo_information[1]
    top_right_corner_lat = geo_information[3] + col * geo_information[4]

    # 右下角经纬度方向投影坐标
    bottom_right_corner_lon = geo_information[0] + col * geo_information[1] + row * geo_information[2]
    bottom_right_corner_lat = geo_information[3] + col * geo_information[4] + row * geo_information[5]

    return top_left_corner_lon, top_left_corner_lat, bottom_left_corner_lon, bottom_left_corner_lat, top_right_corner_lon, top_right_corner_lat, bottom_right_corner_lon, bottom_right_corner_lat

if __name__ == '__main__':
    all_tuple = project_xy('tif路径')
    # 左上角经纬度方向投影坐标(经度方向投影坐标，纬度方向投影坐标)
    top_left = [all_tuple[0],all_tuple[1]]
    print('左上角经纬度方向投影坐标',top_left)
    # 左下角经纬度方向投影坐标(经度方向投影坐标，纬度方向投影坐标)
    top_right = [all_tuple[2], all_tuple[3]]
    print('左下角经纬度方向投影坐标', top_right)
    # 右上角经纬度方向投影坐标(经度方向投影坐标，纬度方向投影坐标)
    bottom_left = [all_tuple[4], all_tuple[5]]
    print('右上角经纬度方向投影坐标', bottom_left)
    # 右下角经纬度方向投影坐标(经度方向投影坐标，纬度方向投影坐标)
    bottom_right = [all_tuple[6], all_tuple[7]]
    print('右下角经纬度方向投影坐标', bottom_right)

2.接下来是将计算出的经纬度方向投影坐标转换为WGS1984格式的经纬度坐标，代码如下（因为只需要用到将投影坐标转换为经纬度坐标的方法，只验证了代码中的geo2lonlat和lonlat2geo方法，其它方法没有验证）

from osgeo import gdal
from osgeo import osr
import numpy as np

def getSRSPair(dataset):
    '''
    获得给定数据的投影参考系和地理参考系
    :param dataset: GDAL地理数据
    :return: 投影参考系和地理参考系
    '''
    prosrs = osr.SpatialReference()
    prosrs.ImportFromWkt(dataset.GetProjection())
    geosrs = prosrs.CloneGeogCS()
    return prosrs, geosrs

def geo2lonlat(dataset, x, y):
    '''
    将投影坐标转为经纬度坐标（具体的投影坐标系由给定数据确定）
    :param dataset: GDAL地理数据
    :param x: 投影坐标x
    :param y: 投影坐标y
    :return: 投影坐标(x, y)对应的经纬度坐标(lon, lat)
    '''
    prosrs, geosrs = getSRSPair(dataset)
    ct = osr.CoordinateTransformation(prosrs, geosrs)
    coords = ct.TransformPoint(x, y)
    return coords[:2]

def lonlat2geo(dataset, lat, lon):
    '''
    将经纬度坐标转为投影坐标（具体的投影坐标系由给定数据确定）
    :param dataset: GDAL地理数据
    :param lon: 地理坐标lon经度
    :param lat: 地理坐标lat纬度
    :return: 经纬度坐标(lat, lon)对应的投影坐标
    '''
    prosrs, geosrs = getSRSPair(dataset)
    ct = osr.CoordinateTransformation(geosrs, prosrs)
    # 该方法需要先输入纬度后输入经度，不然就出错啦
    coords = ct.TransformPoint(lat, lon)
    return coords[:2]

def imagexy2geo(dataset, row, col):
    '''
    根据GDAL的六参数模型将影像图上坐标（行列号）转为投影坐标或地理坐标（根据具体数据的坐标系统转换）
    :param dataset: GDAL地理数据
    :param row: 像素的行号
    :param col: 像素的列号
    :return: 行列号(row, col)对应的投影坐标或地理坐标(x, y)
    '''
    trans = dataset.GetGeoTransform()
    px = trans[0] + col * trans[1] + row * trans[2]
    py = trans[3] + col * trans[4] + row * trans[5]
    return px, py

def geo2imagexy(dataset, x, y):
    '''
    根据GDAL的六 参数模型将给定的投影或地理坐标转为影像图上坐标（行列号）
    :param dataset: GDAL地理数据
    :param x: 投影或地理坐标x
    :param y: 投影或地理坐标y
    :return: 影坐标或地理坐标(x, y)对应的影像图上行列号(row, col)
    '''
    trans = dataset.GetGeoTransform()
    a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])
    b = np.array([x - trans[0], y - trans[3]])
    return np.linalg.solve(a, b)  # 使用numpy的linalg.solve进行二元一次方程的求解

if __name__ == '__main__':
    gdal.AllRegister()
    dataset = gdal.Open(r"1.tif")
    print('数据投影：')
    print(dataset.GetProjection())
    print('数据的大小（行，列）：')
    print('(%s %s)' % (dataset.RasterYSize, dataset.RasterXSize))
	
	# 此处的坐标自行设定，即输入需要转换的坐标值
    # 经度方向投影坐标
    x = 539980.0
    # 纬度方向投影坐标
    y = 2340000.0
    # 纬度
    lat = 21.667277636420355
    # 经度
    lon = 80.99980668837317

    print('投影坐标 -> 经纬度：')
    coords = geo2lonlat(dataset, x, y)
    # coords[0]纬度, coords[1]经度
    print('(%s, %s)->(%s, %s)' % (x, y, coords[0], coords[1]))
    print('经纬度 -> 投影坐标：')
    coords = lonlat2geo(dataset, lat, lon)
    # 经度方向投影坐标coords[0], 纬度方向投影坐标coords[1]
    print('(%s, %s)->(%s, %s)' % (lon, lat, coords[0], coords[1]))

3.最后完整版代码如下：

if __name__ == '__main__':
	fileroot = 'TIFF文件路径'
    # 加载数据
    dataset = gdal.Open(fileroot)
    # 输出图像的投影信息(WGS UTM投影)
    # print(dataset.GetProjectionRef())
    # 先计算左上和右下角的投影坐标
    all_tuple = project_xy(fileroot)
    # 左上角投影坐标
    # 经度方向
    top_left_x = all_tuple[0]
    # 纬度方向
    top_left_y = all_tuple[1]
    # 右下角投影坐标
    # 经度方向
    bottom_right_x = all_tuple[6]
    # 纬度方向
    bottom_right_y = all_tuple[7]

    # 将左上角投影坐标换算为经纬度
    top_left_lon_lat = geo2lonlat(dataset, top_left_x, top_left_y)
    # 左上角经度
    top_left_lon = top_left_lon_lat[1]
    # 左上角纬度
    top_left_lat = top_left_lon_lat[0]

    # 将右下角投影坐标换算为经纬度
    bottom_right_lon_lat = geo2lonlat(dataset, bottom_right_x, bottom_right_y)
    # 右下角经度
    bottom_right_lon = bottom_right_lon_lat[1]
    # 右下角纬度
    bottom_right_lat = bottom_right_lon_lat[0]

参考博文：
https://blog.csdn.net/weixin_43575792/article/details/123620337
https://blog.csdn.net/ChaoChao66666/article/details/124844722
https://blog.csdn.net/qq_32306361/article/details/128553141