多值提取至点Python实现：找到不同栅格相同位置像元

  本文介绍基于Python语言中的gdal模块，对2景不同的遥感影像加以对应位置像素值匹配的方法——即基于一景遥感影像的每一个像元，提取另一景遥感影像中，与之空间位置相同的像元的像素值的方法。

  首先，明确一下本文的需求。现在有2景成像范围不完全一致、但是具有重叠部分的遥感影像，如下图所示；我们就将其称作大遥感影像（成像范围更大的、灰色系的那一景遥感影像）和小遥感影像（成像范围更小的、蓝色系的那一景遥感影像）。这2景遥感影像的成像范围、空间分辨率、空间坐标系等都不一致（当然，也可以一致；而且如果一致的话，后续处理起来也会更方便理解一些）。

  其中，可以很明显地看到，小遥感影像的空间分辨率高于大遥感影像，但其成像范围是小于大遥感影像的；如下图所示。

  我们现在希望，对于小遥感影像中的每一个像元（除了NoData值的像元），找到其在大遥感影像中对应位置处的像元，并将这个大遥感影像对应像元的像素提取出来。可以认为，我们希望得到2个相同大小的二维数组——这2个二维数组的行数、列数就是小遥感影像的行数与列数，而这2个二维数组的值，分别为小遥感影像的像素值，以及大遥感影像在同一空间位置上的像元的像素值。换句话说，这个需求有点类似于ArcGIS的“多值提取到点”这一工具的作用——只不过相当于我们需要对小遥感影像的每一个像元都执行一次“多值提取到点”操作。

  这里需要注意，如果待处理的2景遥感影像一个为地理坐标系，一个为投影坐标系，那么首先需要将2景遥感影像都处理为同一种类型的坐标系（建议都处理为投影坐标系）；具体处理方法，大家可以参考GDAL命令：一行代码转换坐标系这篇文章。

  明确了需求，我们即可开始代码的撰写。本文所用代码如下。

1# -*- coding: utf-8 -*-

2"""

3Created on Sun Feb 18 21:15:08 2024

4

5@author: fkxxgis

6"""

7

8import numpy as np

9from osgeo import gdal

10

11def raster2array(file_path):

12    dataset = gdal.Open(file_path)

13    band = dataset.GetRasterBand(1)

14    array = band.ReadAsArray()

15    dataset = None

16    return array

17

18def get_geotransform(file_path):

19    dataset = gdal.Open(file_path)

20    geotransform = dataset.GetGeoTransform()

21    dataset = None

22    return geotransform

23

24def get_pixel_size(geotransform):

25    pixel_size_x = geotransform[1]

26    pixel_size_y = geotransform[5]

27    return pixel_size_x, pixel_size_y

28

29def pixel2coordinate(geotransform, pixel_x, pixel_y):

30    coordinate_x = geotransform[0] + pixel_x * geotransform[1] + pixel_y * geotransform[2]

31    coordinate_y = geotransform[3] + pixel_x * geotransform[4] + pixel_y * geotransform[5]

32    return coordinate_x, coordinate_y

33

34gf_file_path = r"F:\Data_Reflectance_Rec\若尔盖GF反射率\2021\GF1WFV1.16m.2021001035028.48STA.000000_SR.tiff"

35gf_array = raster2array(gf_file_path)

36gf_geotransform = get_geotransform(gf_file_path)

37gf_pixel_size_x, gf_pixel_size_y = get_pixel_size(gf_geotransform)

38

39type_file_path = r"F:\Data_Reflectance_Rec\Type\result.tif"

40type_array = raster2array(type_file_path)

41type_geotransform = get_geotransform(type_file_path)

42type_pixel_size_x, type_pixel_size_y = get_pixel_size(type_geotransform)

43

44type_new_array = np.empty_like(gf_array)

45

46for row in range(gf_array.shape[0]):

47    for col in range(gf_array.shape[1]):

48        if not gf_array[row, col]:

49            type_new_array[row, col] = 0

50            continue;

51        gf_coordinate_x, gf_coordinate_y = pixel2coordinate(gf_geotransform, col, row)

52        type_col = int((gf_coordinate_x - type_geotransform[0]) / type_pixel_size_x)

53        type_row = int((gf_coordinate_y - type_geotransform[3]) / type_pixel_size_y)

54        type_new_array[row, col] = type_array[type_row, type_col]

  其中，我们首先需要引入必要的库。在本文中，numpy用于处理数组数据，gdal则用于读取栅格数据文件和获取地理转换参数。

  随后，我们定义了几个关键的函数。其中，raster2array()用于将栅格数据文件读取为numpy库的数组，get_geotransform()用于获取栅格数据文件的地理转换参数，get_pixel_size()用于从地理转换参数中提取像素大小，pixel2coordinate()用于将像素坐标转换为地理坐标。

  接下来，我们即可开始读取待处理的数据。在上述代码中，gf_开头的数据就是前文中提到的小遥感影像对应的相关数据，而type_开头的数据就是前文中提到的大遥感影像对应的相关数据。

  首先，我们使用raster2array()函数将小遥感影像读取为数组，并存储在gf_array变量中；随后，使用get_geotransform()函数获取小遥感影像的地理转换参数，并存储在gf_geotransform变量中；接下来，使用get_pixel_size()函数从小遥感影像的地理转换参数中提取像素大小，并分别存储在gf_pixel_size_x和gf_pixel_size_y变量中。

  类似地，对大遥感影像文件同样执行上一段中描述的操作。

  接下来，创建一个与小遥感影像的数组具有相同形状和数据类型的空数组。在这里，我们直接使用np.empty_like()函数创建名为type_new_array的空数组，其形状与gf_array相同。

  随后，遍历小遥感影像的数组（相当于就是按行、按列遍历小遥感影像的全部像元），根据条件进行处理。其中，如果gf_array中的元素为0（也就是我这里小遥感影像的NoData值），则不用再进行后续处理了，直接将type_new_array相应位置的元素也设置为0并继续下一个循环。而如果gf_array中的元素不为0，根据像素坐标和地理转换参数进行计算，从类型数据数组type_array中获取相应位置的值，并将其赋值给type_new_array相应位置的元素。

  执行上述代码后，我们来检查一下代码的运行是否符合预期。因为大遥感影像的空间分辨率低一些，所以我们就用它来验证我们的结果（空间分辨率低一些的话，验证起来反而更方便）。首先，如下图所示，可以看到，我们的代码结果（也就是type_new_array这个数组）显示，当行号为19与20时，其像素值发生了变化。

  我们到ArcGIS中验证一下，将小遥感影像从左上角开始，向下数20行，可以看到对应的像元（如下图中左下角的紫色框内所示）确实位于大遥感影像像元的分界处，且二者的像素值也都和上图中2个二维数组所示的一致。

  至此，大功告成。