有段时间没有更新Python处理栅格、矢量数据了,一部分是因为之前基本上已经把如何使用Python处理地理数据的方法覆盖完了,另一部分是因为最近有其他方面的知识需要学习和巩固。也是赶巧,最近有个项目需要构建模型对影像进行反演需要用到归一化,所以就编了一段代码,今天就和大家分享一下如何使用Python对栅格数据进行归一化的操作。
1 原理
归一化是一种简化计算的方式,即将有量纲的表达式经过变换化为无量纲的表达式,成为标量。在多种计算中都经常用到这种方法。归一化是一种无量纲处理手段,使物理系统数值的绝对值变成某种相对值关系。
我这里使用的是最常规的归一化算法,即最大-最小归一化,将数值归一化至0~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 gdal
def Get_data(filepath):
print("------------------------------------影像基础信息------------------------------------")
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 Normalization_value(filepath, out_path):
print("-------------------------------------归一化处理-------------------------------------")
ds = gdal.Open(filepath) # 打开数据集dataset
ds_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 = None
del array_result
ds_result.SetGeoTransform(ds_geo) # 导入仿射地理变换参数
ds_result.SetProjection(ds_prj) # 导入投影信息
ds_result.GetRasterBand(1).SetNoDataValue(0) # 将无效值设为0
del 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("--------------------------------------程序结束--------------------------------------")
标签:RS,Python,result,归一化,print,array,ds From: https://www.cnblogs.com/RSran/p/17930791.html如果大家在学习Python或者RS时有什么问题,可以随时留言交流!如果大家对批量处理有兴趣同样可以留言给博主,博主会分享相关代码以供学习!