Featured image of post 6S大气校正的Python实现
遥感&算法

6S大气校正的Python实现

基于GDAL，Py6s包实现遥感数据辐射定标和大气校正.

前言

毕业设计做了一个Landsat数据的辐射定标、大气校正、裁剪、NDVI等指数计算的流程化批处理，现在任务上要实现GF-7数据的校正，其中有些内容需要修改，正好趁这个机会把这整个内容梳理一遍。6S大气校正python包名字是Py6s，该包由南安普顿的博士R.T. Wilson开发并维护，具体的介绍和使用可以自行查看，官网内有规范的示例和实际应用案例。目前包内集成了Landsat数据的各种传感器的处理参数，直接可以调用，但是国产卫星的校正参数如增益 (gain) 和偏置值(bias)，只在中国资源卫星应用中心有发布，该网站还发布了校正用的光谱响应函数(SRF)，需要自己整理起来才能进行6s校正处理。至于更多大气叫校正的原理什么的，有需要的话自己查一下论文或者其他博客吧，这里就推荐几个，比如遥感之心、冰柚子 等。

正文

站在巨人的肩膀上，才能更好地前进。国内前人已经有人做过python的6s大气校正工程化模块，作者是武汉大学遥感国家重点实验室的博士，在CSDN做博客记录，在github仓库内查看。该仓库的文件可以通过运行相对应的python文件直接进行校正处理，但是目前脚本也只能够做到单幅处理。但是在该脚本基础上，江北通过将待处理的多幅影像存储到txt文件中，利用txt文件进行BAT脚本批量处理。虽然也变相实现了6s批量大气校正，但个人感觉不太pythonic，所以还是想用自己的想法过一遍批处理的纯python代码。

次要函数

先把几个次要函数，为主函数服务的函数先给上。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42


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)

这两个函数基本没有改动，只改了untar为传参制。

主函数

以下介绍几个主要函数的改动。

辐射定标

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16


def RadiometricCalibration(BandId,SatelliteID,SensorID,Year,config):
    ''' 获取辐射校正增益和偏置值；
    BANDID : 第n波段
    SatelliteID : 卫星平台，存放在json内
    SensorID : 传感器类型，存放在json内
    Year : 定标系数发布时间，存放在json内
    config : 定标参数存储json文件路径 '''
    if SensorID[0:3] == "bwd":
        Gain_ =config["Parameter"][SatelliteID][SensorID][Year]["gain"]
        Bias_ =config["Parameter"][SatelliteID][SensorID][Year]["offset"]
    else:
        Gain_ =config["Parameter"][SatelliteID][SensorID][Year]["gain"][BandId-1]
        Bias_ =config["Parameter"][SatelliteID][SensorID][Year]["offset"][BandId-1]

    return Gain_,Bias_

大体是和原脚本一样的，但是为了用于高分数据不同传感器，这里改了判别函数。未来会加入前视传感器(FWD)的判别和json参数，其实也就是加elif形成多判别就行了。

大气校正

接下来是大气校正：

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114


def AtmosphericCorrection(BandId,metedata,config,SatelliteID,SensorID,GMTEDdem):
    ''' 
    大气校正函数
    '''
    dom = xml.dom.minidom.parse(metedata)

    # 6S模型
    s = SixS()

    # 传感器类型 自定义
    s.geometry = Geometry.User()
    s.geometry.solar_z = float(dom.getElementsByTagName('SunAltitude')[0].firstChild.data)
    s.geometry.solar_a = float(dom.getElementsByTagName('SunAzimuth')[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)
    # 中心经纬度\
        
    LeftTop = dom.getElementsByTagName('LeftTopPoint')[0]
    TopLeftLat= float(LeftTop.getElementsByTagName('Latitude')[0].firstChild.data)    
    TopLeftLon= float(LeftTop.getElementsByTagName('Longtitude')[0].firstChild.data)   
        
    RightTop = dom.getElementsByTagName('RightTopPoint')[0]
    TopRightLat= float(RightTop.getElementsByTagName('Latitude')[0].firstChild.data)    
    TopRightLon= float(RightTop.getElementsByTagName('Longtitude')[0].firstChild.data)   
    
    BottomRight = dom.getElementsByTagName('RightBottomPoint')[0]
    BottomRightLat= float(BottomRight.getElementsByTagName('Latitude')[0].firstChild.data)    
    BottomRightLon= float(BottomRight.getElementsByTagName('Longtitude')[0].firstChild.data)   
    
    Bottomleft = dom.getElementsByTagName('LeftBottomPoint')[0]
    BottomLeftLat= float(Bottomleft.getElementsByTagName('Latitude')[0].firstChild.data)    
    BottomLeftLon= float(Bottomleft.getElementsByTagName('Longtitude')[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, GMTEDdem)) * 0.001

    # 研究区海拔、卫星传感器轨道高度
    s.altitudes = Altitudes()
    s.altitudes.set_target_custom_altitude(meanDEM)
    s.altitudes.set_sensor_satellite_level()

    # 校正波段（根据波段名称）全色440-905
    if BandId == 1:
        SRFband = config["Parameter"]["GF1"]["WFV1"]["SRF"]["1"]
        s.wavelength = Wavelength(0.450,0.520,SRFband)

    elif BandId == 2:
        SRFband = config["Parameter"]["GF1"]["WFV1"]["SRF"]["2"]

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

    elif BandId == 3:
        SRFband = config["Parameter"]["GF1"]["WFV1"]["SRF"]["3"]

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

    elif BandId == 4:
        SRFband = config["Parameter"]["GF1"]["WFV1"]["SRF"]["4"]
        s.wavelength = Wavelength(0.770,0.890,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)

在这里值得一提的是，GF1\GF2中xml的太阳天顶角（SolarZenith）和太阳方位角（SolarAzimuth）参数在高分较新的数据中，xml的标签被改为了SunAltitude和SunAzimuth也就是天顶角被改为了高度角，这样原脚本中

1

    s.geometry.solar_z = 90-float(dom.getElementsByTagName('SolarZenith')[0].firstChild.data)

的写法，就要把90-去掉了具体原理可以查看这里。

另外一点是，在GF6以上的数据中，头文件xml的影像顶点坐标经度、纬度标签被归到了各自的位置标签上，对比原写法和新写法：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

    #原写法
    # 中心经纬度
    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)

    #新写法
    # 中心经纬度  
    LeftTop = dom.getElementsByTagName('LeftTopPoint')[0]
    TopLeftLat= float(LeftTop.getElementsByTagName('Latitude')[0].firstChild.data)    
    TopLeftLon= float(LeftTop.getElementsByTagName('Longtitude')[0].firstChild.data)   
        
    RightTop = dom.getElementsByTagName('RightTopPoint')[0]
    TopRightLat= float(RightTop.getElementsByTagName('Latitude')[0].firstChild.data)    
    TopRightLon= float(RightTop.getElementsByTagName('Longtitude')[0].firstChild.data)   
    
    BottomRight = dom.getElementsByTagName('RightBottomPoint')[0]
    BottomRightLat= float(BottomRight.getElementsByTagName('Latitude')[0].firstChild.data)    
    BottomRightLon= float(BottomRight.getElementsByTagName('Longtitude')[0].firstChild.data)   
    
    Bottomleft = dom.getElementsByTagName('LeftBottomPoint')[0]
    BottomLeftLat= float(Bottomleft.getElementsByTagName('Latitude')[0].firstChild.data)    
    BottomLeftLon= float(Bottomleft.getElementsByTagName('Longtitude')[0].firstChild.data)

其实相当于拐了个弯，但是也要写出来提醒其他开发者，要注意查看影像的xml，才能及时提取对的内容。

分块

有了以上两个最重要的函数后，可以进行分块处理了。这里也是直接用的原脚本中的写法，但是多了几个判别函数，实现自动判别遥感平台和传感器，从而正确调用参数：

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108

def Block(IDataSet,SensorID,Year,SatelliteID,atcfiles,config,metedata,outFileName, GMTEDdem):
    ''' 分块处理函数
    以分块运行效率，2分钟之内可校正一幅GF-7-MUX数据
    '''
    # global cols,rows,atcfiles
    cols = IDataSet.RasterXSize
    rows = IDataSet.RasterYSize

   

    #设置输出波段
    Driver = IDataSet.GetDriver()
    geoTransform1 = IDataSet.GetGeoTransform()
    ListgeoTransform1 = list(geoTransform1)
    ListgeoTransform1[5] = -ListgeoTransform1[5]
    newgeoTransform1 = tuple(ListgeoTransform1)
    proj1 = IDataSet.GetProjection()

    OutRCname = os.path.join(atcfiles,outFileName+".tif")

    outDataset = Driver.Create(OutRCname,cols,rows,4,gdal.GDT_Int32)
    outDataset.SetGeoTransform(newgeoTransform1)
    outDataset.SetProjection(proj1)
    
    #分别读取4个波段，多波段文件校正
    if SensorID=="MUX":
        for m in range(1,5):

            ReadBand = IDataSet.GetRasterBand(m)
            outband = outDataset.GetRasterBand(m)
            outband.SetNoDataValue(-9999)
            #获取对应波段的增益gain和偏移bias
            Gain,Bias = RadiometricCalibration(m,SatelliteID,SensorID,Year,config)

        
            #获取大气校正系数
    
            AtcCofa, AtcCofb, AtcCofc = AtmosphericCorrection( m, metedata, config, SatelliteID, SensorID, GMTEDdem)
            nBlockSize = 2048
            i = 0
            j = 0
            b = cols*rows
        
            print("第%d波段校正："%m)
            try:
                
                while i<rows:
                    while j <cols:
                        #保存分块大小
                        nXBK = nBlockSize
                        nYBK = nBlockSize

                        #最后不够分块的区域，有多少读取多少
                        if i+nBlockSize>rows:
                            nYBK = rows - i
                        if j+nBlockSize>cols:
                            nXBK=cols - j

                        #分块读取影像
                        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)

                        outband.WriteArray(atcImage,j,i)

                        j=j+nXBK
                        time.sleep(1)
                    j=0
                    i=i+nYBK
            except KeyboardInterrupt:
                raise
            
    #单波段文件校正
    elif SensorID == "nad" or SensorID == "bwd" or SensorID == "fwd":
        ReadBand = IDataSet.GetRasterBand(1)
        outband = outDataset.GetRasterBand(1)
        outband.SetNoDataValue(-9999)
        Gain,Bias = RadiometricCalibration(ReadBand,SatelliteID,SensorID,Year,config)
        nBlockSize = 2048
        i = 0
        j = 0
        b = cols*rows
    
        print("单波段校正")

        while i<rows:
            while j <cols:
                #保存分块大小
                nXBK = nBlockSize
                nYBK = nBlockSize

                #最后不够分块的区域，有多少读取多少
                if i+nBlockSize>rows:
                    nYBK = rows - i
                if j+nBlockSize>cols:
                    nXBK=cols - j
                #分块读取影像
                Image = ReadBand.ReadAsArray(j,i,nXBK,nYBK)
                outImage =np.where(Image>0,Image*Gain + Bias,-9999)
                outband.WriteArray(outImage,j,i)
                j=j+nXBK
                time.sleep(1)
            j=0
            i=i+nYBK
        print("已完成辐射定标,不需要大气校正")

这里写法仍然不够严谨，虽然全色波段大气校正效果并没有那么显著，但是由于官方也给了光谱响应函数、增益偏置，因此理应也要进行判别，从而进行辐射和大气校正，但是这里偷了个懒，没有写完所有传感器判别。

主函数

最后为了更加简洁，使参数能够更灵活地输入，这里自己写了一个主函数:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

def read_output_file(input_tarfile,config_file,OutputFilePath,outFileName,GMTEDdem):
    ''' 
    主函数
    input_tarfile: tar压缩文件路径；
    config_file: 辐射校正参数文件路径；
    OutputFilePath: 输出文件保存路径；
    outFileName: 输出文件保存名字；
    GMTEDdem: 校正所需的DEM文件。
    '''
    config = json.load(open(config_file))
    filename = os.path.basename(input_tarfile)
    fileType = filename[0:3]
    filename_split = filename.split("_")
    base_dir = os.path.dirname(input_tarfile)
    
    if fileType == 'GF7':
        #GFType = filename[4:7]
        #GFType = filename_split[5]
        #GF701_011628_E107.2_N22.9_20211202114945_BWD_01_SC0_0001_2112039662
        
        SensorID = filename_split[5]
        Year = filename_split[4][:4]
        SatelliteID = filename_split[0]
        untarfilename = filename[:-7]
        outpath = os.path.join(base_dir,untarfilename)
        atcfiles = os.path.join(OutputFilePath,outFileName)
        untar(input_tarfile, outpath)
        
        tiffFile = glob.glob(os.path.join(outpath,"*.tif"))[0]
        metedata = glob.glob(os.path.join(outpath,"*.xml"))[0]
        rpb_file = glob.glob(os.path.join(outpath,"*_rpc.txt"))[0]

        try:
            os.mkdir(atcfiles)
        except Exception as e:
            pass
    
        #将头文件和几何校正文件拷贝到大气校正结果文件中
        metedata_basename = os.path.basename(metedata)
        copy_metedata = os.path.join(atcfiles,metedata_basename)
        shutil.copy(metedata,copy_metedata)

        rpb_basename=os.path.basename(rpb_file)
        copy_rpb_file = os.path.join(atcfiles,rpb_basename)
        shutil.copy(rpb_file,copy_rpb_file)

        try:
            IDataSet = gdal.Open(tiffFile)
        except Exception as e:
            print("文件%S打开失败" % tiffFile)
   
    elif fileType == 'zy3':  
            #GFType = filename[4:7]
        #GFType = filename_split[5]
        SensorID = filename_split[1]
        Year = filename_split[4][:4]
        SatelliteID = filename_split[0]
        
        untarfilename = filename[:-7]
        outpath = os.path.join(base_dir,untarfilename)
        atcfiles = os.path.join(OutputFilePath,outFileName)
        untar(input_tarfile, outpath)
            
        tiffFile = glob.glob(os.path.join(outpath,"*.tif"))[0]
        metedata = glob.glob(os.path.join(outpath,"*.xml"))[0] 
        rpb_file = glob.glob(os.path.join(outpath,"*_rpc.txt"))[0]

        try:
            os.mkdir(atcfiles)
        except Exception as e:
            pass
        
        #将头文件和几何校正文件拷贝到大气校正结果文件中
        metedata_basename = os.path.basename(metedata)
        copy_metedata = os.path.join(atcfiles,metedata_basename)
        shutil.copy(metedata,copy_metedata)

        rpb_basename=os.path.basename(rpb_file)
        copy_rpb_file = os.path.join(atcfiles,rpb_basename)
        shutil.copy(rpb_file,copy_rpb_file)

        try:
            IDataSet = gdal.Open(tiffFile)
        except Exception as e:
            print("文件%S打开失败" % tiffFile)
    Block(IDataSet,SensorID,Year,SatelliteID,atcfiles,config,metedata,outFileName,GMTEDdem)

写法思路也是参照脚本中的if __main__部分，将其封装成函数以后，直接给出函数的 5 参数，就能方便地在python中进行定制化的批处理或者其他流程化应用。

1
2
3
4
5
6

Dem = R".GMTED2km.tif"
outFileName = ".123"
OutputFilePath=r".outtar"
config_file = r".RadiometricCorrectionParameter.json"
input_tarfile = r".GFXXXXX.tar.gz"
read_output_file(input_tarfile, config_file, OutputFilePath, outFileName, Dem)

辅助数据

大气校正是复杂的物理机制，其中需要的参数是繁多的，6s包给定了大部分参数，仍有一些参数需要自己输入，比如DEM、光谱响应函数，这里介绍辅助数据的使用。

DEM

高程数据就很简单，如果需要高精度的可以直接上高分辨率dem，原脚本中使用的是美国USGS的GMTED高程数据，分辨率为2km，如果有条件，直接用SRTM或者日本ALOS2也可以，但是要保证覆盖完研究区。

辐射定标和光谱响应json文件

辐射定标参数和光谱响应函数所需要的参数都以json格式存放在RadiometricCorrectionParameter.json文件中，其结构如下所示：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

{
    "Parameter": {
        "GF1": {
            "WFV1": {
                "2013": {
                    "gain": [
                        5.851,
                        7.153,
                        8.368,
                        7.474
                    ],
                    "offset": [
                        0.0039,
                        0.0047,
                        0.003,
                        0.0274
                    ]
                },
                "SRF": {
                    "1": [
                        0.46233000000000113,
                        ...
                        0.10741000000000049
                    ],
                    "2": [
                        0.7663399999999977,
                    ....
                    ]
                }
            },

这里要注意，Py6s包的光谱响应函数:

1

s.wavelength = Wavelength(0.770,0.890,SRFband)

其中校正波段范围为0.770~0.890nm，其数据分辨率是2.5nm，而官方给的参数往往是1nm分辨率，因此需要插值到0.5，然后间隔2.5取值，这也是后续需要进行的内容。在卫星中心官网下载到需要的SRF和辐射定标参数后，按json数据格式整理到其中，就能方便调用了，根据不同的传感器和遥感平台，将数据添加并定制化即可。

后记

目前开源的大气校正算法，仅有6s算是比较靠谱的了，而python的封装让大气校正的实现变得简单了起来，因此也对南安普顿威廉森教授致谢，也感谢赵博士的开源。后续考虑加入GF和ZY系列传感器的判别，整理一下json然后commit到github。另一方面，国产卫星的快速大气校正和深度学习大气校正也需要进行研究，目前想法就这些，如果有感兴趣的内容想与我讨论，欢迎发送邮件和我交流。

© 2022 - 2023 jerrychoices 的博客
Author by Jerrychoices
Built with Hugo
主题 Stack 由 Jimmy 设计

本站访客数人次 总访问量次 本文阅读量次