任务背景：需要根据经纬度坐标提取AIG文件（AIG—Arc/Info二进制网格）对应像素值
了解到gdal能够完成这项任务，但是之前没有接触过gdal，所以现在网络上查找资料，发现如下链接所示的教程。
基于GDAL批量提取经纬度/投影坐标对应像元的值
查找gdal支持的数据格式，了解gdal支持AIG数据格式：
gdal文档
具体格式介绍如上，只需知在给予‘hdr.adf'文件的路径的条件下即可打开AIG文件
直接在上述教程进行测试
发现能够顺利读取AIG，但是根据正确坐标返回的坐标为像素值为空（或者在行列计算时就不存在），思考该问题应该是投影系统出现了问题。
打开QGIS对AIG文件进行检查
坐标系统unamed
发现我的AIG文件的坐标系统无法识别，也就是说明没有EPSG编号，但是该文件在QGIS中能够正常加载。
AIG文件投影unamed查看prj文件
打开'prj.adf'，虽然获取了投影信息，但是不知道怎样得到投影定义的表达式。
投影信息获取投影表达的方式
在QGIS中将原本的AIG文件转为tiff格式文件，打开tiff文件源信息：
通过格式转换获取到真实投影点击右侧的投影信息：
image.png可以看到左下角的投影定义语句，感兴趣的同学试一试直接使用左下角WKT信息是否能够成功。我是通过gdal读取tiff文件，然后使用下面代码获取的。
代码语言：txt复制// dataset.GetProjection()获取的投影信息也有了，接下来是对源代码进行个人定制，需要在原始函数上增加一项输入投影信息的参数。
代码实现
代码语言：txt复制// '''
本脚本通过来拾取影像上的像素值，支持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)
