㈠ 基於python編程的Modis地表溫度數據缺失值批量填補——以克里金插值法為例
本文介紹一種基於Python編程的柵格數據缺失值填補方法，利用克里金插值法實現。對於內存有限的筆記本電腦，我們採用逐張處理策略，用戶可根據自身需求進行優化以提高處理速度。以下是實施步驟與代碼示例：
首先，導入所需的庫：
python
import os
import numpy as np
from scipy.interpolate import griddata
from osgeo import gdal
接著，定義插值函數：
python
def interpolate_raster(input_raster, output_raster):
dataset = gdal.Open(input_raster)
band = dataset.GetRasterBand(1)
array = band.ReadAsArray()
nonzero_indices = np.where(array != 0)
points = np.array(list(zip(nonzero_indices[1], nonzero_indices[0])))
values = array[nonzero_indices]
x_range = np.arange(0, array.shape[1], 1)
y_range = np.arange(0, array.shape[0], 1)
grid_x, grid_y = np.meshgrid(x_range, y_range)
interpolated_values = griddata(points, values, (grid_x, grid_y), method='cubic')
driver = gdal.GetDriverByName('GTiff')
output_dataset = driver.Create(output_raster, array.shape[1], array.shape[0], 1, gdal.GDT_Float32)
output_dataset.SetGeoTransform(dataset.GetGeoTransform())
output_dataset.SetProjection(dataset.GetProjection())
output_band = output_dataset.GetRasterBand(1)
output_band.WriteArray(interpolated_values)
output_band.FlushCache()
print("Interpolation completed. Output raster saved as", output_raster)
接下來，定義處理流程：
python
input_folder = '輸入路徑'
output_folder = '輸出路徑'
if not os.path.exists(output_folder):
os.makedirs(output_folder)
for filename in os.listdir(input_folder):
if filename.endswith('.tif'):
input_raster = os.path.join(input_folder, filename)
output_raster = os.path.join(output_folder, filename)
interpolate_raster(input_raster, output_raster)
執行代碼後，盡管處理速度可能較慢，但實現了批量填補缺失值的目標。此方法對於內存有限的環境較為實用，通過調整代碼，用戶可以進一步優化處理效率。