GDAL利用地理坐标读取图像像元值

最近的一个项目需要在电子海图中下载已知水深点，导出点的地理坐标（经纬度）。然后在arcgis中打开这些地理坐标输出为shp，利用GDAL读取不同波段的点对应的像元值，从而构建水深和像元值的对应关系。

其实想要根据经纬度得到像元值有两种方法。

1.通过arcpy选择点输出像元值，这个方法简单快捷，依赖于arcgis。

2.通过GDAL读取经纬度（度分秒），转换成小数形式。使用GDAL的CoordinateTransformation方法将经纬度转换成投影坐标，再转换成行列坐标，最后读取像元值。这里给出官方文档里面的代码以供参考。

补充一点：这里之所以把经纬度先转成投影坐标是因为遥感影像的坐标系是投影坐标系，而仿射六参的转换方法要求待转换点坐标和影像的投影坐标系一致，所以需要先转换成投影坐标。

　　　　　如果遥感影像的坐标系是地理坐标系（即没有投影坐标系），那么这里直接对经纬度坐标进行仿射六参转换即可得到行列号。

　　　　　而按下面的代码此时srs为none，srs与geocs的转换器仍然是对应的，那么经纬度转换后仍然是对应的格式，经过仿射变换最终得到行列号，也同样成立（该步骤可省略）。

 

其中GetGeoTransform获取获取图像的地理空间范围、分辨率信息。第1、4个值是左上角第一个像元中心的坐标。2、6是xy方向的空间分辨率。3、5是旋转系数和平移系数。

SpatialReference()建立空间参考，GetProjection返回投影信息，以wkt格式输出，最后用ImportFromWkt方法读取字符串并创建地理坐标系、基准面、投影方法、分辨率等。

CloneGeoCS方法创建地理坐标（经纬度），CoordinateTransformation方法建立投影坐标系和地理坐标系的转换关系，TranformPoint方法将地理坐标带入转换关系式，得到投影坐标，最后再计算得到行列数。

Tips

1.list out of range 应该是输入的列表格式有问题。

2.此外注意经度和纬度的输入位置。