Python多值提取至点：基于一景栅格的像元提取另一栅格中空间位置匹配的像元_python 根据栅格点找出20km内的另一个栅格格点位置

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

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

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

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

  这里需要注意，如果待处理的2景遥感影像一个为地理坐标系，一个为投影坐标系，那么首先需要将2景遥感影像都处理为同一种类型的坐标系（建议都处理为投影坐标系）；具体处理方法，大家可以参考GDAL一行代码实现投影：将栅格的地理坐标系转为投影坐标系（https://blog.csdn.net/zhebushibiaoshifu/article/details/136200172）这篇文章。

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

# -*- coding: utf-8 -*-
"""
Created on Sun Feb 18 21:15:08 2024

@author: fkxxgis
"""

import numpy as np
from osgeo import gdal

def raster2array(file_path):
    dataset = gdal.Open(file_path)
    band = dataset.GetRasterBand(1)
    array = band.ReadAsArray()
    dataset = None
    return array

def get_geotransform(file_path):
    dataset = gdal.Open(file_path)
    geotransform = dataset.GetGeoTransform()
    dataset = None
    return geotransform

def get_pixel_size(geotransform):
    pixel_size_x = geotransform[1]
    pixel_size_y = geotransform[5]
    return pixel_size_x, pixel_size_y

def pixel2coordinate(geotransform, pixel_x, pixel_y):
    coordinate_x = geotransform[0] + pixel_x * geotransform[1] + pixel_y * geotransform[2]
    coordinate_y = geotransform[3] + pixel_x * geotransform[4] + pixel_y * geotransform[5]
    return coordinate_x, coordinate_y

gf_file_path = r"F:\Data_Reflectance_Rec\GF1WFV1.16m.2021001035028.48STA.000000_SR.tiff"
gf_array = raster2array(gf_file_path)
gf_geotransform = get_geotransform(gf_file_path)
gf_pixel_size_x, gf_pixel_size_y = get_pixel_size(gf_geotransform)

type_file_path = r"F:\Data_Reflectance_Rec\Type\result.tif"
type_array = raster2array(type_file_path)
type_geotransform = get_geotransform(type_file_path)
type_pixel_size_x, type_pixel_size_y = get_pixel_size(type_geotransform)

type_new_array = np.empty_like(gf_array)

for row in range(gf_array.shape[0]):
    for col in range(gf_array.shape[1]):
        if not gf_array[row, col]:
            type_new_array[row, col] = 0
            continue;
        gf_coordinate_x, gf_coordinate_y = pixel2coordinate(gf_geotransform, col, row)
        type_col = int((gf_coordinate_x - type_geotransform[0]) / type_pixel_size_x)
        type_row = int((gf_coordinate_y - type_geotransform[3]) / type_pixel_size_y)
        type_new_array[row, col] = type_array[type_row, type_col]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

  其中，我们首先需要引入必要的库。在本文中，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个二维数组所示的一致。

  至此，大功告成。

欢迎关注：疯狂学习GIS