【Python&RS】基于GDAL修改栅格数据的DN值

        遥感工作者离不开栅格数据，有时候我们可能需要修改栅格数据的值，但ENVI和ArcGIS中并没有直接修改DN值的工具，只有栅格计算器、Band math这些工具去计算整个波段的值，或者Edit Classification Image工具可以修改ENVI分类后的像元值，但这个工具只对分类格式有效，博主整不明白怎么把单波段数据转为分类格式，所以这些工具都无法满足我们的需求。

        今天跟大家分享一下如何使用Python的GDAL库将栅格数据中特定的DN值修改成我们想要的。

一、安装库

import os
import numpy as np
from osgeo import gdal

二、读取栅格数据信息

        这一步可有可无，只是让你了解一下栅格数据的基本信息，如投影信息、长度、宽度等。

def Get_data(filepath):
    ds = gdal.Open(filepath)  # 打开数据集dataset
    ds_width = ds.RasterXSize  # 获取数据宽度
    ds_height = ds.RasterYSize  # 获取数据高度
    ds_bands = ds.RasterCount  # 获取波段数
    ds_geo = ds.GetGeoTransform()  # 获取仿射地理变换参数
    ds_prj = ds.GetProjection()  # 获取投影信息
    print("影像的宽度为：" + str(ds_width))
    print("影像的高度为：" + str(ds_height))
    print("仿射地理变换参数为：" + str(ds_geo))
    print("投影坐标系为：" + str(ds_prj))
    # data = ds.ReadAsArray(0, 0, ds_width, ds_height)  # 以数组的形式读取整个数据集

三、修改DN值

        我这里写的函数，只适用于修改单波段的栅格影像，如果需要修改多波段就自己加一个循环。其实原理都一样，就是将波段读成数组，然后通过索引读取第几行第几列像元的值，然后判断这个值是否为你想要修改的值，如果是，就将其赋予新的值。

def Modify_value(filepath, out_path, original, target):
    print("-----正在进行DN值的修改-----")
    ds = gdal.Open(filepath)  # 打开数据集dataset
    ds_width = ds.RasterXSize  # 获取数据宽度
    ds_height = ds.RasterYSize  # 获取数据高度
    ds_geo = ds.GetGeoTransform()  # 获取仿射地理变换参数
    ds_prj = ds.GetProjection()  # 获取投影信息
    array_band = ds.GetRasterBand(1).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)
    # 读取第一个波段全部
    for row in range(0, ds_height):
        # 循环当前波段的行
        for col in range(0, ds_width):
            # 循环当前波段的列
            if array_band[row][col] == original:
                # 判断第row行第col列的DN值是否为需要修改的值
                array_band[row][col] = target
                # 修改该值
            else:
                continue
    driver = gdal.GetDriverByName('GTiff')  # 载入数据驱动，用于存储内存中的数组
    ds_result = driver.Create(out_path, ds_width, ds_height, bands=1, eType=gdal.GDT_Float64)
    # 创建一个数组，宽高为原始尺寸
    ds_result.SetGeoTransform(ds_geo)  # 导入仿射地理变换参数
    ds_result.SetProjection(ds_prj)  # 导入投影信息
    ds_result.GetRasterBand(1).SetNoDataValue(0)  # 将无效值设为0
    ds_result.GetRasterBand(1).WriteArray(array_band)  # 将结果写入数组
    del ds_result
    # 删除内存中的结果，否则结果不会写入图像中
    print("计算完成")

四、完整代码

# -*- coding: utf-8 -*-
"""
@Time ： 2023/9/8 8:51
@Auth ： RS迷途小书童
@File ：Modifying Raster Data DN Values.py
@IDE ：PyCharm
@Purpose：修改栅格数据DN值
"""
import os
import numpy as np
from osgeo import gdal
 
 
def Get_data(filepath):
    ds = gdal.Open(filepath)  # 打开数据集dataset
    ds_width = ds.RasterXSize  # 获取数据宽度
    ds_height = ds.RasterYSize  # 获取数据高度
    ds_bands = ds.RasterCount  # 获取波段数
    ds_geo = ds.GetGeoTransform()  # 获取仿射地理变换参数
    ds_prj = ds.GetProjection()  # 获取投影信息
    print("影像的宽度为：" + str(ds_width))
    print("影像的高度为：" + str(ds_height))
    print("仿射地理变换参数为：" + str(ds_geo))
    print("投影坐标系为：" + str(ds_prj))
    # data = ds.ReadAsArray(0, 0, ds_width, ds_height)  # 以数组的形式读取整个数据集
 
 
def Modify_value(filepath, out_path, original, target):
    print("-----正在进行DN值的修改-----")
    ds = gdal.Open(filepath)  # 打开数据集dataset
    ds_width = ds.RasterXSize  # 获取数据宽度
    ds_height = ds.RasterYSize  # 获取数据高度
    ds_geo = ds.GetGeoTransform()  # 获取仿射地理变换参数
    ds_prj = ds.GetProjection()  # 获取投影信息
    array_band = ds.GetRasterBand(1).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)
    # 读取第一个波段全部
    for row in range(0, ds_height):
        # 循环当前波段的行
        for col in range(0, ds_width):
            # 循环当前波段的列
            if array_band[row][col] == original:
                # 判断第row行第col列的DN值是否为需要修改的值
                array_band[row][col] = target
                # 修改该值
            else:
                continue
    driver = gdal.GetDriverByName('GTiff')  # 载入数据驱动，用于存储内存中的数组
    ds_result = driver.Create(out_path, ds_width, ds_height, bands=1, eType=gdal.GDT_Float64)
    # 创建一个数组，宽高为原始尺寸
    ds_result.SetGeoTransform(ds_geo)  # 导入仿射地理变换参数
    ds_result.SetProjection(ds_prj)  # 导入投影信息
    ds_result.GetRasterBand(1).SetNoDataValue(0)  # 将无效值设为0
    ds_result.GetRasterBand(1).WriteArray(array_band)  # 将结果写入数组
    del ds_result
    # 删除内存中的结果，否则结果不会写入图像中
    print("计算完成")
 
 
if __name__ == "__main__":
    os.environ['PROJ_LIB'] = 'G:/Anaconda/envs/pyDL/Lib/site-packages/osgeo/data/proj'
    os.environ['GDAL_DATA'] = 'G:/Anaconda/envs/pyDL/Lib/site-packages/osgeo/data'
    # 添加PROJ至环境变量，消除警告
    file_path = r"B:\1m_xiugai.tif"  # 输入的栅格数据路径
    out_path1 = r"B:\1m_xiugai1.tif"  # 导出的文件路径
    data1 = int(input("请输入需要修改的DN值:"))
    data2 = int(input("请输入目标DN值     :"))
    Get_data(file_path)  # 执行函数，获取影像基本信息
    print("\n")
    print("--------------DN值修改--------------")
    Modify_value(file_path, out_path1, data1, data2)  # 执行函数，修改DN值

        今天的分享就到这里了，大家需要注意的是，我这段代码只适用于单波段数据且想要修改的值只有一种时，如你想要将所有DN值等于1的像元全部改成0，就可以直接使用我的点吗改数据路径，然后再输入1和0就可以了（因为我的任务就是将分类数据（DN值为0，1，2，3，4）中分错的部分改成正确的）。不要问为什么不使用工具，因为我手里的分类数据不是ENVI支持的分类格式（泪目）。

        如果大家在学习Python或者RS时有什么问题，可以随时留言交流！如果大家对批量处理有兴趣同样可以留言给博主，博主会分享相关代码以供学习！