用 python 对陆地卫星 8 数据进行辐射和大气校正
最编程
2024-04-04 12:41:53
...
# glob模块参考：
# glob是python自带的一个操作文件的相关模。用它可以查找符合特定规则的文件路径名。使用该模块查找文件，只需要用到： “*”, “?”, “[]”这三个匹配符;
import glob
import os
import sys
import tarfile
import re
from osgeo import gdal
import numpy
from Py6S import *
from osgeo import gdal
import pdb
import shutil
# argparse是一个Python模块：命令行选项、参数和子命令解析器。
import argparse
# 从base.py文件导入MeanDEM函数
from base import MeanDEM
def parse_arguments(argv):
# 此部分参考
# 此部分参考：
# 使用argparse的第一步是创建一个ArgumentParser对象。ArgumentParser对象包含将命令行解析成Python数据类型所需的全部信息。
parser = argparse.ArgumentParser() # 创建一个解析对象
# parser.add_argument() 向该对象中添加你要关注的命令行参数和选项
# name or flags - 一个命名或者一个选项字符串的列表；type - 命令行参数应当被转换成的类型；
# help - 一个此选项作用的简单描述；default - 当参数未在命令行中出现时使用的值。
# parser.add_argument('--Input_dir', type=str, help=r'D:\ProfessionalProfile\LandsatImage\LC08_L1TP_134036_20170808_20170813_01_T1\.', default=None)
# parser.add_argument('--Output_dir', type=str, help=r'D:\ProfessionalProfile\LandsatImage\1_OutputPy-RadiationCalibration0414\.', default=None)
parser.add_argument('--Input_dir', type=str, help='Input dir', default=r'D:\ProfessionalProfile\LandsatImage\LC08_L1TP_134036_20170808_20170813_01_T1')
parser.add_argument('--Output_dir', type=str, help='Output dir', default=r'D:\ProfessionalProfile\LandsatImage\1_OutputPy-RadiationCalibration0414')
# parser.add_argument('--Input_dir', type=str, help='Input dir', default=None)
# parser.add_argument('--Output_dir', type=str, help='Output dir', default=None)
# parse_args(args=None, nampespace=None)，args 参数名称，namespace 赋值。
return parser.parse_args(argv) # 进行解析
# 逐波段辐射定标
def RadiometricCalibration(BandId):
# LandSat8 TM辐射定标参数
global data2,ImgRasterData
parameter_OLI = numpy.zeros((9,2))
#计算辐射亮度参数
parameter_OLI[0,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_1.+',data2)[0]).split("=")[1])
parameter_OLI[1,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_2.+',data2)).split("=")[1])
parameter_OLI[2,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_3.+',data2)).split("=")[1])
parameter_OLI[3,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_4.+',data2)).split("=")[1])
parameter_OLI[4,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_5.+',data2)).split("=")[1])
parameter_OLI[5,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_6.+',data2)).split("=")[1])
parameter_OLI[6,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_7.+',data2)).split("=")[1])
parameter_OLI[7,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_8.+',data2)).split("=")[1])
parameter_OLI[8,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_9.+',data2)).split("=")[1])
parameter_OLI[0,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_1.+',data2)[0]).split("=")[1])
parameter_OLI[1,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_2.+',data2)).split("=")[1])
parameter_OLI[2,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_3.+',data2)).split("=")[1])
parameter_OLI[3,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_4.+',data2)).split("=")[1])
parameter_OLI[4,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_5.+',data2)).split("=")[1])
parameter_OLI[5,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_6.+',data2)).split("=")[1])
parameter_OLI[6,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_7.+',data2)).split("=")[1])
parameter_OLI[7,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_8.+',data2)).split("=")[1])
parameter_OLI[8,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_9.+',data2)).split("=")[1])
Gain = parameter_OLI[int(BandId) - 1,0]
Bias = parameter_OLI[int(BandId) - 1,1]
RaCal = numpy.where(ImgRasterData>0 ,Gain * ImgRasterData + Bias,-9999)
return (RaCal)
# 6s大气校正
def AtmosphericCorrection(BandId):
global data
# 6S模型
s = SixS()
s.geometry = Geometry.User()
s.geometry.solar_z = 90-float(''.join(re.findall('SUN_ELEVATION.+',data2)).split("=")[1])
s.geometry.solar_a = float(''.join(re.findall('SUN_AZIMUTH.+',data2)).split("=")[1])
s.geometry.view_z = 0
s.geometry.view_a = 0
# 日期
Dateparm = ''.join(re.findall('DATE_ACQUIRED.+',data2)).split("=")
Date = Dateparm[1].split('-')
s.geometry.month = int(Date[1])
s.geometry.day = int(Date[2])
# 中心经纬度
point1lat = float(''.join(re.findall('CORNER_UL_LAT_PRODUCT.+',data2)).split("=")[1])
point1lon = float(''.join(re.findall('CORNER_UL_LON_PRODUCT.+',data2)).split("=")[1])
point2lat = float(''.join(re.findall('CORNER_UR_LAT_PRODUCT.+',data2)).split("=")[1])
point2lon = float(''.join(re.findall('CORNER_UR_LON_PRODUCT.+',data2)).split("=")[1])
point3lat = float(''.join(re.findall('CORNER_LL_LAT_PRODUCT.+',data2)).split("=")[1])
point3lon = float(''.join(re.findall('CORNER_LL_LON_PRODUCT.+',data2)).split("=")[1])
point4lat = float(''.join(re.findall('CORNER_LR_LAT_PRODUCT.+',data2)).split("=")[1])
point4lon = float(''.join(re.findall('CORNER_LR_LON_PRODUCT.+',data2)).split("=")[1])
sLongitude = (point1lon + point2lon + point3lon + point4lon) / 4
sLatitude = (point1lat + point2lat + point3lat + point4lat) / 4
# 大气模式类型
if sLatitude > -15 and sLatitude <= 15:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical)
if sLatitude > 15 and sLatitude <= 45:
if s.geometry.month > 4 and s.geometry.month <= 9:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer)
else:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeWinter)
if sLatitude > 45 and sLatitude <= 60:
if s.geometry.month > 4 and s.geometry.month <= 9:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticSummer)
else:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticWinter)
# 气溶胶类型大陆
s.aero_profile = AtmosProfile.PredefinedType(AeroProfile.Continental)
# 目标地物？？？？？？
s.ground_reflectance = GroundReflectance.HomogeneousLambertian(0.36)
# 550nm气溶胶光学厚度,根据日期从MODIS处获取。
#s.visibility=40.0
s.aot550 = 0.14497
# 通过研究去区的范围去求DEM高度。
pointUL = dict()
pointDR = dict()
pointUL["lat"] = point1lat
pointUL["lon"] = point1lon
pointDR["lat"] = point4lat
pointDR["lon"] = point2lon
meanDEM = (MeanDEM(pointUL, pointDR)) * 0.001
# 研究区海拔、卫星传感器轨道高度
s.altitudes = Altitudes()
s.altitudes.set_target_custom_altitude(meanDEM)
s.altitudes.set_sensor_satellite_level()
# 校正波段（根据波段名称）
if BandId == '1':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B1)
elif BandId == '2':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B2)
elif BandId == '3':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B3)
elif BandId == '4':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B4)
elif BandId == '5':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B5)
elif BandId == '6':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B6)
elif BandId == '7':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B7)
elif BandId == '8':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B8)
elif BandId == '9':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B9)
# 下垫面非均一、朗伯体
s.atmos_corr = AtmosCorr.AtmosCorrLambertianFromReflectance(-0.1)
# 运行6s大气模型
s.run()
xa = s.outputs.coef_xa
xb = s.outputs.coef_xb
xc = s.outputs.coef_xc
x = s.outputs.values
return (xa, xb, xc)
if __name__ == '__main__':
#输入数据路径
# sys.argv参考
# sys.argv[]说白了就是一个从程序外部获取参数的桥梁，这个“外部”很关键，所以那些试图从代码来说明它作用的解释一直没看明白。因为我们从外部取得的参数可以是多个，所以获得的是一个列表（list)，也就是说sys.argv其实可以看作是一个列表，所以才能用[]提取其中的元素。其第一个元素是程序本身，随后才依次是外部给予的参数。
RootInputPath = parse_arguments(sys.argv[1:]).Input_dir
RootOutName = parse_arguments(sys.argv[2:]).Output_dir
#创建日志文件
LogFile = open(os.path.join(RootOutName,'log.txt'),'w')
# Python os.walk() 方法 https://www.runoob.com/python/os-walk.html
# os.walk(top[, topdown=True[, onerror=None[, followlinks=False]]])
# top -- 是你所要遍历的目录的地址, 返回的是一个三元组(root,dirs,files)。
# root 所指的是当前正在遍历的这个文件夹的本身的地址
# dirs 是一个 list ，内容是该文件夹中所有的目录的名字(不包括子目录)
# files 同样是 list , 内容是该文件夹中所有的文件(不包括子目录)
for root,dirs,RSFiles in os.walk(RootInputPath):
#判断是否进入最底层
if len(dirs)==0:
#根据输入输出路径建立生成新文件的路径
RootInputPathList = RootInputPath.split(os.path.sep)
RootList = root.split(os.path.sep)
StartList = len(RootInputPathList)
EndList = len(RootList)
outname = RootOutName
for i in range(StartList,EndList):
if os.path.exists(os.path.join(outname,RootList[i]))==False:
os.makedirs(os.path.join(outname,RootList[i]))
outname=os.path.join(outname,RootList[i])
else:
outname=os.path.join(outname,RootList[i])
MeteDatas = glob.glob(os.path.join(root,'*MTL.txt'))
for MeteData in MeteDatas:
pass
f = open(MeteData)
data = f.readlines()
data2 =' '.join(data)
shutil.copyfile(MeteData,os.path.join(outname,os.path.basename(MeteData)))
if len(os.path.basename(MeteData))<10:
RSbands = glob.glob(os.path.join(root,"B0[1-8].tiff"))
else:
RSbands = glob.glob(os.path.join(root,"*B[1-8].TIF"))
print('影像'+root+'开始大气校正')
print(RSbands)
for tifFile in RSbands:
BandId = (os.path.basename(tifFile).split('.')[0])[-1]
#捕捉打开数据出错异常
try:
IDataSet = gdal.Open(tifFile)
except Exception as e:
print("文件%S打开失败" % tifFile)
LogFile.write('\n'+tifFile+'数据打开失败')
if IDataSet == None:
LogFile.write('\n'+tifFile+'数据集读取为空')
continue
else:
#获取行列号
cols = IDataSet.RasterXSize
rows = IDataSet.RasterYSize
ImgBand = IDataSet.GetRasterBand(1)
ImgRasterData = ImgBand.ReadAsArray(0, 0, cols, rows)
if ImgRasterData is None:
LogFile.write('\n'+tifFile+'栅格数据为空')
continue
else:
#设置输出文件路径
outFilename=os.path.join(outname,os.path.basename(tifFile))
#如果文件存在就跳过，进行下一波段操作
if os.path.isfile(outFilename):
print("%s已经完成" % outFilename)
continue
else:
# #辐射校正
RaCalRaster = RadiometricCalibration(BandId)
#大气校正
a, b, c = AtmosphericCorrection(BandId)
y = numpy.where(RaCalRaster!=-9999,a * RaCalRaster - b,-9999)
atc = numpy.where(y!=-9999,(y / (1 + y * c))*10000,-9999)
driver = IDataSet.GetDriver()
#输出栅格数据集
outDataset = driver.Create(outFilename, cols, rows, 1, gdal.GDT_Int16)
# 设置投影信息，与原数据一样
geoTransform = IDataSet.GetGeoTransform()
outDataset.SetGeoTransform(geoTransform)
proj = IDataSet.GetProjection()
outDataset.SetProjection(proj)
outband = outDataset.GetRasterBand(1)
outband.SetNoDataValue(-9999)
outband.WriteArray(atc, 0, 0)
print('第%s波段计算完成'%BandId)
# print(root+'计算完成')
print('\n')
#关闭日志文件
LogFile.close()
# glob是python自带的一个操作文件的相关模。用它可以查找符合特定规则的文件路径名。使用该模块查找文件，只需要用到： “*”, “?”, “[]”这三个匹配符;
import glob
import os
import sys
import tarfile
import re
from osgeo import gdal
import numpy
from Py6S import *
from osgeo import gdal
import pdb
import shutil
# argparse是一个Python模块：命令行选项、参数和子命令解析器。
import argparse
# 从base.py文件导入MeanDEM函数
from base import MeanDEM
def parse_arguments(argv):
# 此部分参考
# 此部分参考：
# 使用argparse的第一步是创建一个ArgumentParser对象。ArgumentParser对象包含将命令行解析成Python数据类型所需的全部信息。
parser = argparse.ArgumentParser() # 创建一个解析对象
# parser.add_argument() 向该对象中添加你要关注的命令行参数和选项
# name or flags - 一个命名或者一个选项字符串的列表；type - 命令行参数应当被转换成的类型；
# help - 一个此选项作用的简单描述；default - 当参数未在命令行中出现时使用的值。
# parser.add_argument('--Input_dir', type=str, help=r'D:\ProfessionalProfile\LandsatImage\LC08_L1TP_134036_20170808_20170813_01_T1\.', default=None)
# parser.add_argument('--Output_dir', type=str, help=r'D:\ProfessionalProfile\LandsatImage\1_OutputPy-RadiationCalibration0414\.', default=None)
parser.add_argument('--Input_dir', type=str, help='Input dir', default=r'D:\ProfessionalProfile\LandsatImage\LC08_L1TP_134036_20170808_20170813_01_T1')
parser.add_argument('--Output_dir', type=str, help='Output dir', default=r'D:\ProfessionalProfile\LandsatImage\1_OutputPy-RadiationCalibration0414')
# parser.add_argument('--Input_dir', type=str, help='Input dir', default=None)
# parser.add_argument('--Output_dir', type=str, help='Output dir', default=None)
# parse_args(args=None, nampespace=None)，args 参数名称，namespace 赋值。
return parser.parse_args(argv) # 进行解析
# 逐波段辐射定标
def RadiometricCalibration(BandId):
# LandSat8 TM辐射定标参数
global data2,ImgRasterData
parameter_OLI = numpy.zeros((9,2))
#计算辐射亮度参数
parameter_OLI[0,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_1.+',data2)[0]).split("=")[1])
parameter_OLI[1,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_2.+',data2)).split("=")[1])
parameter_OLI[2,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_3.+',data2)).split("=")[1])
parameter_OLI[3,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_4.+',data2)).split("=")[1])
parameter_OLI[4,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_5.+',data2)).split("=")[1])
parameter_OLI[5,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_6.+',data2)).split("=")[1])
parameter_OLI[6,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_7.+',data2)).split("=")[1])
parameter_OLI[7,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_8.+',data2)).split("=")[1])
parameter_OLI[8,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_9.+',data2)).split("=")[1])
parameter_OLI[0,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_1.+',data2)[0]).split("=")[1])
parameter_OLI[1,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_2.+',data2)).split("=")[1])
parameter_OLI[2,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_3.+',data2)).split("=")[1])
parameter_OLI[3,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_4.+',data2)).split("=")[1])
parameter_OLI[4,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_5.+',data2)).split("=")[1])
parameter_OLI[5,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_6.+',data2)).split("=")[1])
parameter_OLI[6,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_7.+',data2)).split("=")[1])
parameter_OLI[7,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_8.+',data2)).split("=")[1])
parameter_OLI[8,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_9.+',data2)).split("=")[1])
Gain = parameter_OLI[int(BandId) - 1,0]
Bias = parameter_OLI[int(BandId) - 1,1]
RaCal = numpy.where(ImgRasterData>0 ,Gain * ImgRasterData + Bias,-9999)
return (RaCal)
# 6s大气校正
def AtmosphericCorrection(BandId):
global data
# 6S模型
s = SixS()
s.geometry = Geometry.User()
s.geometry.solar_z = 90-float(''.join(re.findall('SUN_ELEVATION.+',data2)).split("=")[1])
s.geometry.solar_a = float(''.join(re.findall('SUN_AZIMUTH.+',data2)).split("=")[1])
s.geometry.view_z = 0
s.geometry.view_a = 0
# 日期
Dateparm = ''.join(re.findall('DATE_ACQUIRED.+',data2)).split("=")
Date = Dateparm[1].split('-')
s.geometry.month = int(Date[1])
s.geometry.day = int(Date[2])
# 中心经纬度
point1lat = float(''.join(re.findall('CORNER_UL_LAT_PRODUCT.+',data2)).split("=")[1])
point1lon = float(''.join(re.findall('CORNER_UL_LON_PRODUCT.+',data2)).split("=")[1])
point2lat = float(''.join(re.findall('CORNER_UR_LAT_PRODUCT.+',data2)).split("=")[1])
point2lon = float(''.join(re.findall('CORNER_UR_LON_PRODUCT.+',data2)).split("=")[1])
point3lat = float(''.join(re.findall('CORNER_LL_LAT_PRODUCT.+',data2)).split("=")[1])
point3lon = float(''.join(re.findall('CORNER_LL_LON_PRODUCT.+',data2)).split("=")[1])
point4lat = float(''.join(re.findall('CORNER_LR_LAT_PRODUCT.+',data2)).split("=")[1])
point4lon = float(''.join(re.findall('CORNER_LR_LON_PRODUCT.+',data2)).split("=")[1])
sLongitude = (point1lon + point2lon + point3lon + point4lon) / 4
sLatitude = (point1lat + point2lat + point3lat + point4lat) / 4
# 大气模式类型
if sLatitude > -15 and sLatitude <= 15:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical)
if sLatitude > 15 and sLatitude <= 45:
if s.geometry.month > 4 and s.geometry.month <= 9:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer)
else:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeWinter)
if sLatitude > 45 and sLatitude <= 60:
if s.geometry.month > 4 and s.geometry.month <= 9:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticSummer)
else:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticWinter)
# 气溶胶类型大陆
s.aero_profile = AtmosProfile.PredefinedType(AeroProfile.Continental)
# 目标地物？？？？？？
s.ground_reflectance = GroundReflectance.HomogeneousLambertian(0.36)
# 550nm气溶胶光学厚度,根据日期从MODIS处获取。
#s.visibility=40.0
s.aot550 = 0.14497
# 通过研究去区的范围去求DEM高度。
pointUL = dict()
pointDR = dict()
pointUL["lat"] = point1lat
pointUL["lon"] = point1lon
pointDR["lat"] = point4lat
pointDR["lon"] = point2lon
meanDEM = (MeanDEM(pointUL, pointDR)) * 0.001
# 研究区海拔、卫星传感器轨道高度
s.altitudes = Altitudes()
s.altitudes.set_target_custom_altitude(meanDEM)
s.altitudes.set_sensor_satellite_level()
# 校正波段（根据波段名称）
if BandId == '1':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B1)
elif BandId == '2':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B2)
elif BandId == '3':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B3)
elif BandId == '4':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B4)
elif BandId == '5':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B5)
elif BandId == '6':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B6)
elif BandId == '7':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B7)
elif BandId == '8':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B8)
elif BandId == '9':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B9)
# 下垫面非均一、朗伯体
s.atmos_corr = AtmosCorr.AtmosCorrLambertianFromReflectance(-0.1)
# 运行6s大气模型
s.run()
xa = s.outputs.coef_xa
xb = s.outputs.coef_xb
xc = s.outputs.coef_xc
x = s.outputs.values
return (xa, xb, xc)
if __name__ == '__main__':
#输入数据路径
# sys.argv参考
# sys.argv[]说白了就是一个从程序外部获取参数的桥梁，这个“外部”很关键，所以那些试图从代码来说明它作用的解释一直没看明白。因为我们从外部取得的参数可以是多个，所以获得的是一个列表（list)，也就是说sys.argv其实可以看作是一个列表，所以才能用[]提取其中的元素。其第一个元素是程序本身，随后才依次是外部给予的参数。
RootInputPath = parse_arguments(sys.argv[1:]).Input_dir
RootOutName = parse_arguments(sys.argv[2:]).Output_dir
#创建日志文件
LogFile = open(os.path.join(RootOutName,'log.txt'),'w')
# Python os.walk() 方法 https://www.runoob.com/python/os-walk.html
# os.walk(top[, topdown=True[, onerror=None[, followlinks=False]]])
# top -- 是你所要遍历的目录的地址, 返回的是一个三元组(root,dirs,files)。
# root 所指的是当前正在遍历的这个文件夹的本身的地址
# dirs 是一个 list ，内容是该文件夹中所有的目录的名字(不包括子目录)
# files 同样是 list , 内容是该文件夹中所有的文件(不包括子目录)
for root,dirs,RSFiles in os.walk(RootInputPath):
#判断是否进入最底层
if len(dirs)==0:
#根据输入输出路径建立生成新文件的路径
RootInputPathList = RootInputPath.split(os.path.sep)
RootList = root.split(os.path.sep)
StartList = len(RootInputPathList)
EndList = len(RootList)
outname = RootOutName
for i in range(StartList,EndList):
if os.path.exists(os.path.join(outname,RootList[i]))==False:
os.makedirs(os.path.join(outname,RootList[i]))
outname=os.path.join(outname,RootList[i])
else:
outname=os.path.join(outname,RootList[i])
MeteDatas = glob.glob(os.path.join(root,'*MTL.txt'))
for MeteData in MeteDatas:
pass
f = open(MeteData)
data = f.readlines()
data2 =' '.join(data)
shutil.copyfile(MeteData,os.path.join(outname,os.path.basename(MeteData)))
if len(os.path.basename(MeteData))<10:
RSbands = glob.glob(os.path.join(root,"B0[1-8].tiff"))
else:
RSbands = glob.glob(os.path.join(root,"*B[1-8].TIF"))
print('影像'+root+'开始大气校正')
print(RSbands)
for tifFile in RSbands:
BandId = (os.path.basename(tifFile).split('.')[0])[-1]
#捕捉打开数据出错异常
try:
IDataSet = gdal.Open(tifFile)
except Exception as e:
print("文件%S打开失败" % tifFile)
LogFile.write('\n'+tifFile+'数据打开失败')
if IDataSet == None:
LogFile.write('\n'+tifFile+'数据集读取为空')
continue
else:
#获取行列号
cols = IDataSet.RasterXSize
rows = IDataSet.RasterYSize
ImgBand = IDataSet.GetRasterBand(1)
ImgRasterData = ImgBand.ReadAsArray(0, 0, cols, rows)
if ImgRasterData is None:
LogFile.write('\n'+tifFile+'栅格数据为空')
continue
else:
#设置输出文件路径
outFilename=os.path.join(outname,os.path.basename(tifFile))
#如果文件存在就跳过，进行下一波段操作
if os.path.isfile(outFilename):
print("%s已经完成" % outFilename)
continue
else:
# #辐射校正
RaCalRaster = RadiometricCalibration(BandId)
#大气校正
a, b, c = AtmosphericCorrection(BandId)
y = numpy.where(RaCalRaster!=-9999,a * RaCalRaster - b,-9999)
atc = numpy.where(y!=-9999,(y / (1 + y * c))*10000,-9999)
driver = IDataSet.GetDriver()
#输出栅格数据集
outDataset = driver.Create(outFilename, cols, rows, 1, gdal.GDT_Int16)
# 设置投影信息，与原数据一样
geoTransform = IDataSet.GetGeoTransform()
outDataset.SetGeoTransform(geoTransform)
proj = IDataSet.GetProjection()
outDataset.SetProjection(proj)
outband = outDataset.GetRasterBand(1)
outband.SetNoDataValue(-9999)
outband.WriteArray(atc, 0, 0)
print('第%s波段计算完成'%BandId)
# print(root+'计算完成')
print('\n')
#关闭日志文件
LogFile.close()
推荐阅读
-
[姿势估计] 实践记录：使用 Dlib 和 mediapipe 进行人脸姿势估计 - 本文重点介绍方法 2）：方法 1：基于深度学习的方法：。 基于深度学习的方法：基于深度学习的方法利用深度学习模型，如卷积神经网络（CNN）或递归神经网络（RNN），直接从人脸图像中学习姿势估计。这些方法能够学习更复杂的特征表征，并在大规模数据集上取得优异的性能。方法二：基于二维校准信息估计三维姿态信息（计算机视觉 PnP 问题）。 特征点定位：人脸姿态估计的第一步是通过特征点定位来检测和定位人脸的关键点，如眼睛、鼻子和嘴巴。这些关键点提供了人脸的局部结构信息，可用于后续的姿势估计。 旋转表示：常见的旋转表示方法包括欧拉角和旋转矩阵。欧拉角通过三个旋转角度（通常是俯仰、偏航和滚动）描述头部的旋转姿态。旋转矩阵是一个 3x3 矩阵，表示头部从一个坐标系到另一个坐标系的变换。 三维模型重建：根据特征点的定位结果，三维人脸模型可用于姿势估计。通过将人脸的二维图像映射到三维模型上，可以估算出人脸的旋转和平移信息。这就需要建立人脸的三维模型，然后通过优化方法将模型与特征点对齐，从而获得姿势估计结果。 特征点定位 特征点定位是用于检测人脸关键部位的五官基础部分，还有其他更多的特征点表示方法，大家可以参考我上一篇文章中介绍的特征点检测方案实践：人脸校正二次定位操作来解决人脸校正的问题，客户在检测关键点的代码上略有修改，坐标转换部分客户见上图 def get_face_info(image). img_copy = image.copy image.flags.writeable = False image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB) results = face_detection.process(image) # 在图像上绘制人脸检测注释。 image.flags.writeable = True image = cv2.cvtColor(image, cv2.COLOR_RGB2BGR) box_info, facial = None, None if results.detections： for detection in results. for detection in results.detections: mp_drawing.Drawing.detection = 无 mp_drawing.draw_detection(image, detection) 面部 = detection.location_data.relative_keypoints 返回面部 在上述代码中，返回的数据是五官（6 个关键点的坐标），这是用 mediapipe 库实现的，下面我们可以尝试用另一个库：dlib 来实现。 使用 dlib 使用 Dlib 库在 Python 中实现人脸关键点检测的步骤如下： 确保已安装 Dlib 库，可使用以下命令： pip install dlib 导入必要的库： 加载 Dlib 的人脸检测器和关键点检测器模型： 读取图像并将其灰度化： 使用人脸检测器检测图像中的人脸： 对检测到的人脸进行遍历，并使用关键点检测器检测人脸关键点： 显示绘制了关键点的图像： 以下代码将参数 landmarks_part 添加到要返回的关键点坐标中。
-
用 python 对陆地卫星 8 数据进行辐射和大气校正