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, lon, lat):
    '''
    將經緯度座標轉為投影座標（具體的投影座標系由給定資料確定）
    :param dataset: GDAL地理資料
    :param lon: 地理座標lon經度
    :param lat: 地理座標lat緯度
    :return: 經緯度座標(lon, lat)對應的投影座標
    '''

    prosrs, geosrs = getSRSPair(dataset)
    ct 
 
= osr.CoordinateTransformation(geosrs, prosrs)
    coords = ct.TransformPoint(lon, lat)
    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進行二元一次方程的求解