输入行列号，获取遥感影像的对应位置的数值

在以前我在学习监督分类的时候，使用过这个功能。即输入行列号，获取遥感影像的对应位置的数值。

往往在写和遥感相关的代码时，我优先考虑使用python的gdal库去完成，而不是rasterio。这个是是个人的习惯。

rasterio很优秀，它是基于gdal库进行编写、扩张的python库，在安装完rasterio后，我们可以直接去学它的源码，它的源码结构很清晰，比较容易去理解，不过这是另一件事情了，我们还是把目光转向今天的主题：通过行列号获取遥感影像的具体位置的数值

现在使用单波段的遥感影像作为示例，讲一讲如何实现这个功能。

只需要十行代码即可。对于可以扩展的功能，一般会写成一个类而不是一个函数。

from osgeo import gdal
class GetImageXY:
    def __init__(self, image_path):
        self.image_path = image_path
        self.read_image()
    def read_image(self):
        """
        read image and get whole image array
        """
        dataset = gdal.Open(self.image_path)
        self.data = dataset.GetRasterBand(1).ReadAsArray()
    def get_xy_data(self, x, y):
        """
        get x and y from image array
        """
        return self.data[x, y]

以上代码，首先导入了gdal；然后建立了一个GetImageXY类，_init_是 Python 类的一个初始化方法（构造函数），image_path是影像文件的路径。

self.image_path = image_path

是把影像文件路径 赋予给这个类的属性，即在GetImageXY类中的其他方法函数，可以通过self.image_path 获取我们输入的参数，这样就避免重复输入参数。下图是GetImageXY类的属性。

一个是矩阵，一个是字符串，代表着影像的矩阵数据和影像的路径

image-20240918141009742

在点开data，可以看到data的详细信息，如下。

image-20240918141246447

然后建立一个read_image的类方法，用于获取整个影像的第一个波段的数据，并把这个数据赋予给GetImageXY类的属性，这里写为self.data，从此，在这类中self.data就代表了影像的第一个波段的数据。

同时把read_image的类方法写在初始化方法中，这样一旦实例化这个类，自动调用read_image类函数。

最后建立一个get_xy_data的类方法，输入x,y，返回影像中行列号对应的数值。

第一次测试

if __name__ == '__main__':
    image_path = r'LC08_L1TP_118038_20240912_20240912_02_RT_B3.TIF'
    A = GetImageXY(image_path)
    a = A.get_xy_data(2000, 2000)
    print('\n'*10)
    print(a)

    # 在控制台显示10803

说明了在行号为2000、列号为2000的位置上的影像的数值为10803

image-20240918140055163

我们在envi打开该影像，使用cursor value工具，获取2000,2000的数值，结果如下:

image-20240918133527040

验证成功。

获取多波段影像的xy位置的数值

获取多波段影像的xy位置的数值，原理和上面的一致，但需要改两行代码。

from osgeo import gdal
class GetImageXY:
    def __init__(self, image_path):
        self.image_path = image_path
        self.read_image()
    def read_image(self):
        """
        read image and get whole image array
        """
        dataset = gdal.Open(self.image_path)
        self.data = dataset.ReadAsArray()
    def get_xy_data(self, x, y):
        """
        get x and y from image array
        """
        if len(self.data.shape) <3:

            return self.data[x, y]
        else:
            return self.data[:, x, y]

在read_image函数中，改为不指定获取第几个波段的数据，直接使用dataset.ReadAsArray()

在get_xy_data函数中，对data的尺寸做判断，若data的尺寸为（x,y），尺寸小于3，则data是单波段影像

若data的尺寸为（z,x,y），尺寸大于等于3，则data不是单波段影像

所以，整体代码的运行逻辑如下:

image-20240918144517393

第二次测试

这次，我们玩个花的。使用的是普通图片来测试。众所周知，遥感影像只是图片的子集。

所以，这次测试，使用的是由微信截图获取的图片进行测试。

微信截图_20240918140700

改一下程序的输入参数

if __name__ == '__main__':
    image_path = r'LC08_L1TP_118038_20240912_20240912_02_RT_B3.TIF'
    image_path = r'微信截图.png'
    A = GetImageXY(image_path)
    a = A.get_xy_data(200, 200)
    print('\n'*20)
    print(a)

结果如下：image-20240918142432949

还是在envi打开该图片，使用cursor value工具，获取200,200的数值，结果如下:

image-20240918145429351

再放大，结果如下：

image-20240918142758732

第二次测试，结果也一致。

题外话

如果不输入单个行列号，怎么获取对应的位置。

这里分情况：

1.你有一个csv表格，存储了很多xy行列号，这时你需要做的是，使用类似于下面的代码，将行列号信息读取出来

with open(file) as f:
	text = f.readlines()

然后进行循环获取对应的影像数据的数值。

2.你有一个矢量文件，存储了很多点信息，这时你需要做的是，将你矢量点的经纬度提取出来，然后批量将经纬度换算为影像的行列号，最后进行循环获取对应的影像数据的数值。

流程图如下：

image-20240918144651883

《google earth engine：像素坐标和地理坐标的互相转换的技术谈论》