校正GF-2数据后，水体的校正结果与FLAASH有较大差异

[alive9999]

校正后的GF-2数据，植被的地表反射率光谱曲线与FLAASH校正曲线基本吻合，但是水体的光谱曲线差异较大，有些水体的近红外波段校正后的数据偏高。

[kymi014]

你好，请问你的问题解决了吗，能否将你的代码发送给我，我下载这个作者的代码运行的影像没有坐标

[alive9999]

你好，这个问题没有解决。 关于影像坐标的问题，我先给正射校正后的影像加上了投影，然后再拿带有坐标的影像进行的大气校正。 附件是我的代码，供您参考。  

…

------------------ 原始邮件 ------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***>; 发送时间: 2021年12月22日(星期三) 下午4:05 ***@***.***>; ***@***.******@***.***>; 主题: Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 你好，请问你的问题解决了吗，能否将你的代码发送给我，我下载这个作者的代码运行的影像没有坐标 — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***>

[kymi014]

抱歉，我没有看到附件哎

…

------------------ 原始邮件 ------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***>; 发送时间: 2021年12月22日(星期三) 下午4:12 ***@***.***>; ***@***.******@***.***>; 主题: Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 你好，这个问题没有解决。 关于影像坐标的问题，我先给正射校正后的影像加上了投影，然后再拿带有坐标的影像进行的大气校正。 附件是我的代码，供您参考。  

------------------ 原始邮件 ------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***>; 发送时间: 2021年12月22日(星期三) 下午4:05 ***@***.***>; ***@***.******@***.***>; 主题: Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 你好，请问你的问题解决了吗，能否将你的代码发送给我，我下载这个作者的代码运行的影像没有坐标 — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***> — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.Message ID: ***@***.***>

[kymi014]

那请问你用高分原始影像跑的话，在没有正射校正的情况下有坐标吗

…

------------------ 原始邮件 ------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***>; 发送时间: 2021年12月22日(星期三) 下午4:12 ***@***.***>; ***@***.******@***.***>; 主题: Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 你好，这个问题没有解决。 关于影像坐标的问题，我先给正射校正后的影像加上了投影，然后再拿带有坐标的影像进行的大气校正。 附件是我的代码，供您参考。  

------------------ 原始邮件 ------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***>; 发送时间: 2021年12月22日(星期三) 下午4:05 ***@***.***>; ***@***.******@***.***>; 主题: Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 你好，请问你的问题解决了吗，能否将你的代码发送给我，我下载这个作者的代码运行的影像没有坐标 — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***> — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.Message ID: ***@***.***>

[alive9999]

我拿到的影像是经过正射校正后的，不太清楚没有正射校正的影像是什么情况。这次附件里可以看到了吗？  

…

------------------ 原始邮件 ------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***>; 发送时间: 2021年12月22日(星期三) 下午4:15 ***@***.***>; ***@***.******@***.***>; 主题: Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 那请问你用高分原始影像跑的话，在没有正射校正的情况下有坐标吗

------------------ 原始邮件 ------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***>; 发送时间: 2021年12月22日(星期三) 下午4:12 ***@***.***>; ***@***.******@***.***>; 主题: Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 你好，这个问题没有解决。 关于影像坐标的问题，我先给正射校正后的影像加上了投影，然后再拿带有坐标的影像进行的大气校正。 附件是我的代码，供您参考。  

------------------ 原始邮件 ------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***>; 发送时间: 2021年12月22日(星期三) 下午4:05 ***@***.***>; ***@***.******@***.***>; 主题: Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 你好，请问你的问题解决了吗，能否将你的代码发送给我，我下载这个作者的代码运行的影像没有坐标 — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***> — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.Message ID: ***@***.***> — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***>

[kymi014]

抱歉啊，还是看不见，看来屏蔽了，你能看到我的扣扣吗1310118141，如果可以麻烦你可以发到我邮箱吗万分感谢

…

------------------ 原始邮件 ------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***>; 发送时间: 2021年12月22日(星期三) 下午4:22 ***@***.***>; ***@***.******@***.***>; 主题: Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 我拿到的影像是经过正射校正后的，不太清楚没有正射校正的影像是什么情况。这次附件里可以看到了吗？  

------------------ 原始邮件 ------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***>; 发送时间: 2021年12月22日(星期三) 下午4:15 ***@***.***>; ***@***.******@***.***>; 主题: Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 那请问你用高分原始影像跑的话，在没有正射校正的情况下有坐标吗

------------------ 原始邮件 ------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***>; 发送时间: 2021年12月22日(星期三) 下午4:12 ***@***.***>; ***@***.******@***.***>; 主题: Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 你好，这个问题没有解决。 关于影像坐标的问题，我先给正射校正后的影像加上了投影，然后再拿带有坐标的影像进行的大气校正。 附件是我的代码，供您参考。  

------------------ 原始邮件 ------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***>; 发送时间: 2021年12月22日(星期三) 下午4:05 ***@***.***>; ***@***.******@***.***>; 主题: Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 你好，请问你的问题解决了吗，能否将你的代码发送给我，我下载这个作者的代码运行的影像没有坐标 — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***> — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.Message ID: ***@***.***> — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***> — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.Message ID: ***@***.***>

[alive9999]

我直接粘贴到邮件里吧~ #! usr/bin/env python # -*- coding:utf-8 -*- # created by zhaoguanhua 2017/09/04 import glob import os import sys import tarfile&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; #解压缩 import json import numpy as np import gdal import pdb import math import time import xml.dom.minidom&nbsp; &nbsp; #读取xml格式的影像头文件 from tqdm import tqdm&nbsp; &nbsp; &nbsp;#进度条 from Py6S import * import argparse from base import MeanDEM import shutil def parse_arguments(argv): &nbsp; &nbsp; parser = argparse.ArgumentParser() &nbsp; &nbsp; parser.add_argument('--Input_dir',type=str,help='Input dir',default=None) &nbsp; &nbsp; parser.add_argument('--Output_dir',type=str,help='Output dir',default=None) &nbsp; &nbsp; return parser.parse_args(argv) def search_file(path,str): &nbsp; &nbsp; for x in os.listdir(path): &nbsp; &nbsp; &nbsp; &nbsp; next_path = os.path.join(path,x) &nbsp; &nbsp; &nbsp; &nbsp; if os.path.isfile(next_path): &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; if str in x: &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; return x # 解压缩原始文件 def untar(fname, dirs): &nbsp; &nbsp; print("文件路径",fname) &nbsp; &nbsp; try: &nbsp; &nbsp; &nbsp; &nbsp; t = tarfile.open(fname) &nbsp; &nbsp; except Exception as e: &nbsp; &nbsp; &nbsp; &nbsp; print("文件%s打开失败" % fname) &nbsp; &nbsp; t.extractall(path=dirs) def Block(IDataSet,filename_split,outPath_file,ImageType,config,metedata): #def Block(IDataSet,ImageType,config,outFileName): &nbsp; &nbsp; # global cols,rows,atcfiles &nbsp; &nbsp; cols = IDataSet.RasterXSize &nbsp; &nbsp; rows = IDataSet.RasterYSize &nbsp; &nbsp; SatelliteID = filename_split[0] &nbsp; &nbsp; SensorID = filename_split[1] &nbsp; &nbsp; Year = filename_split[4][:4] &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; #设置输出波段 &nbsp; &nbsp; Driver = IDataSet.GetDriver() &nbsp; &nbsp; geoTransform1 = IDataSet.GetGeoTransform() &nbsp; &nbsp; ListgeoTransform1 = list(geoTransform1) &nbsp; &nbsp; ListgeoTransform1[5] = -ListgeoTransform1[5] &nbsp; &nbsp; newgeoTransform1 = tuple(ListgeoTransform1) &nbsp; &nbsp; proj1 = IDataSet.GetProjection() &nbsp; &nbsp; #OutRCname = os.path.join(atcfiles,outFileName+".tiff") &nbsp; &nbsp; #OutRCname = 'data/gf2_0503_rf.tif' &nbsp; &nbsp; OutRCname = outPath_file &nbsp; &nbsp; outDataset = Driver.Create(OutRCname,cols,rows,4,gdal.GDT_Int32) &nbsp; &nbsp; #outDataset.SetGeoTransform(newgeoTransform1) &nbsp; &nbsp; outDataset.SetGeoTransform(geoTransform1) &nbsp; &nbsp; outDataset.SetProjection(proj1) &nbsp; &nbsp; #分别读取4个波段 &nbsp; &nbsp; for m in range(1,5): &nbsp; &nbsp; &nbsp; &nbsp; ReadBand = IDataSet.GetRasterBand(m) &nbsp; &nbsp; &nbsp; &nbsp; outband = outDataset.GetRasterBand(m) &nbsp; &nbsp; &nbsp; &nbsp; outband.SetNoDataValue(-9999) &nbsp; &nbsp; &nbsp; &nbsp; #获取对应波段的增益gain和偏移bias &nbsp; &nbsp; &nbsp; &nbsp; Gain,Bias = RadiometricCalibration(m,SatelliteID,SensorID,Year,ImageType,config) &nbsp; &nbsp; &nbsp; &nbsp; #获取大气校正系数 &nbsp; &nbsp; &nbsp; &nbsp; AtcCofa, AtcCofb, AtcCofc = AtmosphericCorrection(m,metedata,config,SatelliteID,SensorID) &nbsp; &nbsp; &nbsp; &nbsp; nBlockSize = 2048 &nbsp; &nbsp; &nbsp; &nbsp; i = 0 &nbsp; &nbsp; &nbsp; &nbsp; j = 0 &nbsp; &nbsp; &nbsp; &nbsp; b = cols*rows &nbsp; &nbsp; &nbsp; &nbsp; #进度条参数 &nbsp; &nbsp; &nbsp; &nbsp; XBlockcount = math.ceil(cols/nBlockSize) &nbsp; &nbsp; &nbsp; &nbsp; YBlockcount = math.ceil(rows/nBlockSize) &nbsp; &nbsp; &nbsp; &nbsp; print("第%d波段校正："%m) &nbsp; &nbsp; &nbsp; &nbsp; try: &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; with tqdm(total=XBlockcount*YBlockcount,iterable='iterable',desc = '第%i波段:'%m,mininterval=10) as pbar: &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; # with tqdm(total=XBlockcount*YBlockcount) as pbar: &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; # print(pbar) &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; while i<rows: &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; while j <cols: &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; #保存分块大小 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; nXBK = nBlockSize &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; nYBK = nBlockSize &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; #最后不够分块的区域，有多少读取多少 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; if i+nBlockSize&gt;rows: &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; nYBK = rows - i &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; if j+nBlockSize&gt;cols: &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; nXBK=cols - j &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; #分块读取影像 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Image = ReadBand.ReadAsArray(j,i,nXBK,nYBK) &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; outImage =np.where(Image&gt;0,Image*Gain + Bias,-9999) &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; y = np.where(outImage!=-9999,AtcCofa * outImage - AtcCofb,-9999) &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; atcImage = np.where(y!=-9999,(y / (1 + y * AtcCofc))*10000,-9999) &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; outband.WriteArray(atcImage,j,i) &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; j=j+nXBK &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; time.sleep(1) &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; pbar.update(1) &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; j=0 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; i=i+nYBK &nbsp; &nbsp; &nbsp; &nbsp; except KeyboardInterrupt: &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; pbar.close() &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; raise &nbsp; &nbsp; &nbsp; &nbsp; pbar.close() def RadiometricCalibration(BandId,SatelliteID,SensorID,Year,ImageType,config): &nbsp; &nbsp; # global cols,rows,SatelliteID,SensorID,Year,ImageType,config &nbsp; &nbsp; if SensorID[0:3] == "WFV": &nbsp; &nbsp; &nbsp; &nbsp; Gain_ =config["Parameter"][SatelliteID][SensorID][Year]["gain"][BandId-1] &nbsp; &nbsp; &nbsp; &nbsp; Bias_ =config["Parameter"][SatelliteID][SensorID][Year]["offset"][BandId-1] &nbsp; &nbsp; else: &nbsp; &nbsp; &nbsp; &nbsp; Gain_ =config["Parameter"][SatelliteID][SensorID][Year][ImageType]["gain"][BandId-1] &nbsp; &nbsp; &nbsp; &nbsp; Bias_ =config["Parameter"][SatelliteID][SensorID][Year][ImageType]["offset"][BandId-1] &nbsp; &nbsp; return Gain_,Bias_ # 6s大气校正 #def AtmosphericCorrection(BandId,config,SatelliteID,SensorID): def AtmosphericCorrection(BandId,metedata,config,SatelliteID,SensorID): &nbsp; &nbsp; # global metedata,config,SatelliteID,SensorID &nbsp; &nbsp; #读取头文件 &nbsp; &nbsp; dom = xml.dom.minidom.parse(metedata) &nbsp; &nbsp; # 6S模型 &nbsp; &nbsp; s = SixS() &nbsp; &nbsp; # 传感器类型 自定义 &nbsp; &nbsp; s.geometry = Geometry.User() &nbsp; &nbsp; #s.geometry.solar_z = 90-float(dom.getElementsByTagName('SolarZenith')[0].firstChild.data) &nbsp; &nbsp; s.geometry.solar_z = float(dom.getElementsByTagName('SolarZenith')[0].firstChild.data) &nbsp; &nbsp; s.geometry.solar_a = float(dom.getElementsByTagName('SolarAzimuth')[0].firstChild.data) &nbsp; &nbsp; #s.geometry.solar_a = 147.636 &nbsp; &nbsp; #s.geometry.solar_z = 27.5293 &nbsp; &nbsp; #s.geometry.solar_a = 146.242 &nbsp; &nbsp; #s.geometry.solar_z = 27.8728 &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; #s.geometry.view_z = float(dom.getElementsByTagName('SatelliteZenith')[0].firstChild.data) &nbsp; &nbsp; #s.geometry.view_a = float(dom.getElementsByTagName('SatelliteAzimuth')[0].firstChild.data) &nbsp; &nbsp; s.geometry.view_z = 0 &nbsp; &nbsp; s.geometry.view_a = 0 &nbsp; &nbsp; # 日期 &nbsp; &nbsp; DateTimeparm = dom.getElementsByTagName('CenterTime')[0].firstChild.data &nbsp; &nbsp; DateTime = DateTimeparm.split(' ') &nbsp; &nbsp; Date = DateTime[0].split('-') &nbsp; &nbsp; s.geometry.month = int(Date[1]) &nbsp; &nbsp; s.geometry.day = int(Date[2]) &nbsp; &nbsp; #s.geometry.month = 5 &nbsp; &nbsp; #s.geometry.day = 3 &nbsp; &nbsp; # print(s.geometry) &nbsp; &nbsp; # 中心经纬度 &nbsp; &nbsp; TopLeftLat = float(dom.getElementsByTagName('TopLeftLatitude')[0].firstChild.data) &nbsp; &nbsp; TopLeftLon = float(dom.getElementsByTagName('TopLeftLongitude')[0].firstChild.data) &nbsp; &nbsp; TopRightLat = float(dom.getElementsByTagName('TopRightLatitude')[0].firstChild.data) &nbsp; &nbsp; TopRightLon = float(dom.getElementsByTagName('TopRightLongitude')[0].firstChild.data) &nbsp; &nbsp; BottomRightLat = float(dom.getElementsByTagName('BottomRightLatitude')[0].firstChild.data) &nbsp; &nbsp; BottomRightLon = float(dom.getElementsByTagName('BottomRightLongitude')[0].firstChild.data) &nbsp; &nbsp; BottomLeftLat = float(dom.getElementsByTagName('BottomLeftLatitude')[0].firstChild.data) &nbsp; &nbsp; BottomLeftLon = float(dom.getElementsByTagName('BottomLeftLongitude')[0].firstChild.data) &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; #TopLeftLat = 39.848448 &nbsp; &nbsp; #TopLeftLon = 115.244064 &nbsp; &nbsp; #TopRightLat = 39.848448 &nbsp; &nbsp; #TopRightLon = 115.569216 &nbsp; &nbsp; #BottomRightLat = 39.604288 &nbsp; &nbsp; #BottomRightLon = 115.569216 &nbsp; &nbsp; #BottomLeftLat = 39.604288 &nbsp; &nbsp; #BottomLeftLon = 115.244064 &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; # TopLeftLat = 40.15664 &nbsp; &nbsp; # TopLeftLon = 116.033536 &nbsp; &nbsp; # TopRightLat = 40.15664 &nbsp; &nbsp; # TopRightLon = 116.356896 &nbsp; &nbsp; # BottomRightLat = 39.914944 &nbsp; &nbsp; # BottomRightLon = 116.356896 &nbsp; &nbsp; # BottomLeftLat = 39.914944 &nbsp; &nbsp; # BottomLeftLon = 116.033536 &nbsp; &nbsp; ImageCenterLat = (TopLeftLat + TopRightLat + BottomRightLat + BottomLeftLat) / 4 &nbsp; &nbsp; # 大气模式类型 &nbsp; &nbsp; if ImageCenterLat &gt; -15 and ImageCenterLat < 15: &nbsp; &nbsp; &nbsp; &nbsp; s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical) &nbsp; &nbsp; if ImageCenterLat &gt; 15 and ImageCenterLat < 45: &nbsp; &nbsp; &nbsp; &nbsp; if s.geometry.month &gt; 4 and s.geometry.month < 9: &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer) &nbsp; &nbsp; &nbsp; &nbsp; else: &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeWinter) &nbsp; &nbsp; if ImageCenterLat &gt; 45 and ImageCenterLat < 60: &nbsp; &nbsp; &nbsp; &nbsp; if s.geometry.month &gt; 4 and s.geometry.month < 9: &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticSummer) &nbsp; &nbsp; &nbsp; &nbsp; else: &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticWinter) &nbsp; &nbsp; # 气溶胶类型大陆 &nbsp; &nbsp; s.aero_profile = AtmosProfile.PredefinedType(AeroProfile.Continental) &nbsp; &nbsp; # 下垫面类型 &nbsp; &nbsp; s.ground_reflectance = GroundReflectance.HomogeneousLambertian(0.36) &nbsp; &nbsp; #s.ground_reflectance = GroundReflectance.HomogeneousLambertian(GroundReflectance.ClearWater) &nbsp; &nbsp; # 550nm气溶胶光学厚度,对应能见度为40km &nbsp; &nbsp; s.aot550 = 0.14497 &nbsp; &nbsp; #s.aot550 = 40 &nbsp; &nbsp; # 通过研究去区的范围去求DEM高度。 &nbsp; &nbsp; pointUL = dict() &nbsp; &nbsp; pointDR = dict() &nbsp; &nbsp; pointUL["lat"] = max(TopLeftLat,TopRightLat,BottomRightLat,BottomLeftLat) &nbsp; &nbsp; pointUL["lon"] = min(TopLeftLon,TopRightLon,BottomRightLon,BottomLeftLon) &nbsp; &nbsp; pointDR["lat"] = min(TopLeftLat,TopRightLat,BottomRightLat,BottomLeftLat) &nbsp; &nbsp; pointDR["lon"] = max(TopLeftLon,TopRightLon,BottomRightLon,BottomLeftLon) &nbsp; &nbsp; #meanDEM = (MeanDEM(pointUL, pointDR)) * 0.001 &nbsp; &nbsp; meanDEM = 0.0 &nbsp; &nbsp; # 研究区海拔、卫星传感器轨道高度 &nbsp; &nbsp; s.altitudes = Altitudes() &nbsp; &nbsp; s.altitudes.set_target_custom_altitude(meanDEM) &nbsp; &nbsp; s.altitudes.set_sensor_satellite_level() &nbsp; &nbsp; # 校正波段（根据波段名称） &nbsp; &nbsp; if BandId == 1: &nbsp; &nbsp; &nbsp; &nbsp; SRFband = config["Parameter"][SatelliteID][SensorID]["SRF"]["1"] &nbsp; &nbsp; &nbsp; &nbsp; s.wavelength = Wavelength(0.450,0.520,SRFband) &nbsp; &nbsp; elif BandId == 2: &nbsp; &nbsp; &nbsp; &nbsp; SRFband = config["Parameter"][SatelliteID][SensorID]["SRF"]["2"] &nbsp; &nbsp; &nbsp; &nbsp; s.wavelength = Wavelength(0.520,0.590,SRFband) &nbsp; &nbsp; elif BandId == 3: &nbsp; &nbsp; &nbsp; &nbsp; SRFband = config["Parameter"][SatelliteID][SensorID]["SRF"]["3"] &nbsp; &nbsp; &nbsp; &nbsp; s.wavelength = Wavelength(0.630,0.690,SRFband) &nbsp; &nbsp; elif BandId == 4: &nbsp; &nbsp; &nbsp; &nbsp; SRFband = config["Parameter"][SatelliteID][SensorID]["SRF"]["4"] &nbsp; &nbsp; &nbsp; &nbsp; s.wavelength = Wavelength(0.770,0.890,SRFband) &nbsp; &nbsp; s.atmos_corr = AtmosCorr.AtmosCorrLambertianFromReflectance(-0.1) &nbsp; &nbsp; # 运行6s大气模型 &nbsp; &nbsp; s.run() &nbsp; &nbsp; xa = s.outputs.coef_xa &nbsp; &nbsp; xb = s.outputs.coef_xb &nbsp; &nbsp; xc = s.outputs.coef_xc &nbsp; &nbsp; # x = s.outputs.values &nbsp; &nbsp; return (xa, xb, xc) if __name__ == '__main__': &nbsp; &nbsp; a=time.time() &nbsp; &nbsp; script_path = os.path.split(os.path.realpath(__file__))[0] &nbsp; &nbsp; #读取辐射校正和大气校正所需参数:增益、偏移和光谱响应函数 &nbsp; &nbsp; config_file = os.path.join(script_path,"RadiometricCorrectionParameter.json") &nbsp; &nbsp; config = json.load(open(config_file)) &nbsp; &nbsp; #输入数据路径 &nbsp; &nbsp; #InputFilePath = parse_arguments(sys.argv[1:]).Input_dir &nbsp; &nbsp; #OutputFilePath = parse_arguments(sys.argv[2:]).Output_dir &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; #头文件文件夹 &nbsp; &nbsp; xml_path = r'F:\2-Project-data\BeijingEco\xml' &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; #遍历文件夹 &nbsp; &nbsp; input_path = r'F:\2-Project-data\BeijingEco\output_proj' &nbsp; &nbsp; #input_path = r'F:\2-Project-data\BeijingEco\test_ac' &nbsp; &nbsp; files = os.listdir(input_path) &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; output_path = r'F:\2-Project-data\BeijingEco\atmospheric_correction' &nbsp; &nbsp; #file = 'GF2_PMS2_E116.2_N40.0_20210503_L1A0005629120-MSS2_ORTHO_MS.tif' &nbsp; &nbsp; for file in files: &nbsp; &nbsp; &nbsp; &nbsp; #file:'GF2_PMS1_E115.6_N39.9_20200602_L1A0004842146-MSS1_ORTHO_MS_proj.tif'，含有tif &nbsp; &nbsp; &nbsp; &nbsp; if os.path.splitext(file)[1].upper() == ".TIF": &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; print(file) &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; tiffFile = input_path + os.sep + file &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; IDataSet = gdal.Open(tiffFile) &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; ImageType = 'MSS' &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; filename = os.path.splitext(file)[0]&nbsp; #'GF2_PMS1_E115.6_N39.9_20200602_L1A0004842146-MSS1_ORTHO_MS_proj' &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; fileType = filename[0:2] &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; filename_split = filename.split("_") &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; outFileName = filename[:-4]+'RF.tif' &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; outPath_file = os.path.join(output_path, outFileName) &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; GFType = filename_split[1][:3] #GF2 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; #读取头文件 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; productId = filename_split[5] &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; xml_file = search_file(xml_path,productId) &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; metedata = glob.glob(os.path.join(xml_path,xml_file))[0] &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; #metedata = glob.glob(os.path.join(outname,"*MSS*.xml"))[0] &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; #outFileName = os.path.splitext(file)[0] &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; #outFileName = output_path+os.sep+ &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; #Block(IDataSet,filename_split,outPath_file,ImageType,config,metedata) &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Block(IDataSet,filename_split,outPath_file,ImageType,config,metedata) &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp; b=time.time() &nbsp; &nbsp; print("总时间:",b-a) &nbsp; &nbsp;&nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp;&nbsp; &nbsp;

…

------------------&nbsp;原始邮件&nbsp;------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***&gt;; 发送时间:&nbsp;2021年12月22日(星期三) 下午4:24 ***@***.***&gt;; ***@***.******@***.***&gt;; 主题:&nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 抱歉啊，还是看不见，看来屏蔽了，你能看到我的扣扣吗1310118141，如果可以麻烦你可以发到我邮箱吗万分感谢

------------------&amp;nbsp;原始邮件&amp;nbsp;------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***&amp;gt;; 发送时间:&amp;nbsp;2021年12月22日(星期三) 下午4:22 ***@***.***&amp;gt;; ***@***.******@***.***&amp;gt;; 主题:&amp;nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 我拿到的影像是经过正射校正后的，不太清楚没有正射校正的影像是什么情况。这次附件里可以看到了吗？ &amp;amp;nbsp;

------------------&amp;amp;nbsp;原始邮件&amp;amp;nbsp;------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***&amp;amp;gt;; 发送时间:&amp;amp;nbsp;2021年12月22日(星期三) 下午4:15 ***@***.***&amp;amp;gt;; ***@***.******@***.***&amp;amp;gt;; 主题:&amp;amp;nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 那请问你用高分原始影像跑的话，在没有正射校正的情况下有坐标吗

------------------&amp;amp;amp;nbsp;原始邮件&amp;amp;amp;nbsp;------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***&amp;amp;amp;gt;; 发送时间:&amp;amp;amp;nbsp;2021年12月22日(星期三) 下午4:12 ***@***.***&amp;amp;amp;gt;; ***@***.******@***.***&amp;amp;amp;gt;; 主题:&amp;amp;amp;nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 你好，这个问题没有解决。 关于影像坐标的问题，我先给正射校正后的影像加上了投影，然后再拿带有坐标的影像进行的大气校正。 附件是我的代码，供您参考。 &amp;amp;amp;amp;nbsp;

------------------&amp;amp;amp;amp;nbsp;原始邮件&amp;amp;amp;amp;nbsp;------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***&amp;amp;amp;amp;gt;; 发送时间:&amp;amp;amp;amp;nbsp;2021年12月22日(星期三) 下午4:05 ***@***.***&amp;amp;amp;amp;gt;; ***@***.******@***.***&amp;amp;amp;amp;gt;; 主题:&amp;amp;amp;amp;nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 你好，请问你的问题解决了吗，能否将你的代码发送给我，我下载这个作者的代码运行的影像没有坐标 — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***&amp;amp;amp;amp;gt; — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.Message ID: ***@***.***&amp;amp;amp;gt; — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***&amp;amp;gt; — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.Message ID: ***@***.***&amp;gt; — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***&gt;

[kymi014]

抱歉，是我的号被屏蔽了吗？我这边看你复制的内容很大部分都被屏蔽了

…

------------------&nbsp;原始邮件&nbsp;------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***&gt;; 发送时间:&nbsp;2021年12月22日(星期三) 下午4:26 ***@***.***&gt;; ***@***.******@***.***&gt;; 主题:&nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 我直接粘贴到邮件里吧~ #! usr/bin/env python # -*- coding:utf-8 -*- # created by zhaoguanhua 2017/09/04 import glob import os import sys import tarfile&amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; #解压缩 import json import numpy as np import gdal import pdb import math import time import xml.dom.minidom&amp;nbsp; &amp;nbsp; #读取xml格式的影像头文件 from tqdm import tqdm&amp;nbsp; &amp;nbsp; &amp;nbsp;#进度条 from Py6S import * import argparse from base import MeanDEM import shutil def parse_arguments(argv): &amp;nbsp; &amp;nbsp; parser = argparse.ArgumentParser() &amp;nbsp; &amp;nbsp; parser.add_argument('--Input_dir',type=str,help='Input dir',default=None) &amp;nbsp; &amp;nbsp; parser.add_argument('--Output_dir',type=str,help='Output dir',default=None) &amp;nbsp; &amp;nbsp; return parser.parse_args(argv) def search_file(path,str): &amp;nbsp; &amp;nbsp; for x in os.listdir(path): &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; next_path = os.path.join(path,x) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; if os.path.isfile(next_path): &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; if str in x: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; return x # 解压缩原始文件 def untar(fname, dirs): &amp;nbsp; &amp;nbsp; print("文件路径",fname) &amp;nbsp; &amp;nbsp; try: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; t = tarfile.open(fname) &amp;nbsp; &amp;nbsp; except Exception as e: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; print("文件%s打开失败" % fname) &amp;nbsp; &amp;nbsp; t.extractall(path=dirs) def Block(IDataSet,filename_split,outPath_file,ImageType,config,metedata): #def Block(IDataSet,ImageType,config,outFileName): &amp;nbsp; &amp;nbsp; # global cols,rows,atcfiles &amp;nbsp; &amp;nbsp; cols = IDataSet.RasterXSize &amp;nbsp; &amp;nbsp; rows = IDataSet.RasterYSize &amp;nbsp; &amp;nbsp; SatelliteID = filename_split[0] &amp;nbsp; &amp;nbsp; SensorID = filename_split[1] &amp;nbsp; &amp;nbsp; Year = filename_split[4][:4] &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; #设置输出波段 &amp;nbsp; &amp;nbsp; Driver = IDataSet.GetDriver() &amp;nbsp; &amp;nbsp; geoTransform1 = IDataSet.GetGeoTransform() &amp;nbsp; &amp;nbsp; ListgeoTransform1 = list(geoTransform1) &amp;nbsp; &amp;nbsp; ListgeoTransform1[5] = -ListgeoTransform1[5] &amp;nbsp; &amp;nbsp; newgeoTransform1 = tuple(ListgeoTransform1) &amp;nbsp; &amp;nbsp; proj1 = IDataSet.GetProjection() &amp;nbsp; &amp;nbsp; #OutRCname = os.path.join(atcfiles,outFileName+".tiff") &amp;nbsp; &amp;nbsp; #OutRCname = 'data/gf2_0503_rf.tif' &amp;nbsp; &amp;nbsp; OutRCname = outPath_file &amp;nbsp; &amp;nbsp; outDataset = Driver.Create(OutRCname,cols,rows,4,gdal.GDT_Int32) &amp;nbsp; &amp;nbsp; #outDataset.SetGeoTransform(newgeoTransform1) &amp;nbsp; &amp;nbsp; outDataset.SetGeoTransform(geoTransform1) &amp;nbsp; &amp;nbsp; outDataset.SetProjection(proj1) &amp;nbsp; &amp;nbsp; #分别读取4个波段 &amp;nbsp; &amp;nbsp; for m in range(1,5): &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; ReadBand = IDataSet.GetRasterBand(m) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; outband = outDataset.GetRasterBand(m) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; outband.SetNoDataValue(-9999) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; #获取对应波段的增益gain和偏移bias &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; Gain,Bias = RadiometricCalibration(m,SatelliteID,SensorID,Year,ImageType,config) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; #获取大气校正系数 &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; AtcCofa, AtcCofb, AtcCofc = AtmosphericCorrection(m,metedata,config,SatelliteID,SensorID) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; nBlockSize = 2048 &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; i = 0 &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; j = 0 &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; b = cols*rows &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; #进度条参数 &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; XBlockcount = math.ceil(cols/nBlockSize) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; YBlockcount = math.ceil(rows/nBlockSize) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; print("第%d波段校正："%m) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; try: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; with tqdm(total=XBlockcount*YBlockcount,iterable='iterable',desc = '第%i波段:'%m,mininterval=10) as pbar: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; # with tqdm(total=XBlockcount*YBlockcount) as pbar: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; # print(pbar) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; while i<rows: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; while j <cols: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; #保存分块大小 &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; nXBK = nBlockSize &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; nYBK = nBlockSize &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; #最后不够分块的区域，有多少读取多少 &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; if i+nBlockSize&amp;gt;rows: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; nYBK = rows - i &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; if j+nBlockSize&amp;gt;cols: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; nXBK=cols - j &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; #分块读取影像 &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; Image = ReadBand.ReadAsArray(j,i,nXBK,nYBK) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; outImage =np.where(Image&amp;gt;0,Image*Gain + Bias,-9999) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; y = np.where(outImage!=-9999,AtcCofa * outImage - AtcCofb,-9999) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; atcImage = np.where(y!=-9999,(y / (1 + y * AtcCofc))*10000,-9999) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; outband.WriteArray(atcImage,j,i) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; j=j+nXBK &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; time.sleep(1) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; pbar.update(1) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; j=0 &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; i=i+nYBK &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; except KeyboardInterrupt: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; pbar.close() &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; raise &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; pbar.close() def RadiometricCalibration(BandId,SatelliteID,SensorID,Year,ImageType,config): &amp;nbsp; &amp;nbsp; # global cols,rows,SatelliteID,SensorID,Year,ImageType,config &amp;nbsp; &amp;nbsp; if SensorID[0:3] == "WFV": &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; Gain_ =config["Parameter"][SatelliteID][SensorID][Year]["gain"][BandId-1] &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; Bias_ =config["Parameter"][SatelliteID][SensorID][Year]["offset"][BandId-1] &amp;nbsp; &amp;nbsp; else: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; Gain_ =config["Parameter"][SatelliteID][SensorID][Year][ImageType]["gain"][BandId-1] &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; Bias_ =config["Parameter"][SatelliteID][SensorID][Year][ImageType]["offset"][BandId-1] &amp;nbsp; &amp;nbsp; return Gain_,Bias_ # 6s大气校正 #def AtmosphericCorrection(BandId,config,SatelliteID,SensorID): def AtmosphericCorrection(BandId,metedata,config,SatelliteID,SensorID): &amp;nbsp; &amp;nbsp; # global metedata,config,SatelliteID,SensorID &amp;nbsp; &amp;nbsp; #读取头文件 &amp;nbsp; &amp;nbsp; dom = xml.dom.minidom.parse(metedata) &amp;nbsp; &amp;nbsp; # 6S模型 &amp;nbsp; &amp;nbsp; s = SixS() &amp;nbsp; &amp;nbsp; # 传感器类型 自定义 &amp;nbsp; &amp;nbsp; s.geometry = Geometry.User() &amp;nbsp; &amp;nbsp; #s.geometry.solar_z = 90-float(dom.getElementsByTagName('SolarZenith')[0].firstChild.data) &amp;nbsp; &amp;nbsp; s.geometry.solar_z = float(dom.getElementsByTagName('SolarZenith')[0].firstChild.data) &amp;nbsp; &amp;nbsp; s.geometry.solar_a = float(dom.getElementsByTagName('SolarAzimuth')[0].firstChild.data) &amp;nbsp; &amp;nbsp; #s.geometry.solar_a = 147.636 &amp;nbsp; &amp;nbsp; #s.geometry.solar_z = 27.5293 &amp;nbsp; &amp;nbsp; #s.geometry.solar_a = 146.242 &amp;nbsp; &amp;nbsp; #s.geometry.solar_z = 27.8728 &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; #s.geometry.view_z = float(dom.getElementsByTagName('SatelliteZenith')[0].firstChild.data) &amp;nbsp; &amp;nbsp; #s.geometry.view_a = float(dom.getElementsByTagName('SatelliteAzimuth')[0].firstChild.data) &amp;nbsp; &amp;nbsp; s.geometry.view_z = 0 &amp;nbsp; &amp;nbsp; s.geometry.view_a = 0 &amp;nbsp; &amp;nbsp; # 日期 &amp;nbsp; &amp;nbsp; DateTimeparm = dom.getElementsByTagName('CenterTime')[0].firstChild.data &amp;nbsp; &amp;nbsp; DateTime = DateTimeparm.split(' ') &amp;nbsp; &amp;nbsp; Date = DateTime[0].split('-') &amp;nbsp; &amp;nbsp; s.geometry.month = int(Date[1]) &amp;nbsp; &amp;nbsp; s.geometry.day = int(Date[2]) &amp;nbsp; &amp;nbsp; #s.geometry.month = 5 &amp;nbsp; &amp;nbsp; #s.geometry.day = 3 &amp;nbsp; &amp;nbsp; # print(s.geometry) &amp;nbsp; &amp;nbsp; # 中心经纬度 &amp;nbsp; &amp;nbsp; TopLeftLat = float(dom.getElementsByTagName('TopLeftLatitude')[0].firstChild.data) &amp;nbsp; &amp;nbsp; TopLeftLon = float(dom.getElementsByTagName('TopLeftLongitude')[0].firstChild.data) &amp;nbsp; &amp;nbsp; TopRightLat = float(dom.getElementsByTagName('TopRightLatitude')[0].firstChild.data) &amp;nbsp; &amp;nbsp; TopRightLon = float(dom.getElementsByTagName('TopRightLongitude')[0].firstChild.data) &amp;nbsp; &amp;nbsp; BottomRightLat = float(dom.getElementsByTagName('BottomRightLatitude')[0].firstChild.data) &amp;nbsp; &amp;nbsp; BottomRightLon = float(dom.getElementsByTagName('BottomRightLongitude')[0].firstChild.data) &amp;nbsp; &amp;nbsp; BottomLeftLat = float(dom.getElementsByTagName('BottomLeftLatitude')[0].firstChild.data) &amp;nbsp; &amp;nbsp; BottomLeftLon = float(dom.getElementsByTagName('BottomLeftLongitude')[0].firstChild.data) &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; #TopLeftLat = 39.848448 &amp;nbsp; &amp;nbsp; #TopLeftLon = 115.244064 &amp;nbsp; &amp;nbsp; #TopRightLat = 39.848448 &amp;nbsp; &amp;nbsp; #TopRightLon = 115.569216 &amp;nbsp; &amp;nbsp; #BottomRightLat = 39.604288 &amp;nbsp; &amp;nbsp; #BottomRightLon = 115.569216 &amp;nbsp; &amp;nbsp; #BottomLeftLat = 39.604288 &amp;nbsp; &amp;nbsp; #BottomLeftLon = 115.244064 &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; # TopLeftLat = 40.15664 &amp;nbsp; &amp;nbsp; # TopLeftLon = 116.033536 &amp;nbsp; &amp;nbsp; # TopRightLat = 40.15664 &amp;nbsp; &amp;nbsp; # TopRightLon = 116.356896 &amp;nbsp; &amp;nbsp; # BottomRightLat = 39.914944 &amp;nbsp; &amp;nbsp; # BottomRightLon = 116.356896 &amp;nbsp; &amp;nbsp; # BottomLeftLat = 39.914944 &amp;nbsp; &amp;nbsp; # BottomLeftLon = 116.033536 &amp;nbsp; &amp;nbsp; ImageCenterLat = (TopLeftLat + TopRightLat + BottomRightLat + BottomLeftLat) / 4 &amp;nbsp; &amp;nbsp; # 大气模式类型 &amp;nbsp; &amp;nbsp; if ImageCenterLat &amp;gt; -15 and ImageCenterLat < 15: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical) &amp;nbsp; &amp;nbsp; if ImageCenterLat &amp;gt; 15 and ImageCenterLat < 45: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; if s.geometry.month &amp;gt; 4 and s.geometry.month < 9: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; else: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeWinter) &amp;nbsp; &amp;nbsp; if ImageCenterLat &amp;gt; 45 and ImageCenterLat < 60: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; if s.geometry.month &amp;gt; 4 and s.geometry.month < 9: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticSummer) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; else: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticWinter) &amp;nbsp; &amp;nbsp; # 气溶胶类型大陆 &amp;nbsp; &amp;nbsp; s.aero_profile = AtmosProfile.PredefinedType(AeroProfile.Continental) &amp;nbsp; &amp;nbsp; # 下垫面类型 &amp;nbsp; &amp;nbsp; s.ground_reflectance = GroundReflectance.HomogeneousLambertian(0.36) &amp;nbsp; &amp;nbsp; #s.ground_reflectance = GroundReflectance.HomogeneousLambertian(GroundReflectance.ClearWater) &amp;nbsp; &amp;nbsp; # 550nm气溶胶光学厚度,对应能见度为40km &amp;nbsp; &amp;nbsp; s.aot550 = 0.14497 &amp;nbsp; &amp;nbsp; #s.aot550 = 40 &amp;nbsp; &amp;nbsp; # 通过研究去区的范围去求DEM高度。 &amp;nbsp; &amp;nbsp; pointUL = dict() &amp;nbsp; &amp;nbsp; pointDR = dict() &amp;nbsp; &amp;nbsp; pointUL["lat"] = max(TopLeftLat,TopRightLat,BottomRightLat,BottomLeftLat) &amp;nbsp; &amp;nbsp; pointUL["lon"] = min(TopLeftLon,TopRightLon,BottomRightLon,BottomLeftLon) &amp;nbsp; &amp;nbsp; pointDR["lat"] = min(TopLeftLat,TopRightLat,BottomRightLat,BottomLeftLat) &amp;nbsp; &amp;nbsp; pointDR["lon"] = max(TopLeftLon,TopRightLon,BottomRightLon,BottomLeftLon) &amp;nbsp; &amp;nbsp; #meanDEM = (MeanDEM(pointUL, pointDR)) * 0.001 &amp;nbsp; &amp;nbsp; meanDEM = 0.0 &amp;nbsp; &amp;nbsp; # 研究区海拔、卫星传感器轨道高度 &amp;nbsp; &amp;nbsp; s.altitudes = Altitudes() &amp;nbsp; &amp;nbsp; s.altitudes.set_target_custom_altitude(meanDEM) &amp;nbsp; &amp;nbsp; s.altitudes.set_sensor_satellite_level() &amp;nbsp; &amp;nbsp; # 校正波段（根据波段名称） &amp;nbsp; &amp;nbsp; if BandId == 1: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; SRFband = config["Parameter"][SatelliteID][SensorID]["SRF"]["1"] &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; s.wavelength = Wavelength(0.450,0.520,SRFband) &amp;nbsp; &amp;nbsp; elif BandId == 2: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; SRFband = config["Parameter"][SatelliteID][SensorID]["SRF"]["2"] &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; s.wavelength = Wavelength(0.520,0.590,SRFband) &amp;nbsp; &amp;nbsp; elif BandId == 3: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; SRFband = config["Parameter"][SatelliteID][SensorID]["SRF"]["3"] &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; s.wavelength = Wavelength(0.630,0.690,SRFband) &amp;nbsp; &amp;nbsp; elif BandId == 4: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; SRFband = config["Parameter"][SatelliteID][SensorID]["SRF"]["4"] &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; s.wavelength = Wavelength(0.770,0.890,SRFband) &amp;nbsp; &amp;nbsp; s.atmos_corr = AtmosCorr.AtmosCorrLambertianFromReflectance(-0.1) &amp;nbsp; &amp;nbsp; # 运行6s大气模型 &amp;nbsp; &amp;nbsp; s.run() &amp;nbsp; &amp;nbsp; xa = s.outputs.coef_xa &amp;nbsp; &amp;nbsp; xb = s.outputs.coef_xb &amp;nbsp; &amp;nbsp; xc = s.outputs.coef_xc &amp;nbsp; &amp;nbsp; # x = s.outputs.values &amp;nbsp; &amp;nbsp; return (xa, xb, xc) if __name__ == '__main__': &amp;nbsp; &amp;nbsp; a=time.time() &amp;nbsp; &amp;nbsp; script_path = os.path.split(os.path.realpath(__file__))[0] &amp;nbsp; &amp;nbsp; #读取辐射校正和大气校正所需参数:增益、偏移和光谱响应函数 &amp;nbsp; &amp;nbsp; config_file = os.path.join(script_path,"RadiometricCorrectionParameter.json") &amp;nbsp; &amp;nbsp; config = json.load(open(config_file)) &amp;nbsp; &amp;nbsp; #输入数据路径 &amp;nbsp; &amp;nbsp; #InputFilePath = parse_arguments(sys.argv[1:]).Input_dir &amp;nbsp; &amp;nbsp; #OutputFilePath = parse_arguments(sys.argv[2:]).Output_dir &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; #头文件文件夹 &amp;nbsp; &amp;nbsp; xml_path = r'F:\2-Project-data\BeijingEco\xml' &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; #遍历文件夹 &amp;nbsp; &amp;nbsp; input_path = r'F:\2-Project-data\BeijingEco\output_proj' &amp;nbsp; &amp;nbsp; #input_path = r'F:\2-Project-data\BeijingEco\test_ac' &amp;nbsp; &amp;nbsp; files = os.listdir(input_path) &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; output_path = r'F:\2-Project-data\BeijingEco\atmospheric_correction' &amp;nbsp; &amp;nbsp; #file = 'GF2_PMS2_E116.2_N40.0_20210503_L1A0005629120-MSS2_ORTHO_MS.tif' &amp;nbsp; &amp;nbsp; for file in files: &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; #file:'GF2_PMS1_E115.6_N39.9_20200602_L1A0004842146-MSS1_ORTHO_MS_proj.tif'，含有tif &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; if os.path.splitext(file)[1].upper() == ".TIF": &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; print(file) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; tiffFile = input_path + os.sep + file &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; IDataSet = gdal.Open(tiffFile) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; ImageType = 'MSS' &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; filename = os.path.splitext(file)[0]&amp;nbsp; #'GF2_PMS1_E115.6_N39.9_20200602_L1A0004842146-MSS1_ORTHO_MS_proj' &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; fileType = filename[0:2] &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; filename_split = filename.split("_") &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; outFileName = filename[:-4]+'RF.tif' &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; outPath_file = os.path.join(output_path, outFileName) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; GFType = filename_split[1][:3] #GF2 &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; #读取头文件 &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; productId = filename_split[5] &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; xml_file = search_file(xml_path,productId) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; metedata = glob.glob(os.path.join(xml_path,xml_file))[0] &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; #metedata = glob.glob(os.path.join(outname,"*MSS*.xml"))[0] &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; #outFileName = os.path.splitext(file)[0] &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; #outFileName = output_path+os.sep+ &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; #Block(IDataSet,filename_split,outPath_file,ImageType,config,metedata) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; Block(IDataSet,filename_split,outPath_file,ImageType,config,metedata) &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp; b=time.time() &amp;nbsp; &amp;nbsp; print("总时间:",b-a) &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp; &amp;nbsp;&amp;nbsp; &amp;nbsp;

------------------&amp;nbsp;原始邮件&amp;nbsp;------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***&amp;gt;; 发送时间:&amp;nbsp;2021年12月22日(星期三) 下午4:24 ***@***.***&amp;gt;; ***@***.******@***.***&amp;gt;; 主题:&amp;nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 抱歉啊，还是看不见，看来屏蔽了，你能看到我的扣扣吗1310118141，如果可以麻烦你可以发到我邮箱吗万分感谢

------------------&amp;amp;nbsp;原始邮件&amp;amp;nbsp;------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***&amp;amp;gt;; 发送时间:&amp;amp;nbsp;2021年12月22日(星期三) 下午4:22 ***@***.***&amp;amp;gt;; ***@***.******@***.***&amp;amp;gt;; 主题:&amp;amp;nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 我拿到的影像是经过正射校正后的，不太清楚没有正射校正的影像是什么情况。这次附件里可以看到了吗？ &amp;amp;amp;nbsp;

------------------&amp;amp;amp;nbsp;原始邮件&amp;amp;amp;nbsp;------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***&amp;amp;amp;gt;; 发送时间:&amp;amp;amp;nbsp;2021年12月22日(星期三) 下午4:15 ***@***.***&amp;amp;amp;gt;; ***@***.******@***.***&amp;amp;amp;gt;; 主题:&amp;amp;amp;nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 那请问你用高分原始影像跑的话，在没有正射校正的情况下有坐标吗

------------------&amp;amp;amp;amp;nbsp;原始邮件&amp;amp;amp;amp;nbsp;------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***&amp;amp;amp;amp;gt;; 发送时间:&amp;amp;amp;amp;nbsp;2021年12月22日(星期三) 下午4:12 ***@***.***&amp;amp;amp;amp;gt;; ***@***.******@***.***&amp;amp;amp;amp;gt;; 主题:&amp;amp;amp;amp;nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 你好，这个问题没有解决。 关于影像坐标的问题，我先给正射校正后的影像加上了投影，然后再拿带有坐标的影像进行的大气校正。 附件是我的代码，供您参考。 &amp;amp;amp;amp;amp;nbsp;

------------------&amp;amp;amp;amp;amp;nbsp;原始邮件&amp;amp;amp;amp;amp;nbsp;------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***&amp;amp;amp;amp;amp;gt;; 发送时间:&amp;amp;amp;amp;amp;nbsp;2021年12月22日(星期三) 下午4:05 ***@***.***&amp;amp;amp;amp;amp;gt;; ***@***.******@***.***&amp;amp;amp;amp;amp;gt;; 主题:&amp;amp;amp;amp;amp;nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 你好，请问你的问题解决了吗，能否将你的代码发送给我，我下载这个作者的代码运行的影像没有坐标 — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***&amp;amp;amp;amp;amp;gt; — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.Message ID: ***@***.***&amp;amp;amp;amp;gt; — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***&amp;amp;amp;gt; — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.Message ID: ***@***.***&amp;amp;gt; — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***&amp;gt; — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.Message ID: ***@***.***&gt;

[alive9999]

我加你QQ吧 &nbsp;

…

------------------&nbsp;原始邮件&nbsp;------------------ 发件人: "Zhaoguanhua/AtmosphericCorrection" ***@***.***&gt;; 发送时间:&nbsp;2021年12月22日(星期三) 下午4:31 ***@***.***&gt;; ***@***.******@***.***&gt;; 主题:&nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 抱歉，是我的号被屏蔽了吗？我这边看你复制的内容很大部分都被屏蔽了

------------------&amp;nbsp;原始邮件&amp;nbsp;------------------ 发件人:"Zhaoguanhua/AtmosphericCorrection"***@***.***&amp;gt;; 发送时间:&amp;nbsp;2021年12月22日(星期三) 下午4:26 ***@***.***&amp;gt;; ***@***.******@***.***&amp;gt;; 主题:&amp;nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 我直接粘贴到邮件里吧~ #! usr/bin/env python # -*- coding:utf-8 -*- # created by zhaoguanhua 2017/09/04 import glob import os import sys import tarfile&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #解压缩 import json import numpy as np import gdal import pdb import math import time import xml.dom.minidom&amp;amp;nbsp; &amp;amp;nbsp; #读取xml格式的影像头文件 from tqdm import tqdm&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp;#进度条 from Py6S import * import argparse from base import MeanDEM import shutil def parse_arguments(argv): &amp;amp;nbsp; &amp;amp;nbsp; parser = argparse.ArgumentParser() &amp;amp;nbsp; &amp;amp;nbsp; parser.add_argument('--Input_dir',type=str,help='Input dir',default=None) &amp;amp;nbsp; &amp;amp;nbsp; parser.add_argument('--Output_dir',type=str,help='Output dir',default=None) &amp;amp;nbsp; &amp;amp;nbsp; return parser.parse_args(argv) def search_file(path,str): &amp;amp;nbsp; &amp;amp;nbsp; for x in os.listdir(path): &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; next_path = os.path.join(path,x) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; if os.path.isfile(next_path): &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; if str in x: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; return x # 解压缩原始文件 def untar(fname, dirs): &amp;amp;nbsp; &amp;amp;nbsp; print("文件路径",fname) &amp;amp;nbsp; &amp;amp;nbsp; try: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; t = tarfile.open(fname) &amp;amp;nbsp; &amp;amp;nbsp; except Exception as e: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; print("文件%s打开失败" % fname) &amp;amp;nbsp; &amp;amp;nbsp; t.extractall(path=dirs) def Block(IDataSet,filename_split,outPath_file,ImageType,config,metedata): #def Block(IDataSet,ImageType,config,outFileName): &amp;amp;nbsp; &amp;amp;nbsp; # global cols,rows,atcfiles &amp;amp;nbsp; &amp;amp;nbsp; cols = IDataSet.RasterXSize &amp;amp;nbsp; &amp;amp;nbsp; rows = IDataSet.RasterYSize &amp;amp;nbsp; &amp;amp;nbsp; SatelliteID = filename_split[0] &amp;amp;nbsp; &amp;amp;nbsp; SensorID = filename_split[1] &amp;amp;nbsp; &amp;amp;nbsp; Year = filename_split[4][:4] &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #设置输出波段 &amp;amp;nbsp; &amp;amp;nbsp; Driver = IDataSet.GetDriver() &amp;amp;nbsp; &amp;amp;nbsp; geoTransform1 = IDataSet.GetGeoTransform() &amp;amp;nbsp; &amp;amp;nbsp; ListgeoTransform1 = list(geoTransform1) &amp;amp;nbsp; &amp;amp;nbsp; ListgeoTransform1[5] = -ListgeoTransform1[5] &amp;amp;nbsp; &amp;amp;nbsp; newgeoTransform1 = tuple(ListgeoTransform1) &amp;amp;nbsp; &amp;amp;nbsp; proj1 = IDataSet.GetProjection() &amp;amp;nbsp; &amp;amp;nbsp; #OutRCname = os.path.join(atcfiles,outFileName+".tiff") &amp;amp;nbsp; &amp;amp;nbsp; #OutRCname = 'data/gf2_0503_rf.tif' &amp;amp;nbsp; &amp;amp;nbsp; OutRCname = outPath_file &amp;amp;nbsp; &amp;amp;nbsp; outDataset = Driver.Create(OutRCname,cols,rows,4,gdal.GDT_Int32) &amp;amp;nbsp; &amp;amp;nbsp; #outDataset.SetGeoTransform(newgeoTransform1) &amp;amp;nbsp; &amp;amp;nbsp; outDataset.SetGeoTransform(geoTransform1) &amp;amp;nbsp; &amp;amp;nbsp; outDataset.SetProjection(proj1) &amp;amp;nbsp; &amp;amp;nbsp; #分别读取4个波段 &amp;amp;nbsp; &amp;amp;nbsp; for m in range(1,5): &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; ReadBand = IDataSet.GetRasterBand(m) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; outband = outDataset.GetRasterBand(m) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; outband.SetNoDataValue(-9999) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #获取对应波段的增益gain和偏移bias &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; Gain,Bias = RadiometricCalibration(m,SatelliteID,SensorID,Year,ImageType,config) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #获取大气校正系数 &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; AtcCofa, AtcCofb, AtcCofc = AtmosphericCorrection(m,metedata,config,SatelliteID,SensorID) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; nBlockSize = 2048 &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; i = 0 &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; j = 0 &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; b = cols*rows &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #进度条参数 &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; XBlockcount = math.ceil(cols/nBlockSize) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; YBlockcount = math.ceil(rows/nBlockSize) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; print("第%d波段校正："%m) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; try: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; with tqdm(total=XBlockcount*YBlockcount,iterable='iterable',desc = '第%i波段:'%m,mininterval=10) as pbar: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; # with tqdm(total=XBlockcount*YBlockcount) as pbar: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; # print(pbar) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; while i<rows: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; while j <cols: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #保存分块大小 &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; nXBK = nBlockSize &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; nYBK = nBlockSize &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #最后不够分块的区域，有多少读取多少 &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; if i+nBlockSize&amp;amp;gt;rows: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; nYBK = rows - i &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; if j+nBlockSize&amp;amp;gt;cols: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; nXBK=cols - j &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #分块读取影像 &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; Image = ReadBand.ReadAsArray(j,i,nXBK,nYBK) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; outImage =np.where(Image&amp;amp;gt;0,Image*Gain + Bias,-9999) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; y = np.where(outImage!=-9999,AtcCofa * outImage - AtcCofb,-9999) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; atcImage = np.where(y!=-9999,(y / (1 + y * AtcCofc))*10000,-9999) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; outband.WriteArray(atcImage,j,i) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; j=j+nXBK &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; time.sleep(1) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; pbar.update(1) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; j=0 &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; i=i+nYBK &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; except KeyboardInterrupt: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; pbar.close() &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; raise &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; pbar.close() def RadiometricCalibration(BandId,SatelliteID,SensorID,Year,ImageType,config): &amp;amp;nbsp; &amp;amp;nbsp; # global cols,rows,SatelliteID,SensorID,Year,ImageType,config &amp;amp;nbsp; &amp;amp;nbsp; if SensorID[0:3] == "WFV": &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; Gain_ =config["Parameter"][SatelliteID][SensorID][Year]["gain"][BandId-1] &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; Bias_ =config["Parameter"][SatelliteID][SensorID][Year]["offset"][BandId-1] &amp;amp;nbsp; &amp;amp;nbsp; else: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; Gain_ =config["Parameter"][SatelliteID][SensorID][Year][ImageType]["gain"][BandId-1] &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; Bias_ =config["Parameter"][SatelliteID][SensorID][Year][ImageType]["offset"][BandId-1] &amp;amp;nbsp; &amp;amp;nbsp; return Gain_,Bias_ # 6s大气校正 #def AtmosphericCorrection(BandId,config,SatelliteID,SensorID): def AtmosphericCorrection(BandId,metedata,config,SatelliteID,SensorID): &amp;amp;nbsp; &amp;amp;nbsp; # global metedata,config,SatelliteID,SensorID &amp;amp;nbsp; &amp;amp;nbsp; #读取头文件 &amp;amp;nbsp; &amp;amp;nbsp; dom = xml.dom.minidom.parse(metedata) &amp;amp;nbsp; &amp;amp;nbsp; # 6S模型 &amp;amp;nbsp; &amp;amp;nbsp; s = SixS() &amp;amp;nbsp; &amp;amp;nbsp; # 传感器类型 自定义 &amp;amp;nbsp; &amp;amp;nbsp; s.geometry = Geometry.User() &amp;amp;nbsp; &amp;amp;nbsp; #s.geometry.solar_z = 90-float(dom.getElementsByTagName('SolarZenith')[0].firstChild.data) &amp;amp;nbsp; &amp;amp;nbsp; s.geometry.solar_z = float(dom.getElementsByTagName('SolarZenith')[0].firstChild.data) &amp;amp;nbsp; &amp;amp;nbsp; s.geometry.solar_a = float(dom.getElementsByTagName('SolarAzimuth')[0].firstChild.data) &amp;amp;nbsp; &amp;amp;nbsp; #s.geometry.solar_a = 147.636 &amp;amp;nbsp; &amp;amp;nbsp; #s.geometry.solar_z = 27.5293 &amp;amp;nbsp; &amp;amp;nbsp; #s.geometry.solar_a = 146.242 &amp;amp;nbsp; &amp;amp;nbsp; #s.geometry.solar_z = 27.8728 &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #s.geometry.view_z = float(dom.getElementsByTagName('SatelliteZenith')[0].firstChild.data) &amp;amp;nbsp; &amp;amp;nbsp; #s.geometry.view_a = float(dom.getElementsByTagName('SatelliteAzimuth')[0].firstChild.data) &amp;amp;nbsp; &amp;amp;nbsp; s.geometry.view_z = 0 &amp;amp;nbsp; &amp;amp;nbsp; s.geometry.view_a = 0 &amp;amp;nbsp; &amp;amp;nbsp; # 日期 &amp;amp;nbsp; &amp;amp;nbsp; DateTimeparm = dom.getElementsByTagName('CenterTime')[0].firstChild.data &amp;amp;nbsp; &amp;amp;nbsp; DateTime = DateTimeparm.split(' ') &amp;amp;nbsp; &amp;amp;nbsp; Date = DateTime[0].split('-') &amp;amp;nbsp; &amp;amp;nbsp; s.geometry.month = int(Date[1]) &amp;amp;nbsp; &amp;amp;nbsp; s.geometry.day = int(Date[2]) &amp;amp;nbsp; &amp;amp;nbsp; #s.geometry.month = 5 &amp;amp;nbsp; &amp;amp;nbsp; #s.geometry.day = 3 &amp;amp;nbsp; &amp;amp;nbsp; # print(s.geometry) &amp;amp;nbsp; &amp;amp;nbsp; # 中心经纬度 &amp;amp;nbsp; &amp;amp;nbsp; TopLeftLat = float(dom.getElementsByTagName('TopLeftLatitude')[0].firstChild.data) &amp;amp;nbsp; &amp;amp;nbsp; TopLeftLon = float(dom.getElementsByTagName('TopLeftLongitude')[0].firstChild.data) &amp;amp;nbsp; &amp;amp;nbsp; TopRightLat = float(dom.getElementsByTagName('TopRightLatitude')[0].firstChild.data) &amp;amp;nbsp; &amp;amp;nbsp; TopRightLon = float(dom.getElementsByTagName('TopRightLongitude')[0].firstChild.data) &amp;amp;nbsp; &amp;amp;nbsp; BottomRightLat = float(dom.getElementsByTagName('BottomRightLatitude')[0].firstChild.data) &amp;amp;nbsp; &amp;amp;nbsp; BottomRightLon = float(dom.getElementsByTagName('BottomRightLongitude')[0].firstChild.data) &amp;amp;nbsp; &amp;amp;nbsp; BottomLeftLat = float(dom.getElementsByTagName('BottomLeftLatitude')[0].firstChild.data) &amp;amp;nbsp; &amp;amp;nbsp; BottomLeftLon = float(dom.getElementsByTagName('BottomLeftLongitude')[0].firstChild.data) &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #TopLeftLat = 39.848448 &amp;amp;nbsp; &amp;amp;nbsp; #TopLeftLon = 115.244064 &amp;amp;nbsp; &amp;amp;nbsp; #TopRightLat = 39.848448 &amp;amp;nbsp; &amp;amp;nbsp; #TopRightLon = 115.569216 &amp;amp;nbsp; &amp;amp;nbsp; #BottomRightLat = 39.604288 &amp;amp;nbsp; &amp;amp;nbsp; #BottomRightLon = 115.569216 &amp;amp;nbsp; &amp;amp;nbsp; #BottomLeftLat = 39.604288 &amp;amp;nbsp; &amp;amp;nbsp; #BottomLeftLon = 115.244064 &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; # TopLeftLat = 40.15664 &amp;amp;nbsp; &amp;amp;nbsp; # TopLeftLon = 116.033536 &amp;amp;nbsp; &amp;amp;nbsp; # TopRightLat = 40.15664 &amp;amp;nbsp; &amp;amp;nbsp; # TopRightLon = 116.356896 &amp;amp;nbsp; &amp;amp;nbsp; # BottomRightLat = 39.914944 &amp;amp;nbsp; &amp;amp;nbsp; # BottomRightLon = 116.356896 &amp;amp;nbsp; &amp;amp;nbsp; # BottomLeftLat = 39.914944 &amp;amp;nbsp; &amp;amp;nbsp; # BottomLeftLon = 116.033536 &amp;amp;nbsp; &amp;amp;nbsp; ImageCenterLat = (TopLeftLat + TopRightLat + BottomRightLat + BottomLeftLat) / 4 &amp;amp;nbsp; &amp;amp;nbsp; # 大气模式类型 &amp;amp;nbsp; &amp;amp;nbsp; if ImageCenterLat &amp;amp;gt; -15 and ImageCenterLat < 15: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical) &amp;amp;nbsp; &amp;amp;nbsp; if ImageCenterLat &amp;amp;gt; 15 and ImageCenterLat < 45: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; if s.geometry.month &amp;amp;gt; 4 and s.geometry.month < 9: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; else: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeWinter) &amp;amp;nbsp; &amp;amp;nbsp; if ImageCenterLat &amp;amp;gt; 45 and ImageCenterLat < 60: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; if s.geometry.month &amp;amp;gt; 4 and s.geometry.month < 9: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticSummer) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; else: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticWinter) &amp;amp;nbsp; &amp;amp;nbsp; # 气溶胶类型大陆 &amp;amp;nbsp; &amp;amp;nbsp; s.aero_profile = AtmosProfile.PredefinedType(AeroProfile.Continental) &amp;amp;nbsp; &amp;amp;nbsp; # 下垫面类型 &amp;amp;nbsp; &amp;amp;nbsp; s.ground_reflectance = GroundReflectance.HomogeneousLambertian(0.36) &amp;amp;nbsp; &amp;amp;nbsp; #s.ground_reflectance = GroundReflectance.HomogeneousLambertian(GroundReflectance.ClearWater) &amp;amp;nbsp; &amp;amp;nbsp; # 550nm气溶胶光学厚度,对应能见度为40km &amp;amp;nbsp; &amp;amp;nbsp; s.aot550 = 0.14497 &amp;amp;nbsp; &amp;amp;nbsp; #s.aot550 = 40 &amp;amp;nbsp; &amp;amp;nbsp; # 通过研究去区的范围去求DEM高度。 &amp;amp;nbsp; &amp;amp;nbsp; pointUL = dict() &amp;amp;nbsp; &amp;amp;nbsp; pointDR = dict() &amp;amp;nbsp; &amp;amp;nbsp; pointUL["lat"] = max(TopLeftLat,TopRightLat,BottomRightLat,BottomLeftLat) &amp;amp;nbsp; &amp;amp;nbsp; pointUL["lon"] = min(TopLeftLon,TopRightLon,BottomRightLon,BottomLeftLon) &amp;amp;nbsp; &amp;amp;nbsp; pointDR["lat"] = min(TopLeftLat,TopRightLat,BottomRightLat,BottomLeftLat) &amp;amp;nbsp; &amp;amp;nbsp; pointDR["lon"] = max(TopLeftLon,TopRightLon,BottomRightLon,BottomLeftLon) &amp;amp;nbsp; &amp;amp;nbsp; #meanDEM = (MeanDEM(pointUL, pointDR)) * 0.001 &amp;amp;nbsp; &amp;amp;nbsp; meanDEM = 0.0 &amp;amp;nbsp; &amp;amp;nbsp; # 研究区海拔、卫星传感器轨道高度 &amp;amp;nbsp; &amp;amp;nbsp; s.altitudes = Altitudes() &amp;amp;nbsp; &amp;amp;nbsp; s.altitudes.set_target_custom_altitude(meanDEM) &amp;amp;nbsp; &amp;amp;nbsp; s.altitudes.set_sensor_satellite_level() &amp;amp;nbsp; &amp;amp;nbsp; # 校正波段（根据波段名称） &amp;amp;nbsp; &amp;amp;nbsp; if BandId == 1: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; SRFband = config["Parameter"][SatelliteID][SensorID]["SRF"]["1"] &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; s.wavelength = Wavelength(0.450,0.520,SRFband) &amp;amp;nbsp; &amp;amp;nbsp; elif BandId == 2: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; SRFband = config["Parameter"][SatelliteID][SensorID]["SRF"]["2"] &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; s.wavelength = Wavelength(0.520,0.590,SRFband) &amp;amp;nbsp; &amp;amp;nbsp; elif BandId == 3: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; SRFband = config["Parameter"][SatelliteID][SensorID]["SRF"]["3"] &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; s.wavelength = Wavelength(0.630,0.690,SRFband) &amp;amp;nbsp; &amp;amp;nbsp; elif BandId == 4: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; SRFband = config["Parameter"][SatelliteID][SensorID]["SRF"]["4"] &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; s.wavelength = Wavelength(0.770,0.890,SRFband) &amp;amp;nbsp; &amp;amp;nbsp; s.atmos_corr = AtmosCorr.AtmosCorrLambertianFromReflectance(-0.1) &amp;amp;nbsp; &amp;amp;nbsp; # 运行6s大气模型 &amp;amp;nbsp; &amp;amp;nbsp; s.run() &amp;amp;nbsp; &amp;amp;nbsp; xa = s.outputs.coef_xa &amp;amp;nbsp; &amp;amp;nbsp; xb = s.outputs.coef_xb &amp;amp;nbsp; &amp;amp;nbsp; xc = s.outputs.coef_xc &amp;amp;nbsp; &amp;amp;nbsp; # x = s.outputs.values &amp;amp;nbsp; &amp;amp;nbsp; return (xa, xb, xc) if __name__ == '__main__': &amp;amp;nbsp; &amp;amp;nbsp; a=time.time() &amp;amp;nbsp; &amp;amp;nbsp; script_path = os.path.split(os.path.realpath(__file__))[0] &amp;amp;nbsp; &amp;amp;nbsp; #读取辐射校正和大气校正所需参数:增益、偏移和光谱响应函数 &amp;amp;nbsp; &amp;amp;nbsp; config_file = os.path.join(script_path,"RadiometricCorrectionParameter.json") &amp;amp;nbsp; &amp;amp;nbsp; config = json.load(open(config_file)) &amp;amp;nbsp; &amp;amp;nbsp; #输入数据路径 &amp;amp;nbsp; &amp;amp;nbsp; #InputFilePath = parse_arguments(sys.argv[1:]).Input_dir &amp;amp;nbsp; &amp;amp;nbsp; #OutputFilePath = parse_arguments(sys.argv[2:]).Output_dir &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #头文件文件夹 &amp;amp;nbsp; &amp;amp;nbsp; xml_path = r'F:\2-Project-data\BeijingEco\xml' &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #遍历文件夹 &amp;amp;nbsp; &amp;amp;nbsp; input_path = r'F:\2-Project-data\BeijingEco\output_proj' &amp;amp;nbsp; &amp;amp;nbsp; #input_path = r'F:\2-Project-data\BeijingEco\test_ac' &amp;amp;nbsp; &amp;amp;nbsp; files = os.listdir(input_path) &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; output_path = r'F:\2-Project-data\BeijingEco\atmospheric_correction' &amp;amp;nbsp; &amp;amp;nbsp; #file = 'GF2_PMS2_E116.2_N40.0_20210503_L1A0005629120-MSS2_ORTHO_MS.tif' &amp;amp;nbsp; &amp;amp;nbsp; for file in files: &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #file:'GF2_PMS1_E115.6_N39.9_20200602_L1A0004842146-MSS1_ORTHO_MS_proj.tif'，含有tif &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; if os.path.splitext(file)[1].upper() == ".TIF": &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; print(file) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; tiffFile = input_path + os.sep + file &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; IDataSet = gdal.Open(tiffFile) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; ImageType = 'MSS' &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; filename = os.path.splitext(file)[0]&amp;amp;nbsp; #'GF2_PMS1_E115.6_N39.9_20200602_L1A0004842146-MSS1_ORTHO_MS_proj' &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; fileType = filename[0:2] &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; filename_split = filename.split("_") &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; outFileName = filename[:-4]+'RF.tif' &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; outPath_file = os.path.join(output_path, outFileName) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; GFType = filename_split[1][:3] #GF2 &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #读取头文件 &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; productId = filename_split[5] &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; xml_file = search_file(xml_path,productId) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; metedata = glob.glob(os.path.join(xml_path,xml_file))[0] &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #metedata = glob.glob(os.path.join(outname,"*MSS*.xml"))[0] &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #outFileName = os.path.splitext(file)[0] &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #outFileName = output_path+os.sep+ &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; #Block(IDataSet,filename_split,outPath_file,ImageType,config,metedata) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; Block(IDataSet,filename_split,outPath_file,ImageType,config,metedata) &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp; b=time.time() &amp;amp;nbsp; &amp;amp;nbsp; print("总时间:",b-a) &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp; &amp;amp;nbsp;&amp;amp;nbsp; &amp;amp;nbsp;

------------------&amp;amp;nbsp;原始邮件&amp;amp;nbsp;------------------ 发件人:"Zhaoguanhua/AtmosphericCorrection"***@***.***&amp;amp;gt;; 发送时间:&amp;amp;nbsp;2021年12月22日(星期三) 下午4:24 ***@***.***&amp;amp;gt;; ***@***.******@***.***&amp;amp;gt;; 主题:&amp;amp;nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 抱歉啊，还是看不见，看来屏蔽了，你能看到我的扣扣吗1310118141，如果可以麻烦你可以发到我邮箱吗万分感谢

------------------&amp;amp;amp;nbsp;原始邮件&amp;amp;amp;nbsp;------------------ 发件人:"Zhaoguanhua/AtmosphericCorrection"***@***.***&amp;amp;amp;gt;; 发送时间:&amp;amp;amp;nbsp;2021年12月22日(星期三) 下午4:22 ***@***.***&amp;amp;amp;gt;; ***@***.******@***.***&amp;amp;amp;gt;; 主题:&amp;amp;amp;nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 我拿到的影像是经过正射校正后的，不太清楚没有正射校正的影像是什么情况。这次附件里可以看到了吗？ &amp;amp;amp;amp;nbsp;

------------------&amp;amp;amp;amp;nbsp;原始邮件&amp;amp;amp;amp;nbsp;------------------ 发件人:"Zhaoguanhua/AtmosphericCorrection"***@***.***&amp;amp;amp;amp;gt;; 发送时间:&amp;amp;amp;amp;nbsp;2021年12月22日(星期三) 下午4:15 ***@***.***&amp;amp;amp;amp;gt;; ***@***.******@***.***&amp;amp;amp;amp;gt;; 主题:&amp;amp;amp;amp;nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 那请问你用高分原始影像跑的话，在没有正射校正的情况下有坐标吗

------------------&amp;amp;amp;amp;amp;nbsp;原始邮件&amp;amp;amp;amp;amp;nbsp;------------------ 发件人:"Zhaoguanhua/AtmosphericCorrection"***@***.***&amp;amp;amp;amp;amp;gt;; 发送时间:&amp;amp;amp;amp;amp;nbsp;2021年12月22日(星期三) 下午4:12 ***@***.***&amp;amp;amp;amp;amp;gt;; ***@***.******@***.***&amp;amp;amp;amp;amp;gt;; 主题:&amp;amp;amp;amp;amp;nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 你好，这个问题没有解决。 关于影像坐标的问题，我先给正射校正后的影像加上了投影，然后再拿带有坐标的影像进行的大气校正。 附件是我的代码，供您参考。 &amp;amp;amp;amp;amp;amp;nbsp;

------------------&amp;amp;amp;amp;amp;amp;nbsp;原始邮件&amp;amp;amp;amp;amp;amp;nbsp;------------------ 发件人:"Zhaoguanhua/AtmosphericCorrection"***@***.***&amp;amp;amp;amp;amp;amp;gt;; 发送时间:&amp;amp;amp;amp;amp;amp;nbsp;2021年12月22日(星期三) 下午4:05 ***@***.***&amp;amp;amp;amp;amp;amp;gt;; ***@***.******@***.***&amp;amp;amp;amp;amp;amp;gt;; 主题:&amp;amp;amp;amp;amp;amp;nbsp;Re: [Zhaoguanhua/AtmosphericCorrection] 校正GF-2数据后，水体的校正结果与FLAASH有较大差异 (Issue #9) 你好，请问你的问题解决了吗，能否将你的代码发送给我，我下载这个作者的代码运行的影像没有坐标 — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***&amp;amp;amp;amp;amp;amp;gt; — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.Message ID: ***@***.***&amp;amp;amp;amp;amp;gt; — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***&amp;amp;amp;amp;gt; — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.Message ID: ***@***.***&amp;amp;amp;gt; — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***&amp;amp;gt; — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.Message ID: ***@***.***&amp;gt; — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: ***@***.***&gt;