python实现栅格数据缺失值/nodata值/空白值填充

使用python实现栅格数据中缺失值的填充

weixin_52325473

1104人浏览 · 2024-10-22 20:42:16

weixin_52325473 · 2024-10-22 20:42:16 发布

声明：创作内容仅供个人记录使用，如需转载，请注明出处。

1.问题背景

在下载dem时，难免会遇到dem有缺失值的情况：

像是上面这种情况，我们可以通过ArcGIS中的【识别】工具查看这些位置的像素值，一般是nodata值

2.解决思路

我们要做的就是让这些位置获得一个相对正确的像素值，这里我们想起来，地理现象、地理要素是具有空间连续性的，因此我们可以通过这个像素点周围的像素值来推断这点的像素值。

这里我采用的是领域分析求均值的方法，举个例子：

这是含有空缺值的栅格数据打开后的样子，我这里把nodata替换为nan了。

我们要做的是，以nan值为中心，定义一个窗口，比如5*5，使用这个窗口中有效数据的平均值作为这个像元处的像素值。如图：

图中蓝色的区域就是5*5的大小，有效数据为：

valid_values=[1087. 1085. 1084. 1083. 1081. 1087. 1086. 1084. 1084. 1083. 1088. 1087., 1090. 1087. 1086. 1085. 1082. 1090. 1089. 1088. 1086. 1084.]

因此我们计算mean(valid_values)=1085.7273

然后重复上述操作。

注意：在栅格数据边界处也会有nodata值，这些值可能是由于裁剪等出现，这些边缘的数据我们不需要处理，因此在遍历栅格数组时，我们会忽略边界区域的nodata。

如果你想要处理边缘处的nodata，可以尝试ArcGIS中的外推

3.代码实现

import time
import numpy as np
from osgeo import gdal


def write_dem(dem, array, out_path):
    """
    这个函数用于保存填充结果，
    :param dem: 原始dem，用于获取投影等信息
    :param array: 填充后的数组
    :param out_path: 输出路径
    :return: 
    """
    print("开始写入数据")
    dataset = gdal.Open(dem)
    img_width = dataset.RasterXSize
    img_height = dataset.RasterYSize
    adf_GeoTransform = dataset.GetGeoTransform()
    im_Proj = dataset.GetProjection()
    band = dataset.GetRasterBand(1)
    nodata = band.GetNoDataValue()
    dtype = band.DataType

    driver = gdal.GetDriverByName("GTiff")
    datasetnew = driver.Create(out_path, img_width, img_height, 1, dtype)
    datasetnew.SetGeoTransform(adf_GeoTransform)
    datasetnew.SetProjection(im_Proj)
    band = datasetnew.GetRasterBand(1)
    band.SetNoDataValue(nodata)

    band.WriteArray(array)
    datasetnew.FlushCache()
    print("写入数据完成")


def fill_nodata(dempath, outpath, windowsize):
    """
    填充主函数，使用时只需要调用这个函数
    :param dempath: 含有缺失值dem的路径
    :param outpath: 输出结果路径
    :param windowsize: 领域分析的参数，必须是奇数，一般在5-17比较合适，想知道最适值可以使用均值变点法
    :return: 
    """
    print("开始填充栅格数据")
    start = time.perf_counter()
    dataset = gdal.Open(dempath)
    band = dataset.GetRasterBand(1)
    before_array = band.ReadAsArray().astype(np.float32)
    nodata = band.GetNoDataValue()
    before_array[before_array == nodata] = np.nan # 将栅格数组中的nodata值替换为nan

    after_array = np.copy(before_array) # 先将before_array复制为after_array

    if windowsize % 2 == 0:
        print("输入的窗口大小应该为奇数")
        return

    half_window = windowsize // 2

    for i in range(half_window, before_array.shape[0] - half_window):
        for j in range(half_window, before_array.shape[1] - half_window): # 对于栅格数据边界处半个windowsize的区域不遍历
            if np.isnan(before_array[i, j]):
                window = before_array[i - half_window:i + half_window + 1, j - half_window:j + half_window + 1]

                valid_values = window[~np.isnan(window)] # 只需要窗口中的有效值
                if len(valid_values) > 0:
                    after_array[i, j] = np.nanmean(valid_values) # 对after_array中的nan值替换
    after_array = np.nan_to_num(after_array, nan=nodata)
    write_dem(dempath, after_array, outpath)
    end = time.perf_counter()
    print("填充栅格数据完成，计算时间为 %.3f 秒\n" % (end - start))


if __name__ == '__main__':
    # 使用方法示例：
    dempath=r"E:\myPaper\Data\fill\填洞测试.tif"
    out_path=r"E:\myPaper\Data\30m\fill_result.tif"
    fill_nodata(dempath, out_path, 5)

使用说明：

如果你想要使用这段代码，全部复制之后，修改：

dempath为你的栅格数据地址，out_path为你要保存的位置

然后点击运行即可