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

最新推荐文章于 2024-05-19 15:24:27 发布

sleepy-R

最新推荐文章于 2024-05-19 15:24:27 发布

阅读量2k

点赞数 3

CC 4.0 BY-SA版权

文章标签： python arcgis 开发语言

本文链接：https://blog.youkuaiyun.com/qq_40960899/article/details/129947548

本文介绍如何利用Python和GDAL库，在WGS UTM投影坐标系下计算图像四个角的投影坐标，并转换为WGS1984格式的经纬度坐标。通过GetGeoTransform方法获取地理信息，然后进行坐标转换。内容包括代码示例和验证过程。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

本文参考了文末的三个博客并进行修改和总结，在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