通过GDAL实现tif图片坐标信息转换为WGS84

参考链接：

http://t.csdnimg.cn/uRIe4

一、实现地图坐标的转换

1.打开tif文件，通过gdal的GetGeoTransform读取其图像信息，其中：

# geoTransform[0]：左上角像素纬度
# geoTransform[1]：影像宽度方向上的分辨率(经度范围/像素个数)
# geoTransform[2]：x像素旋转, 0表示上面为北方
# geoTransform[3]：左上角像素经度
# geoTransform[4]：y像素旋转, 0表示上面为北方
# geoTransform[5]：影像宽度方向上的分辨率(纬度范围/像素个数)

2.通过如上参数计算投影坐标

3.其中图像中源数据为UTM格式，GetProjection获取tif图像中的投影坐标系信息

pcs = osr.SpatialReference()
pcs.ImportFromEPSG(4326)

通过如上代码来确定要转换到的坐标系信息，其中4326代表WGS1984坐标系；

通过TransformPoint函数将得到的经纬度坐标其转为GPS。

定义一个坐标系类：

class Coordinate():
    def __init__(self, extent, im_width, im_height):
        # 栅格数据的六参数。
        # geoTransform[0]：左上角像素纬度
        # geoTransform[1]：影像宽度方向上的分辨率(经度范围/像素个数)
        # geoTransform[2]：x像素旋转, 0表示上面为北方
        # geoTransform[3]：左上角像素经度
        # geoTransform[4]：y像素旋转, 0表示上面为北方
        # geoTransform[5]：影像宽度方向上的分辨率(纬度范围/像素个数)
        self.extent = extent
        self.lat_left = self.extent[0]
        self.x_resolution = self.extent[1]
        self.x_rotation = self.extent[2]
        self.lon_left = self.extent[3]
        self.y_resolution = self.extent[5]
        self.y_rotation = self.extent[4]
        self.im_width = im_width
        self.im_height = im_height
        self.lat_right = self.lat_left + self.x_resolution * self.im_width + self.im_height * self.x_rotation
        self.lon_right = self.lon_left + self.y_resolution * self.im_height + self.im_width * self.y_rotation
        print(self.x_resolution * self.im_width)
        # print(self.lat_left + self.x_resolution * self.im_width + self.im_height * self.x_rotation, self.y_resolution)
        # print(self.lat_right, self.lon_right)
    def coordinate_Lower_Right(self):
        return self.lat_right, self.lon_right

    def coordinate_Upper_Left(self):
        return self.lat_left, self.lon_left

通过如下代码进行转换： 

from osgeo import gdal, osr
from pyproj import Transformer
import os
import sys
from coordinate import Coordinate


filename = "D:/odm/project0/odm_orthophoto/odm_orthophoto.original.tif"
def GetTifInfo():
    dataset = gdal.Open(filename)  # 打开文件
    im_width = dataset.RasterXSize  # 栅格矩阵的列数
    im_height = dataset.RasterYSize  # 栅格矩阵的行数
    print(im_width, im_height)
    im_bands = dataset.RasterCount  # 波段数
    dem = dataset.GetRasterBand(1).ReadAsArray()
    # print(im_bands)
    # for i in dem:
    #     for j in i:
    #         if(j > 0):
    #             print(j)
    extent = dataset.GetGeoTransform()  # 仿射矩阵，左上角像素的大地坐标和像素分辨率
    print(extent)
    # var = dataset.GetProjection()  # 栅格数据的投影
    gps = Coordinate(extent, im_width, im_height)
    lat, lon = gps.coordinate_Upper_Left()
    latr, lonr = gps.coordinate_Lower_Right()
    print(latr, lonr)
    # osr.SpatialReference 提供描绘和转换坐标系统的服务 地理坐标(用经纬度表示)；投影坐标(如 UTM ，用米等单位量度来定位)。
    pcs = osr.SpatialReference()
    pcs.ImportFromEPSG(4326)
    # print(dataset.GetSpatialRef())
    # pcs.SetWellKnownGeogCS('WGS84')
    transform = osr.CoordinateTransformation(dataset.GetSpatialRef(), pcs)
    lat, lon, _ = transform.TransformPoint(lat, lon, 0.0)
    latr, lonr, _ = transform.TransformPoint(latr, lonr, 0.0)
    lat1, lon1,_ = transform.TransformPoint(409819.8955701462 - extent[1] * 683, 3208389.9585364168 - extent[5] * 455, 0.0)
    print(lat, lon)
    print(latr, lonr)
    print(lat1, lon1)
    # pcs.SetWellKnownGeogCS('EPSG:4326')
    # ImportFromWkt()函数可以把 WKT坐标系统设置到OGRSpatialReference中
    pcs.ImportFromWkt(dataset.GetProjection())

    gcs = pcs.CloneGeogCS()  # 地理空间坐标系
    shape = (dataset.RasterYSize, dataset.RasterXSize)
    img = dataset.GetRasterBand(1).ReadAsArray()  # (height, width)
    # print(img)
    # img(ndarray), gdal数据集、地理空间坐标系、投影坐标系、栅格影像大小
    return img, dataset, gcs, pcs, extent, shape


os.environ['PROJ_LIB'] = r'D:\Anaconda\envs\odm\Library\share\proj'
os.environ['GDAL_DATA'] = r'D:\Anaconda\envs\odm\Library\share'
GetTifInfo()