【代码】基于6s模型的影像大气校正工程

这个工程，不是我写的，是赵前辈在github的开源代码， 调用py6s接口，自动读取影像头文件信息，对遥感影像进行大气校正批处理。

地址是：https://github.com/Zhaoguanhua/AtmosphericCorrection

各位有空可以去细读这个工程。

针对GF1、2影像， landsat8影像， Sentinel影像， 已经可以工程化使用。

市面上的大气校正方法大都是基于 6S 。ENVI是独一家开发了 FLAASH大气校正。

在不完全的个人统计下， FLAASH 比6s的效果会好一些。

下面我讲解一下，怎么去理解这份工程。

6s原理的快速理解

这里不细讲6s原理，如果需要研究6S大气校正的细节，请下载原论文。

此外，对于国产卫星， 需要去获取校正参数如 增益 ( gain ) 和 偏置值( bias ) ，只在 中国资源卫星应用中心 有发布，该网站还发布了校正用的 光谱响应函数 ( SRF )，需要自己整理起来才能进行6s校正处理。

代码逻辑

核心逻辑是：

输入参数->6S->得到a\b\c三个参数 ->线性公式 ->地表反射率

文字逻辑是：

通过6s获得三个参数分别是a,b,c

公 式可以表示为：

辐 射 率 = DN * gain + bias

y = 辐 射 率 *a-b

地表反射率 = y/(1+y*c）

请看下面的具体代码：

Image = ReadBand.ReadAsArray(j,i,nXBK,nYBK)outImage =np.where(Image>0,Image*Gain + Bias,-9999)y = np.where(outImage!=-9999,AtcCofa * outImage - AtcCofb,-9999)atcImage = np.where(y!=-9999,(y / (1 + y * AtcCofc))*10000,-9999)

这段代码是在做遥感图像的大气校正，具体来说，它使用了一个线性大气校正公式和一个大气校正系数来修正图像的反射率。 这段代码的大致流程如下：

1. 从遥感图像中读取一个块（band）。这个块的大小由nXBK和nYBK决定，而i和j决定了这个块的位置。

2. 对读取的图像块进行初步的校正。其中，Gain和Bias是用来进行线性校正的参数，它们的值通常由遥感图像的元数据给出。这个校正的目的是将图像的数字数值（DN）转换为反射率。

3. 对校正后的图像进行大气校正。其中，AtcCofa、AtcCofb和AtcCofc是大气校正的系数，它们的值通常由大气校正模型（如6S模型）计算得到。这个大气校正的公式是一个线性公式，它的目的是去除大气对图像反射率的影响。

4. 最后，将大气校正后的图像的反射率乘以10000，以便于存储和处理。这是因为反射率是一个小于1的小数，乘以10000后可以转换为整数，便于存储和处理。

在这个过程中，对于图像中的无效值（如云层或无数据区域），用-9999来表示。这是一个常用的方式，可以方便地识别和处理这些无效值。

所以你知道为什么经过大气校正后，数值会扩大10000倍了么？其实envi的Flaash也会 扩大 1 0000 倍来保存结果。

我们现在来讲解一下代码。

代码讲解

辅助函数代码，自动解压文件，自动根据经纬度获得DEM的高度

def MeanDEM(pointUL, pointDR, GMTEDdem):    '''    计算影像所在区域的平均高程.    '''    try:        DEMIDataSet = gdal.Open(GMTEDdem)    except Exception as e:        pass
    DEMBand = DEMIDataSet.GetRasterBand(1)    geotransform = DEMIDataSet.GetGeoTransform()    # DEM分辨率    pixelWidth = geotransform[1]    pixelHight = geotransform[5]
    # DEM起始点：左上角，X：经度，Y：纬度    originX = geotransform[0]    originY = geotransform[3]
    # 研究区左上角在DEM矩阵中的位置    yoffset1 = int((originY - pointUL['lat']) / pixelWidth)    xoffset1 = int((pointUL['lon'] - originX) / (-pixelHight))
    # 研究区右下角在DEM矩阵中的位置    yoffset2 = int((originY - pointDR['lat']) / pixelWidth)    xoffset2 = int((pointDR['lon'] - originX) / (-pixelHight))
    # 研究区矩阵行列数




    
    xx = xoffset2 - xoffset1    yy = yoffset2 - yoffset1    # 读取研究区内的数据，并计算高程    DEMRasterData = DEMBand.ReadAsArray(xoffset1, yoffset1, xx, yy)
    MeanAltitude = np.mean(DEMRasterData)    return MeanAltitude
def untar(fname, dirs):    ''' 解压缩tar文件函数 '''
    t = tarfile.open(fname)    t.extractall(path=dirs)

这里定义了两个函数， MeanDEM 和 untar 。

1. MeanDEM 函数用于计算给定区域的平均高程。函数首先打开一个DEM（数字高程模型）数据集，然后根据输入的经纬度信息（ pointUL 和 pointDR ），计算出这个区域在DEM数据集中的位置和大小。然后，函数读取这个区域的DEM数据，并计算其平均值。这个平均值就是这个区域的平均高程。函数的输入参数包括：

• pointUL ：研究区域的左上角的经纬度。

• pointDR ：研究区域的右下角的经纬度。

• GMTEDdem ：DEM数据集的文件路径。

• fname ：需要解压缩的tar文件的文件路径。

• dirs ：解压缩的目标目录。

这两个函数都使用了 try - except 语句来处理可能出现的异常，以确保程序在出现错误时不会立即崩溃，而是可以继续执行下去。

这个是获取6S的参数。

def AtmosphericCorrection(BandId,metedata,config,SatelliteID,SensorID):    # global metedata,config,SatelliteID,SensorID    #读取头文件    dom = xml.dom.minidom.parse(metedata)
    # 6S模型    s = SixS()
    # 传感器类型 自定义    s.geometry = Geometry.User()    s.geometry.solar_z = 90-float(dom.getElementsByTagName('SolarZenith')[0].firstChild.data)    s.geometry.solar_a = float(dom.getElementsByTagName('SolarAzimuth')[0].firstChild.data)    # s.geometry.view_z = float(dom.getElementsByTagName('SatelliteZenith')[0].firstChild.data)    # s.geometry.view_a = float(dom.getElementsByTagName('SatelliteAzimuth')[0].firstChild.data)    s.geometry.view_z = 0    s.geometry.view_a = 0    # 日期    DateTimeparm = dom.getElementsByTagName('CenterTime')[0].firstChild.data    DateTime = DateTimeparm.split(' ')    Date = DateTime[0].split('-')    s.geometry.month = int(Date[1])    s.geometry.day = int(Date[2])
    # print(s.geometry)    # 中心经纬度    TopLeftLat = float(dom.getElementsByTagName('TopLeftLatitude')[0].firstChild.data)    TopLeftLon = float(dom.getElementsByTagName('TopLeftLongitude')[0].firstChild.data)    TopRightLat = float(dom.getElementsByTagName('TopRightLatitude')[0].firstChild.data)    TopRightLon = float(dom.getElementsByTagName('TopRightLongitude')[0].firstChild.data)    BottomRightLat = float(dom.getElementsByTagName('BottomRightLatitude')[0].firstChild.data)    BottomRightLon = float(dom.getElementsByTagName('BottomRightLongitude')[0].firstChild.data)    BottomLeftLat = float(dom.getElementsByTagName('BottomLeftLatitude')[0].firstChild.data)    BottomLeftLon = float(dom.getElementsByTagName('BottomLeftLongitude')[0].firstChild.data)
    ImageCenterLat = (TopLeftLat + TopRightLat + BottomRightLat + BottomLeftLat) / 4
    # 大气模式类型    if ImageCenterLat > -15 and ImageCenterLat < 15:        s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical)
    if ImageCenterLat > 15 and ImageCenterLat < 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 ImageCenterLat > 45 and ImageCenterLat < 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气溶胶光学厚度,对应能见度为40km    s.aot550 = 0.14497
    # 通过研究去区的范围去求DEM高度。    pointUL = dict()    pointDR = dict()    pointUL["lat"] = max(TopLeftLat,TopRightLat,BottomRightLat,BottomLeftLat)    pointUL["lon"] = min(TopLeftLon,TopRightLon,BottomRightLon,BottomLeftLon)    pointDR["lat"] = min(TopLeftLat,TopRightLat,BottomRightLat,BottomLeftLat)    pointDR["lon"] = max(TopLeftLon,TopRightLon,BottomRightLon,BottomLeftLon)    meanDEM = (MeanDEM(pointUL, pointDR)) * 0.001
    # 研究区海拔、卫星传感器轨道高度