python 遥感影像像元坐标读取到excel文件

# -*- coding: utf-8 -*-
# @Time : 2020/10/22 17:16
# @Author : ZhangMin
import os
import pandas as pd
from osgeo import gdal, osr

os.environ['PROJ_LIB'] = r'e:\Anaconda3\envs\cloneTF21\Library\share\proj'


class GRID:

    # 读影像文件
    def read_img(self, filename):
        dataset = gdal.Open(filename)  # 打开文件
        im_width = dataset.RasterXSize  # 栅格矩阵的列数
        im_height = dataset.RasterYSize  # 栅格矩阵的行数
        im_geotrans = dataset.GetGeoTransform()  # 仿射矩阵
        im_proj = dataset.GetProjection()  # 地图投影信息
        # 近红外波段
        im_data = dataset.ReadAsArray(0, 0, im_width, im_height)

        # 将转换为经纬度坐标
        # get CRS from dataset
        crs = osr.SpatialReference()
        crs.ImportFromWkt(dataset.GetProjectionRef())
        # create lat/long crs with WGS84 datum
        crsGeo = osr.SpatialReference()
        crsGeo.ImportFromEPSG(4326)  # 4326 is the EPSG id of lat/long crs   EPSG 权威编码   OpenGIS WKT坐标系统定义
        # get CRS from dataset
        t = osr.CoordinateTransformation(crs, crsGeo)

        del dataset
        return im_data, im_width, im_height, im_geotrans, im_proj, t


if __name__ == '__main__':
    run = GRID()
    filename = r"L:\2019-2020changjinagkou3\6SWIR_TV\TV.tif"

    im_data, im_width, im_height, im_geotrans, im_proj, t = run.read_img(filename)
    # (392020.0, 10.0, 0.0, 3499560.0, 0.0, -10.0)
    # im_geotrans[1]是像元宽度(影像在宽度上的分辨率)；
    # im_geotrans[5]是像元高度(影像在高度上的分辨率)
    # 这六个参数包括  左上角坐标，像元X、Y方向大小，旋转等信息。 要注意，Y方向的像元大小为负值。
    print(im_geotrans)
    # 左上角坐标
    print(im_geotrans[0])
    print(im_geotrans[3])
    print(im_proj)

    arrSlope = []
    lat_long_coor = [] # 用于存储每个像素的（X，Y）坐标
    for i in range(im_height):
        row = []
        for j in range(im_width):
            '''
            有了下面的公式， 可以把栅格的每个点代入公式求出其在地球上的实际位置 
            （当然，求出的数值的单位是在坐标系统中定义的那个，而非经纬度坐标， 要求经纬度坐标，还要通过一步转换
            '''

            px = im_geotrans[0] + i * im_geotrans[1] + j * im_geotrans[2]
            py = im_geotrans[3] + i * im_geotrans[4] + j * im_geotrans[5]

            # 转换为经纬度坐标
            (long, lat, z) = t.TransformPoint(px, py)
            # print(long, lat)
            lon_lat = [long,lat]
            row.append(lon_lat)
        

        lat_long_coor.append(row)

    df = pd.DataFrame(lat_long_coor)
    df.to_excel("lat_long_coor.xlsx",index=False,header=False)

    print(len(lat_long_coor))