??本文介紹基于Python語言中的gdal模塊，對2景不同的遙感影像加以對應位置像素值匹配的方法——即基于一景遙感影像的每一個像元，提取另一景遙感影像中，與之空間位置相同的像元的像素值的方法。
??首先，明確一下本文的需求。現(xiàn)在有2景成像范圍不完全一致、但是具有重疊部分的遙感影像，如下圖所示；我們就將其稱作大遙感影像（成像范圍更大的、灰色系的那一景遙感影像）和小遙感影像（成像范圍更小的、藍色系的那一景遙感影像）。這2景遙感影像的成像范圍、空間分辨率、空間坐標系等都不一致（當然，也可以一致；而且如果一致的話，后續(xù)處理起來也會更方便理解一些）。
??其中，可以很明顯地看到，小遙感影像的空間分辨率高于大遙感影像，但其成像范圍是小于大遙感影像的；如下圖所示。
??我們現(xiàn)在希望，對于小遙感影像中的每一個像元（除了NoData值的像元），找到其在大遙感影像中對應位置處的像元，并將這個大遙感影像對應像元的像素提取出來?？梢哉J為，我們希望得到2個相同大小的二維數(shù)組——這2個二維數(shù)組的行數(shù)、列數(shù)就是小遙感影像的行數(shù)與列數(shù)，而這2個二維數(shù)組的值，分別為小遙感影像的像素值，以及大遙感影像在同一空間位置上的像元的像素值。換句話說，這個需求有點類似于ArcGIS的“多值提取到點”這一工具的作用——只不過相當于我們需要對小遙感影像的每一個像元都執(zhí)行一次“多值提取到點”操作。
??這里需要注意，如果待處理的2景遙感影像一個為地理坐標系，一個為投影坐標系，那么首先需要將2景遙感影像都處理為同一種類型的坐標系（建議都處理為投影坐標系）；具體處理方法，大家可以參考地理坐標系變換至投影坐標系：GDAL命令行這篇文章。
??明確了需求，我們即可開始代碼的撰寫。本文所用代碼如下。
# -*- 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\若爾蓋GF反射率\2021\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]
??其中，我們首先需要引入必要的庫。在本文中，numpy用于處理數(shù)組數(shù)據(jù)，gdal則用于讀取柵格數(shù)據(jù)文件和獲取地理轉換參數(shù)。
??隨后，我們定義了幾個關鍵的函數(shù)。其中，raster2array()用于將柵格數(shù)據(jù)文件讀取為numpy庫的數(shù)組，get_geotransform()用于獲取柵格數(shù)據(jù)文件的地理轉換參數(shù)，get_pixel_size()用于從地理轉換參數(shù)中提取像素大小，pixel2coordinate()用于將像素坐標轉換為地理坐標。
??接下來，我們即可開始讀取待處理的數(shù)據(jù)。在上述代碼中，gf_開頭的數(shù)據(jù)就是前文中提到的小遙感影像對應的相關數(shù)據(jù)，而type_開頭的數(shù)據(jù)就是前文中提到的大遙感影像對應的相關數(shù)據(jù)。
??首先，我們使用raster2array()函數(shù)將小遙感影像讀取為數(shù)組，并存儲在gf_array變量中；隨后，使用get_geotransform()函數(shù)獲取小遙感影像的地理轉換參數(shù)，并存儲在gf_geotransform變量中；接下來，使用get_pixel_size()函數(shù)從小遙感影像的地理轉換參數(shù)中提取像素大小，并分別存儲在gf_pixel_size_x和gf_pixel_size_y變量中。
??類似地，對大遙感影像文件同樣執(zhí)行上一段中描述的操作。
??接下來，創(chuàng)建一個與小遙感影像的數(shù)組具有相同形狀和數(shù)據(jù)類型的空數(shù)組。在這里，我們直接使用np.empty_like()函數(shù)創(chuàng)建名為type_new_array的空數(shù)組，其形狀與gf_array相同。
??隨后，遍歷小遙感影像的數(shù)組（相當于就是按行、按列遍歷小遙感影像的全部像元），根據(jù)條件進行處理。其中，如果gf_array中的元素為0（也就是我這里小遙感影像的NoData值），則不用再進行后續(xù)處理了，直接將type_new_array相應位置的元素也設置為0并繼續(xù)下一個循環(huán)。而如果gf_array中的元素不為0，根據(jù)像素坐標和地理轉換參數(shù)進行計算，從類型數(shù)據(jù)數(shù)組type_array中獲取相應位置的值，并將其賦值給type_new_array相應位置的元素。
??執(zhí)行上述代碼后，我們來檢查一下代碼的運行是否符合預期。因為大遙感影像的空間分辨率低一些，所以我們就用它來驗證我們的結果（空間分辨率低一些的話，驗證起來反而更方便）。首先，如下圖所示，可以看到，我們的代碼結果（也就是type_new_array這個數(shù)組）顯示，當行號為19與20時，其像素值發(fā)生了變化。
??我們到ArcGIS中驗證一下，將小遙感影像從左上角開始，向下數(shù)20行，可以看到對應的像元（如下圖中左下角的紫色框內所示）確實位于大遙感影像像元的分界處，且二者的像素值也都和上圖中2個二維數(shù)組所示的一致。
??至此，大功告成。
