|
Landsat遥感影像大气校正是逐波段进行的，因此数据存放方式分别为B1、B2、B3、B4、B5、B6、B7、B8，在应用前习惯上把各个波段合成及将各个波段的数值转换为地表真实的反射率数据（为了减少校正结果存储空间，程序中将大气校正的结果放大了10000倍），话不多说直接上代码：  图1 大气校正完后各波段# -*- coding: utf-8 -*-
"""
@autho:YuanRuiQiang
@date:2023/3/31 15:24
@Location:BeiJing_PieSat
@breif:实现遥感影像波段合成及转换地表真实反射率
"""
from osgeo import gdal, osr
import os
import glob
import numpy as np
def Read_img(filename):
# 数据读取
dataset = gdal.Open(filename)
width = dataset.RasterXSize
height = dataset.RasterYSize
band = dataset.RasterCount
im_data = dataset.ReadAsArray(0, 0, width, height) # 必要是输出影像数组，进行返回
geotrans = dataset.GetGeoTransform()
proj = dataset.GetProjection()
# nodata = dataset.GetRasterBand(1).GetNoDataValue()
return im_data, width, height, proj, geotrans
def write_tiff(filename, proj, geotrans, data):
# 数据写出
# gdal数据类型包括
# gdal.GDT_Byte,
# gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
# gdal.GDT_Float32, gdal.GDT_Float64
# 判断栅格数据的数据类型
if 'int8' in data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
# 判读数组维数
if len(data.shape) == 3:
bands, height, width = data.shape
else:
bands = 1
height, width = data.shape
# 创建文件
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(filename, width, height, bands, datatype)
dataset.SetGeoTransform(geotrans)
dataset.SetProjection(proj)
if bands == 1:
dataset.GetRasterBand(1).WriteArray(data)
else:
for i in range(bands):
dataset.GetRasterBand(i + 1).WriteArray(data[i])
del dataset
# ############################ 遥感影像波段合成 #########################
def LayerStacking(inputPath,outfile_dir):
"""
将多光谱波段合成一个tif
:param dstlist: 输入但波段影像列表
:param outfile_dir: 影像的输出文件夹
"""
Rasters = glob.glob(os.path.join(inputPath + '\\' + "*B[1-7].TIF")) # *B[1-8].TIF
raster = Rasters[1]
a1, bandname = os.path.split(raster)
(name, ext) = os.path.splitext(bandname)
bandname = name[0:-6] # 提取文件名
dataset_init = gdal.Open(raster)
# 创建待输出的图
gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create(outfile_dir + "\\" + bandname + '.tif', dataset_init.RasterXSize, dataset_init.RasterYSize, 7, gdal.GDT_Float32)
# 本文合成的对象是4个波段，按照自己的需要改变
out_ds.SetProjection(dataset_init.GetProjection())
out_ds.SetGeoTransform(dataset_init.GetGeoTransform()) # 获得原始波段的地理信息
# 往栅格数据集中填各波段数值
for i in range(len(Rasters)):
dataset = gdal.Open(Rasters[i])
print('影像' + Rasters[i] + '开始波段合成')
band_temp = dataset.GetRasterBand(1)
band_temp.SetNoDataValue(-9999)
out_ds.GetRasterBand(i+1).WriteArray(band_temp.ReadAsArray()) # 注意band从1开始，所以要+1
# dataset.GetRasterBand(i + 1).WriteArray(data[i])
del out_ds
def GetSurfaceReflectance (coastal, blue, green, red, nir, swir1, swir2):
"""
得到地表发射率数据
:param coastal: landsat8的B1波段,蓝色波段blue
:param blue: landsat8的B2波段,蓝色波段blue
:param green: landsat8的B3波段,绿色波段green
:param red: landsat8的B4波段,红外波段red
:param nir: landsat8的B5波段,近红外波段nir
:param swir1: landsat8的B6波段,短波红外1 swir1
:param swir2: landsat8的B7波段,短波红外2 swir2
:return: 计算结果
"""
coastal = coastal * 0.0001 # 大气校正后的数据将各波段的值放大了10000倍数
blue = blue * 0.0001
green = green * 0.0001
red = red * 0.0001
nir = nir * 0.0001
swir1 = swir1 * 0.0001
swir2 = swir2 * 0.0001
data = np.array([coastal, blue, green, red, nir, swir1, swir2])
return data
if __name__ == '__main__':
print('开始Landsat波段合成')
InputPath = r"D:\algorithm\GEO\RSEI_code\data\Python_6s_AtmosphericCorrection"
Out = r"D:\algorithm\GEO\RSEI_code\data\Python_6s_AtmosphericCorrection\layerStacking"
totalout = r"D:\algorithm\GEO\RSEI_code\data\Python_6s_AtmosphericCorrection\layerStacking\SurfaceReflectance.tif"
LayerStacking(InputPath, Out)
# 读取波段数据
print('开始转换真实影像反射率')
Rasters = glob.glob(os.path.join(Out + '\\' + "*.tif"))
raster = Rasters[0]
im_data, width, height, proj, geotrans = Read_img(raster)
coastal = im_data[0]
blue = im_data[1]
green = im_data[2]
red = im_data[3]
nir = im_data[4]
swir1 = im_data[5]
swir2 = im_data[6]
data = GetSurfaceReflectance(coastal, blue, green, red, nir, swir1, swir2)
write_tiff(totalout, proj, geotrans, data)
print("影像波段合成及真是地表反射率计算完成")
|
图1 大气校正完后各波段# -*- coding: utf-8 -*-
"""
@autho:YuanRuiQiang
@date:2023/3/31 15:24
@Location:BeiJing_PieSat
@breif:实现遥感影像波段合成及转换地表真实反射率
"""
from osgeo import gdal, osr
import os
import glob
import numpy as np
def Read_img(filename):
# 数据读取
dataset = gdal.Open(filename)
width = dataset.RasterXSize
height = dataset.RasterYSize
band = dataset.RasterCount
im_data = dataset.ReadAsArray(0, 0, width, height) # 必要是输出影像数组，进行返回
geotrans = dataset.GetGeoTransform()
proj = dataset.GetProjection()
# nodata = dataset.GetRasterBand(1).GetNoDataValue()
return im_data, width, height, proj, geotrans
def write_tiff(filename, proj, geotrans, data):
# 数据写出
# gdal数据类型包括
# gdal.GDT_Byte,
# gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
# gdal.GDT_Float32, gdal.GDT_Float64
# 判断栅格数据的数据类型
if 'int8' in data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
# 判读数组维数
if len(data.shape) == 3:
bands, height, width = data.shape
else:
bands = 1
height, width = data.shape
# 创建文件
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(filename, width, height, bands, datatype)
dataset.SetGeoTransform(geotrans)
dataset.SetProjection(proj)
if bands == 1:
dataset.GetRasterBand(1).WriteArray(data)
else:
for i in range(bands):
dataset.GetRasterBand(i + 1).WriteArray(data[i])
del dataset
# ############################ 遥感影像波段合成 #########################
def LayerStacking(inputPath,outfile_dir):
"""
将多光谱波段合成一个tif
:param dstlist: 输入但波段影像列表
:param outfile_dir: 影像的输出文件夹
"""
Rasters = glob.glob(os.path.join(inputPath + '\\' + "*B[1-7].TIF")) # *B[1-8].TIF
raster = Rasters[1]
a1, bandname = os.path.split(raster)
(name, ext) = os.path.splitext(bandname)
bandname = name[0:-6] # 提取文件名
dataset_init = gdal.Open(raster)
# 创建待输出的图
gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create(outfile_dir + "\\" + bandname + '.tif', dataset_init.RasterXSize, dataset_init.RasterYSize, 7, gdal.GDT_Float32)
# 本文合成的对象是4个波段，按照自己的需要改变
out_ds.SetProjection(dataset_init.GetProjection())
out_ds.SetGeoTransform(dataset_init.GetGeoTransform()) # 获得原始波段的地理信息
# 往栅格数据集中填各波段数值
for i in range(len(Rasters)):
dataset = gdal.Open(Rasters[i])
print('影像' + Rasters[i] + '开始波段合成')
band_temp = dataset.GetRasterBand(1)
band_temp.SetNoDataValue(-9999)
out_ds.GetRasterBand(i+1).WriteArray(band_temp.ReadAsArray()) # 注意band从1开始，所以要+1
# dataset.GetRasterBand(i + 1).WriteArray(data[i])
del out_ds
def GetSurfaceReflectance (coastal, blue, green, red, nir, swir1, swir2):
"""
得到地表发射率数据
:param coastal: landsat8的B1波段,蓝色波段blue
:param blue: landsat8的B2波段,蓝色波段blue
:param green: landsat8的B3波段,绿色波段green
:param red: landsat8的B4波段,红外波段red
:param nir: landsat8的B5波段,近红外波段nir
:param swir1: landsat8的B6波段,短波红外1 swir1
:param swir2: landsat8的B7波段,短波红外2 swir2
:return: 计算结果
"""
coastal = coastal * 0.0001 # 大气校正后的数据将各波段的值放大了10000倍数
blue = blue * 0.0001
green = green * 0.0001
red = red * 0.0001
nir = nir * 0.0001
swir1 = swir1 * 0.0001
swir2 = swir2 * 0.0001
data = np.array([coastal, blue, green, red, nir, swir1, swir2])
return data
if __name__ == '__main__':
print('开始Landsat波段合成')
InputPath = r"D:\algorithm\GEO\RSEI_code\data\Python_6s_AtmosphericCorrection"
Out = r"D:\algorithm\GEO\RSEI_code\data\Python_6s_AtmosphericCorrection\layerStacking"
totalout = r"D:\algorithm\GEO\RSEI_code\data\Python_6s_AtmosphericCorrection\layerStacking\SurfaceReflectance.tif"
LayerStacking(InputPath, Out)
# 读取波段数据
print('开始转换真实影像反射率')
Rasters = glob.glob(os.path.join(Out + '\\' + "*.tif"))
raster = Rasters[0]
im_data, width, height, proj, geotrans = Read_img(raster)
coastal = im_data[0]
blue = im_data[1]
green = im_data[2]
red = im_data[3]
nir = im_data[4]
swir1 = im_data[5]
swir2 = im_data[6]
data = GetSurfaceReflectance(coastal, blue, green, red, nir, swir1, swir2)
write_tiff(totalout, proj, geotrans, data)
print("影像波段合成及真是地表反射率计算完成")