高分影像批处理第二回——辐射定标与大气校正

最新推荐文章于 2026-05-15 05:11:08 发布

原创最新推荐文章于 2026-05-15 05:11:08 发布 · 4.9k 阅读

51 ·

本内容遵循CC 4.0 BY-SA版权协议

标签

#python #vscode #windows #其他

GF影像处理专栏收录该内容

4 篇文章

订阅专栏

Python3.8

Python 是一种高级、解释型、通用的编程语言，以其简洁易读的语法而闻名，适用于广泛的应用，包括Web开发、数据分析、人工智能和自动化脚本

辐射定标

遥感影像的辐射定标就是将传感器接收的无效的DN（digital number）值转化为我们所需要的辐射亮度值或者大气表观反射率。要达到这个目的，我们需要知道传感器的增益和偏移值。高分影像的gain和offset都可以在官网查到。

$RF = DN*Gain + Offset$

定标系数http://218.247.138.119/CN/Downloads/dbcs/index.shtml

如图所示，即可得到不同年份不同传感器的定标值。

Dataset = gdal.Open(r'xxx.tif')
outDataset = Driver.Create(OutRCname,cols,rows,m,gdal.GDT_Int32)
outDataset.SetGeoTransform(newgeoTransform1)
outDataset.SetProjection(proj1)
m = Dataset.RasterCount
for i in range(1,m):
    ReadBand = Dataset.GetRasterBand(i)
    outband = outDataset.GetRasterBand(m)
    outband.SetNoDataValue(-9999)
    Image = ReadBand.ReadAsArray()
    outImage =np.where(Image>0,Image*Gain + Bias,-9999)
    outband.WriteArray(outImage)

如上图代码所示，只需要简易的代码处理，即可实现传感器的辐射定标工作。

大气校正

辐射在传播的过程中会受到许多现象的干扰，其中大气的散射、折射、反射影响尤为重要，为了得到更加真实的地表辐射或者辐射率，需要进行大气校正。假设一种大气状况，构建辐射的传输模型，最后定量的模拟大气的影响以还原地表的辐射。

平常使用最多的大气校正工具就是ENVI的flaash和quik大气校正模型，但是批量工具如果想调用就必须要IDL程序设计语言。因此只能另辟蹊径，使用另外一种开源的大气校正模型——6s大气校正模型。关于6s的大气校正模型介绍可查阅相关的论文。

(11条消息) Py6S安装教程_qiaofeng0419的博客-CSDN博客_安装py6shttps://blog.csdn.net/qiaofeng0419/article/details/104136012?spm=1001.2101.3001.6650.6&utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7EBlogCommendFromBaidu%7ERate-6.pc_relevant_paycolumn_v3&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7EBlogCommendFromBaidu%7ERate-6.pc_relevant_paycolumn_v3&utm_relevant_index=6py6s（pythonAPI）官方文档：

Installation — Py6S 1.9.0 documentationhttps://py6s.readthedocs.io/en/latest/installation.html；除此之外，还有比较关键的波普响应函数，话说，传感器毕竟是人制作的，由于材料的不同，传感器对不同波段的辐射的感受程度也是不一样的。这个参数也可以在高分中心查到。

光谱响应函数http://218.247.138.119/CN/Downloads/gpxyhs/index.shtml作为一个缝合怪，更多的信息我就不多说了，细节可以查看上下述文档。

高分六号卫星影像预处理踩坑记录（二）——如何使用python代码预处理卫星影像 - 知乎 (zhihu.com)https://zhuanlan.zhihu.com/p/364889134 跳转中... (github.com)https://github.com/Zhaoguanhua/AtmosphericCorrection

def AtmosphericCorrection(BandId,GFType,CameraType,XmlPath,SATPARAMS):
    # global metedata,config,SatelliteID,SensorID
    #读取头文件
    dom = xml.dom.minidom.parse(XmlPath)
    if BandId =='PAN':
        BandId = 0
    else:
        BandId = eval(BandId[1])
    # 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
    # ImageCenterLat = float(dom.getElementsByTagName('CenterLatitude')[0].firstChild.data)
    # 大气模式类型
    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 自动计算DEM校正
    if global_var.get_value('autoDEMDownload'): 
    # 通过研究区的范围去求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
    
    
    meanDEM = global_var.get_value('meanDEM')
    
    # 研究区海拔、卫星传感器轨道高度
    s.altitudes = Altitudes()
    s.altitudes.set_target_custom_altitude(meanDEM)
    s.altitudes.set_sensor_satellite_level()

    # 校正波段（根据波段名称）
    if BandId == 'B1':
        SRFband = SATPARAMS["Parameter"][GFType][CameraType]["SRF"]["B1"]
        s.wavelength = Wavelength(0.450,0.520,SRFband)

    elif BandId == 'B2':
        SRFband = SATPARAMS["Parameter"][GFType][CameraType]["SRF"]["B2"]

        s.wavelength = Wavelength(0.520,0.590,SRFband)

    elif BandId == 'B3':
        SRFband = SATPARAMS["Parameter"][GFType][CameraType]["SRF"]["B3"]

        s.wavelength = Wavelength(0.630,0.690,SRFband)

    elif BandId == 'B4':
        SRFband = SATPARAMS["Parameter"][GFType][CameraType]["SRF"]["B4"]
        s.wavelength = Wavelength(0.770,0.890,SRFband)
    
    elif BandId == 'PAN':
        SRFband = SATPARAMS["Parameter"][GFType][CameraType]["SRF"]["PAN"]
        
        s.wavelength = Wavelength(0.450,0.900,SRFband)
        
    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)

通过上述的代码过程可以获得大气校正的纠正参数，在利用下面的公式即可进行大气校正。

$y=xa*(measured radiance)-xb; acr=y/(1.+xc*y)$

y = np.where(outImage!=-9999,AtcCofa * outImage - AtcCofb,-9999)
atcImage = np.where(y!=-9999,(y / (1 + y * AtcCofc))*10000,-9999)

这样最终就得到了6s大气校正后的地表反射率。

PS: 由于工作的不断深入，代码在不断的优化，因此完整的代码在章节结束时挂上。

您可能感兴趣的与本文相关的镜像

Python3.8

Conda

Python

Python 是一种高级、解释型、通用的编程语言，以其简洁易读的语法而闻名，适用于广泛的应用，包括Web开发、数据分析、人工智能和自动化脚本