【PythonRS】基于Python对栅格数据进行归一化（统一量纲至0~1）

有段时间没有更新Python处理栅格、矢量数据了，一部分是因为之前基本上已经把如何使用Python处理地理数据的方法覆盖完了，另一部分是因为最近有其他方面的知识需要学习和巩固。也是赶巧，最近有个项目需要构建模型对影像进行反演需要用到归一化，所以就编了一段代码，今天就和大家分享一下如何使用Python对栅格数据进行归一化的操作。

1 原理

归一化是一种简化计算的方式，即将有量纲的表达式经过变换化为无量纲的表达式，成为标量。在多种计算中都经常用到这种方法。归一化是一种无量纲处理手段，使物理系统数值的绝对值变成某种相对值关系。

我这里使用的是最常规的归一化算法，即最大-最小归一化，公式如下：

Y = (X - Xmin) / (Xmax - Xmin)

2 Python代码

既然原理搞懂了，接下来就是代码。代码部分也很简单，和之前计算NDVI差不多。

# -*- coding: utf-8 -*-
"""
@Time ： 2023/12/26 11:52
@Auth ： RS迷途小书童
@File ：Normalization of Raster Data.py
@IDE ：PyCharm
@Purpose：栅格数据归一化
@Web：博客地址:https://blog.csdn.net/m0_56729804
"""
import numpy as np
from osgeo import gdaldef Get_data(filepath):print("------------------------------------影像基础信息------------------------------------")ds = gdal.Open(filepath)  # 打开数据集datasetds_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 Normalization_value(filepath, out_path):print("-------------------------------------归一化处理-------------------------------------")ds = gdal.Open(filepath)  # 打开数据集datasetds_width = ds.RasterXSize  # 获取数据宽度ds_height = ds.RasterYSize  # 获取数据高度ds_geo = ds.GetGeoTransform()  # 获取仿射地理变换参数ds_prj = ds.GetProjection()  # 获取投影信息driver = gdal.GetDriverByName('GTiff')  # 载入数据驱动，用于存储内存中的数组ds_result = driver.Create(out_path, ds_width, ds_height, bands=ds.RasterCount, eType=gdal.GDT_Float64)print("正在遍历所有波段......")for i in range(1, ds.RasterCount):  # 遍历所有波段array_band = ds.GetRasterBand(i).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)min_value = np.min(array_band)  # 获取最小值max_value = np.max(array_band)  # 获取最大值array_result = (array_band-min_value)/(max_value-min_value)ds_result.GetRasterBand(i).WriteArray(array_result)  # 将结果写入数组array_result = Nonedel array_resultds_result.SetGeoTransform(ds_geo)  # 导入仿射地理变换参数ds_result.SetProjection(ds_prj)  # 导入投影信息ds_result.GetRasterBand(1).SetNoDataValue(0)  # 将无效值设为0del ds_result# 删除内存中的结果，否则结果不会写入图像中print("归一化处理已完成......")if __name__ == "__main__":print("\n")# 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"彭俊喜/MS600_0009_Multi.tif"  # 输入的栅格数据路径out_path1 = r"彭俊喜\test2.tif"  # 导出的文件路径Get_data(file_path)  # 执行函数，获取影像基本信息Normalization_value(file_path, out_path1)print("--------------------------------------程序结束--------------------------------------")