Skip to content

使用 QGIS、GDAL 和 Python 修复栅格数据缺失值

来源: https://spatialthoughts.com/2020/06/17/fix-nodata-values/

处理栅格数据时，有时需要处理数据缺失（data gaps）。这可能是传感器故障、处理错误或数据损坏的结果。下图是航空影像中数据缺失（即无数据值）的示例。

来源图像：© Commission for Lands (COLA) ; Revolutionary Government of Zanzibar (RGoZ)，下载自 OpenAerialMap。（注意：数据缺失是使用 Python 脚本模拟的，并非原始数据集的一部分）

如果数据缺失较小，可以通过对相邻像素的值进行插值来有效解决。我将概述两种修复方法。第一种使用 QGIS，另一种使用纯 Python。

这里展示的方法使用 gdal_fillnodata 工具应用反距离加权插值和平滑。正如 文档 中指出的，这适用于填充连续栅格数据（如高程）中的缺失区域。它也适用于变化数据（如航空影像）中的非常小的空隙。如果你想对点数据进行插值以创建栅格，则应改用 gdal_grid 工具。在 QGIS 中，这可以通过 处理 → 工具箱 → 栅格分析 → 网格 (Grid) … 找到几个单独的算法。

在 QGIS 中修复数据缺失

GDAL 自带一个 gdal_fillnodata 工具，可以在 QGIS 的处理工具箱中使用。

如果源栅格已设置 NoData 值，且该值与缺失数据值相同，则可以跳过此步骤。否则，第一步是将栅格的 NoData 值设置为数据缺失处的像素值。在 处理 → 工具箱 中，搜索并定位 转换（通用格式转换）(Translate (convert format)) 工具。

在我们的示例中，NoData 像素值为 0。在 为输出波段指定 NoData 值 (Assign a specified nodata value to output bands) 选项中设置值 0，并输入转换后栅格的文件名。

现在我们可以运行处理工具箱中的 填充 NoData (Fill nodata) 工具。

此工具一次处理 1 个波段。选择 波段 1 (Red)。将 搜索插值值的最大距离 (Maximum distance to search out for values to interpolate) 设置为 1，因为我们只有 1 个像素的间隙。将输出保存为 01_red.tif 并点击运行。使用 01_ 这样的前缀保存文件很重要，因为下一步将按文件名的字母顺序合并波段。

对波段 2 (Green) 和波段 3 (Blue) 重复此过程，为它们选择合适的文件名。你应该有 3 个已填充无数据值的独立栅格。现在我们可以将它们合并为一个文件。在处理工具箱中搜索并定位 合并 (Merge) 工具。

在合并工具中，选择所有 3 个单独的栅格。勾选 将每个输入文件放入单独的波段 (Place each input file into a separate band) 复选框。输入输出文件名并点击运行。

生成的合并栅格将有 3 个波段，无数据空隙将被相邻像素的插值填充。

这里是展示修复前后版本的动画。

使用 Python 修复数据缺失

通常，当你的数据集中出现此类问题时，它会在所有源图像中持续存在。如果你有大量数据存在此类数据缺失，在 QGIS 中手动修复是不可行的。所以这里有另一种使用 Python 的方法。

上过我的 Python 空间分析基础 (Python Foundation for Spatial Analysis) 课程的同学可能知道，正是这个特定问题促使我学习 Python，因为当时没有现成的解决方案。下面是一个脚本，展示了如何借助 rasterio 和 numpy 库在 Python 中解决这个问题。

# 填充缺失行的数据
# 此脚本演示如何读取某些行存在缺失数据（即 0）的图像，并用相邻行的平均值填充。
import rasterio
import numpy as np

filename = 'original.tif'
dataset = rasterio.open(filename)
metadata = dataset.meta

red = dataset.read(1)
green = dataset.read(2)
blue = dataset.read(3)

# 下面的代码返回一个数组，其中每一项根据条件 `red==0` 为 True/False
result = np.all(red == 0, axis=1)

def average_rows(array, index):
    # 计算上一行和下一行的平均值，并四舍五入
    result = np.round(np.mean( np.array([array[index[0]-1], array[index[0]+1] ]), axis=0 ))
    array[index] = result

for index, x in np.ndenumerate(result):
    # 确保不是第一行也不是最后一行
    if (x and index[0] != 0 and index[0] != (dataset.height - 1)):
        average_rows(red, index)
        average_rows(blue, index)
        average_rows(green, index)

output_filename = 'fixed.tif'
metadata.update({'driver': 'GTiff'})

# 写入修复后的数据
rgb_dataset = rasterio.open(output_filename, 'w', **metadata)
rgb_dataset.write(red, 1)
rgb_dataset.write(green, 2)
rgb_dataset.write(blue, 3)
rgb_dataset.close()

如果你想尝试这个例子，可以 下载演示数据集，其中包含航空影像和带有 Python 代码的 Jupyter Notebook。