首页
学习
活动
专区
工具
TVP
发布
精选内容/技术社群/优惠产品,尽在小程序
立即前往

基于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包来提取像元值。至于哪种方法好,那就取决于大家的偏好楼,反正小编喜欢这篇文章中的方法。

  • 发表于:
  • 原文链接https://kuaibao.qq.com/s/20180716G1FPCB00?refer=cp_1026
  • 腾讯「腾讯云开发者社区」是腾讯内容开放平台帐号(企鹅号)传播渠道之一,根据《腾讯内容开放平台服务协议》转载发布内容。
  • 如有侵权,请联系 cloudcommunity@tencent.com 删除。

扫码

添加站长 进交流群

领取专属 10元无门槛券

私享最新 技术干货

扫码加入开发者社群
领券