基于GDAL批量提取经纬度/投影坐标对应像元的值

我们知道栅格数据就是由大量紧凑排列的格子组成，其中，每一个格子都由一个固定的数值填充。用专业术语来说，“格子”被称为像元，因此要获取某个像元（格子）的值，只需要知道像元的位置（行列号）即可。“根据经纬度坐标、投影坐标来提取栅格图像的像元值”是学习研究中经常遇到的一个需求，其原理就是：将经纬度坐标或投影坐标转换为图像坐标（可理解为行列号），然后根据行列号索引来获取像元的值。本文，小编将讲解一下，如何基于GDAL，来编写Python脚本程序提取像元值。

任务说明

本文并没有为了解决某个实际任务来编写代码，但代码确实具有一定的实用意义。通过本文，您将收获两方面的技能知识（要注意的是，所谓“批量”，一般而言就是重复执行同一个功能，然后输出多个值，具体到代码里，就是在for循环或者while循环里执行同一段代码多次）：

掌握如何编写Python脚本来实现经纬度坐标、投影坐标和图像坐标之间的批量转换；（文中的代码已经编写好相应的函数）

掌握如何根据经纬度坐标或者投影坐标来批量地提取对应的像元值。

实践操作

完成上述任务的代码如下图所示：

代码注释：

上述代码并没有涉及到复杂的Python语法，因此这里不再做代码的注释，只简单说明一下几个GDAL函数的含义：

（1）osr.SpatialReference()建立了一个空间参考，dataset.GetProjection()函数是获取GDAL打开的数据集的投影信息返回值是一个WKT字符串，ImportFromWkt()则可以根据WKT字符串投影信息来创建投影坐标系相关的信息（包括地理坐标系、基准面、投影方法、分辨率等）

（2）pcs.CloneGeogCS()可以根据投影信息来创建地理坐标系，要知道先有了地理坐标系，而后才能创建投影坐标系，因此，根据投影坐标系很容易获得地理坐标系相关信息。

（3）dataset.GetGeoTransform()函数则是获取图像的地理空间范围和分辨率信息，返回值是包括了图像的左上角第一个像元中心的坐标（第1和第4个值）、X和Y方向的空间分辨率信息（第2和第6个值）的元组。其中，第3和第5个值分别对应这旋转系数和平移系数，一般如果遥感影像所在的条带严格指向地球的正北，则这两个系数均为0。如本文中所使用的影像的空间范围信息如下图所示：

结果展示

代码执行结果如下图所示，可以看出，无论是使用经纬度坐标，还是投影坐标，亦或是图像坐标，所提取的像元值均是相同的。当然，代码中使用的投影坐标和图像坐标均并不是通过其他软件计算得到，而是根据小编编写的代码依据经纬度坐标转换而来的，因此这一份代码的运行结果还可以证明，根据本文中的方法，大家还可以批量的将经纬度坐标转换为投影坐标。

为了检验程序的正确性，这里可以对比一下ArcMap下对应点的实际值（不是拉伸值！！！），如下图，小编在ArcMap中使用“Identify”工具，随机地找了6个点，并获取其像元值。与上面代码运行的结果图对比，可以看出，程序所提取的6个点的像元值完全一致，因此程序是无误的。在下一节文章中，小编将介绍如何使用ArcGIS的ArcPy包来提取像元值。至于哪种方法好，那就取决于大家的偏好楼，反正小编喜欢这篇文章中的方法。