【Python&RS】遥感影像的像素坐标转地理坐标（仿射变换）

GDAL(Geospatial Data Abstraction Library)是一个在X/MIT许可协议下的开源栅格空间数据转换库。它利用抽象数据模型来表达所支持的各种文件格式。它还有一系列命令行工具来进行数据转换和处理。

Python的GDAL库作为栅格数据的处理转换库，其支持几百种栅格数据格式，如常见的TIFF、ENVI、HFA、HDF4等。因为遥感影像大部分都是栅格数据，所以GDAL库非常适合处理遥感影像、如光谱指数计算、波段合成、批量下载、栅格转面等。

本次介绍如何通过遥感影像的仿射地理变换参数将像素坐标转为地理/投影坐标，在ENVI或者ArcGIS中都是自动转换，所以你在点击影像的某一像元时，它就会显示对应的地理/投影坐标。但如果你想用Python去实现，就需要使用影像的放射地理参数。

一、获取仿射地理变换参数

之前发布【Python&RS】GDAL计算遥感影像光谱指数（如NDVI、NDWI、EVI等）时就已经用过这部分代码，没啥好说的，就是获取影像的基本参数。代码中除了仿射地理变换参数其他没啥用，但我习惯全部输出。

def Get_data(filepath):

    """

    :param filepath: 输入影像的路径

    :return: 返回仿射地理变换参数

    """

    ds = gdal.Open(filepath)  # 打开数据集dataset

    ds_width = ds.RasterXSize  # 获取数据宽度

    ds_height = ds.RasterYSize  # 获取数据高度

    ds_bands = ds.RasterCount  # 获取波段数

    ds_geo = ds.GetGeoTransform()  # 获取仿射地理变换参数

    ds_prj = ds.GetProjection()  # 获取投影信息

    print("影像的宽度为：" + str(ds_width))

    print("影像的高度为：" + str(ds_height))

    print("仿射地理变换参数为：" + str(ds_geo))

    print("投影坐标系为：" + str(ds_prj))

    return ds_geo

    # data = ds.ReadAsArray(0, 0, ds_width, ds_height)  # 以数组的形式读取整个数据集

二、像素坐标转地理/投影坐标

这里先介绍一下仿射地理变换参数(3497294, 0.1, 0.0, 3438148, 0.0, -0.1)，这样大家可以更直观的明白代码的意思。

0：图像左上角的X坐标

1：图像东西方向分辨率

2：旋转角度，如果图像北方朝上，该值为0

3：图像左上角的Y坐标

4：旋转角度，如果图像北方朝上，该值为0

5：图像南北方向分辨率

def Pixel_Coordinate(x, y, ds_geo):

    """

    :param x: 输入x像素坐标

    :param y: 输入y像素坐标

    :param ds_geo: 输入仿射地理变换参数

    :return: 返回x,y的地理/投影坐标

    """

    x_geo = ds_geo[0]+ds_geo[1]*x + y*ds_geo[2]

    y_geo = ds_geo[3]+ds_geo[4]*x + y*ds_geo[5]

    print(x_geo, y_geo)

    return x_geo, y_geo

从代码中我们可以看到使用仿射地理变换参数将像素坐标转成地理坐标时，就是将左上角的x,y值加上x,y的像素坐标乘以分辨率（即实际距离），然后再加上偏转就行了。这和我之前写的通过图片中心点地理坐标推算整张图片所有像素点地理坐标的原理是一样的，只不过这里的仿射地理变换参数是影像自带的，而之前写的【Python&GIS】根据像素坐标计算图片某点的地理/投影坐标中的偏转角度是自己算出来的。

三、地理/投影坐标转像素坐标

def Coordinate_Pixel(x_geo, y_geo, ds_geo):

    """

    :param x_geo: 输入x地理坐标

    :param y_geo: 输入y地理坐标

    :param ds_geo: 输入仿射地理变换参数

    :return: 返回x,y像素坐标

    """

    y = ((y_geo - ds_geo[3] - ds_geo[4] / ds_geo[1] * x_geo + ds_geo[4] / ds_geo[1] * ds_geo[

        0]) / (ds_geo[5] - ds_geo[4] / ds_geo[1] * ds_geo[2]))

    x = ((x_geo - ds_geo[0] - y * ds_geo[2]) / ds_geo[1])

    return int(x), int(y)