GDAL对缺失投影定义的AIG文件根据经纬度坐标提取像元值

// dataset.GetProjection()

// '''
本脚本通过来拾取影像上的像素值，支持gdal可读的所有格式，支持读取方式：
1. input（文件+自设坐标信息) 仅当文件格式特殊且坐标系统没有EPSG编号时
2. input(文件) 基本情况下通用
'''

import numpy as np
from osgeo import gdal
from osgeo import osr
from tqdm import tqdm

def get_file_info(in_file_path, in_prj_config=None):
    '''
    v.1 根据指定的图像文件路径，以只读的方式打开图像。（仅支持Tif格式）
    v.2 读取原始的AIG—Arc/Info二进制网格，由于投影文件读取错误会导致坐标转换失败，
    事先获取坐标系统定义语句，用于保留投影信息
    v.3 预处理得到全国WGS84坐标系统tif影像，不再需要定义投影语句

    :param in_file_path:AIG—Arc/Info二进制网格路径
    :param in_prj_config:自设投影定义，示例：'PROJCS["Krasovsky_1940_Albers",GEOGCS["Unknown datum based upon the Krassowsky 1940 ellipsoid",DATUM["Not_specified_based_on_Krassowsky_1940_ellipsoid",SPHEROID["Krassowsky 1940",6378245,298.3,AUTHORITY["EPSG","7024"]],AUTHORITY["EPSG","6024"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",105],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",47],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'
    :return：影像数据集、空间坐标系、投影坐标系、行列像元数（像元宽度）
    '''
    #
    pcs = None
    gcs = None
    shape = None
    #

    gdal.AllRegister()
    dataset = gdal.Open(in_file_path)
    pcs = osr.SpatialReference()
    if in_prj_config != None:
        pcs.ImportFromWkt(in_prj_config)
    else:
        pcs.ImportFromWkt(dataset.GetProjection())
        #
    gcs = pcs.CloneGeogCS()
        #
    extend = dataset.GetGeoTransform()
        #
    shape = (dataset.RasterXSize, dataset.RasterYSize)

    return dataset, gcs, pcs, extend, shape

def lonlat_to_xy(gcs, pcs, lon, lat):
    '''
    经纬度坐标转换为投影坐标
    :param gcs:地理空间坐标信息，可由get_file_info()函数获取
    :param pcs:投影坐标信息，可由get_file_info()函数获取
    :param lon:经度坐标
    :param lat:纬度坐标
    :return:地理空间坐标对应的投影坐标
    '''
    #
    ct = osr.CoordinateTransformation(gcs, pcs)
    coordinates = ct.TransformPoint(lon, lat)
    return coordinates[0], coordinates[1], coordinates[2]

def xy_to_lonlat(gcs, pcs, x, y):
    '''
    投影坐标转换为经纬度坐标
    :param gcs:地理空间坐标信息，可由get_file_info()函数获取
    :param pcs:投影坐标信息，可由get_file_info()函数获取
    :param x:像元的行号
    :param y:像元的列号
    :return:投影坐标对应的地理空间坐标
    '''
    #
    ct = osr.CoordinateTransformation(pcs, gcs)
    lon, lat, _ = ct.TransformPoint(x, y)
    #
    return lon, lat

def xy_to_rowcol(extend, x, y):
    '''
    根据GDAL的六参数模型将给定的投影或地理坐标转为影像图上坐标（行列号）

    :param extend:图像的空间范围
    :param x:投影坐标x
    :param y:投影坐标y
    :return:投影坐标(x,y)对应的影像图像行列号(row,col)
    '''

    a = np.array([[extend[1], extend[2]], [extend[4], extend[5]]])
    b = np.array([x - extend[0], y - extend[3]])
    row_col = np.linalg.solve(a, b)
    try:
        row = int(np.floor(row_col[1]))
        col = int(np.floor(row_col[0]))
        return row, col

    except:
        return -1, -1

def rowcol_to_xy(extend, row, col):
    '''
    图像坐标转换为投影坐标

    根据GDAL的六参数模型将给定的影像图上坐标（行列号）转为投影或地理坐标（根据具体数据的坐标系统转换）
    :param extend:图像的空间范围
    :param row:像元的行号
    :param col:像元的列号
    :return:影像图像行列号(row,col)对应的投影坐标(x,y)
    '''
    #
    x = extend[0] + row * extend[1] + col * extend[2]
    y = extend[3] + row * extend[4] + col + extend[5]
    #
    return x, y

#根据单个坐标对获取
def get_single_value_by_coordinates(file_path, coordinates, prj_config=None):
   '''
   根据单个图像坐标，或者依据GDAL的六参数模型将给定的投影、地理坐标转为影像图上坐标后，返回对应像元的像素值
    :param file_path: 图像的文件路径
    :param coordinates: 坐标、一维列表，【地理空间坐标】，分别为经度、纬度
    :param prj_iconfig: 自设投影定义
    :return: 列表形式，单个坐标的像素值
   '''
   dataset, gcs, pcs, extend, shape = get_file_info(file_path, prj_config)
   img = dataset.GetRasterBand(1).ReadAsArray()

   x, y, _ = lonlat_to_xy(gcs, pcs, coordinates[0], coordinates[1])
   row, col = xy_to_rowcol(extend, x, y)
   try:
       return img[row, col]
   except:
       return 0

#根据多个坐标点对获取
def get_multi_values_by_coordinates(file_path, coordinates, prj_config=None):
    '''
    根据多个图像坐标，或者依据GDAL的六参数模型将给定的投影、地理坐标转为影像图上坐标后，返回对应像元的像素值
    :param file_path: 图像的文件路径
    :param coordinates: 坐标、二维列表，第二维为【地理空间坐标】
    :param prj_config: 自设投影定义
    :return: 列表形式，多个坐标的像素值
    '''
    dataset, gcs, pcs, extend, shape = get_file_info(file_path, prj_config)
    img = dataset.GetRasterBand(1).ReadAsArray()
    values = []
    for coord in tqdm(coordinates):
        x, y, _ = lonlat_to_xy(gcs, pcs, coord[0], coord[1])
        row, col = xy_to_rowcol(extend, x, y)
        try:
            values.append(img[row, col])
        except:
            values.append(0)
    return(values)

# if __name__ == "__main__":
#     coordinates_tuple = [[116.4, 39.9], [116.63, 40.31]]
#
#     #case1 普通格式
#     file_path1 = 'all_province.tif'
#     pixel_value1 = get_multi_values_by_coordinates(file_path1, coordinates_tuple)
#     print(pixel_value1)
#
#     #case2 特殊格式，自设投影
#     prj_config = 'PROJCS["Krasovsky_1940_Albers",GEOGCS["Unknown datum based upon the Krassowsky 1940 ellipsoid",DATUM["Not_specified_based_on_Krassowsky_1940_ellipsoid",SPHEROID["Krassowsky 1940",6378245,298.3,AUTHORITY["EPSG","7024"]],AUTHORITY["EPSG","6024"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",105],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",47],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'
#     file_path2 = 'beijing/beijing/hdr.adf'
#     pixel_value2 = get_multi_values_by_coordinates(file_path2, coordinates_tuple, prj_config)
#     print(pixel_value2)