遥感+python 1.2 辐射定标

QYD2001

已于 2023-08-02 14:43:44 修改

阅读量1.7k

点赞数 3

分类专栏：遥感+python 文章标签： python gdal

于 2023-06-16 17:02:16 首次发布

本文链接：https://blog.csdn.net/Harry233/article/details/131248391

版权

遥感+python 专栏收录该内容

7 篇文章 10 订阅

订阅专栏

本文介绍了遥感图像辐射定标的概念，包括绝对定标和相对定标，并重点讲解了如何使用Python进行辐射定标的代码实现，涉及数据准备、参数获取、空值判断和辐射校正的步骤。示例代码展示了GF-1卫星数据的辐射定标过程。

摘要由CSDN通过智能技术生成

遥感+python 1.2 辐射定标

一、辐射定标概念

辐射定标是用户需要计算地物的光谱反射率或光谱辐射亮度时，或者需要对不同时间、不同传感器获取的图像进行比较时，都必须将图像的亮度灰度值转换为绝对的辐射亮度，这个过程就是辐射定标¹。

绝对定标

通过各种标准辐射源，在不同波谱段建立成像光谱仪入瞳处的光谱辐射亮度值与成像光谱仪输出的数字量化值之间的定量关系。

相对定标

确定场景中各像元之间、各探测器之间、各波谱之间以及不同时间测得的辐射量的相对值。

在这两类辐射定标之中，常见的为绝对定标，本文也将着重实现绝对定标。

二、辐射定标原理

将记录的原始DN值转换为表观反射率²；将图像的亮度灰度值转换为绝对辐射亮度。
$L_{\lambda} = Gain*DN+Offset.$
式中，L为辐射亮度，Gain为辐射增益，Offset为辐射偏置。

其中辐射增益和辐射偏置均为官方提供数据，可以从卫星官网等检索。
例如下所示：

卫星载荷	波段号	辐射增益	辐射偏置
GF-1 PMS	Band1	0.2082	4.6186
GF-1 PMS	Band2	0.1672	4.8768
GF-1 PMS	Band3	0.1748	4.8924
GF-1 PMS	Band4	0.1883	-9.4771

三、代码实现

1、数据准备

在此利用.json文件作为数据对辐射增益等信息进行记录。
获取方式如下：
{}Parameter >{}GF1 >{} WFV1 > {}2013 >[ ]gain
具体如下所示：

{
    "Parameter": {
      "GF1": {
        "WFV1": {
          "2013": {
            "gain": [
              5.851,
              7.153,
              8.368,
              7.474
            ],
            "offset": [
              0.0039,
              0.0047,
              0.003,
              0.0274
            ]
          },
          "SR": {
            "1": [
              0.450,
              0.520
            ],
            "2": [
              0.520,
              0.590
            ],
            "3": [
              0.630,
              0.690
            ],
            "4": [
              0.770,
              0.890
            ]
          }
        }
      }
    }
  }

获取此文件，通常用Python的json库进行数据获取。利用卫星的载荷平台、传感器类型、年份、影像类型。

def RadiometricCalibration(config_file, **kwargs):
    config = json.load(open(config_file, encoding='utf-8-sig'))
    SatelliteID = getKwargs("SatelliteID", kwargs)
    SensorID = getKwargs("SensorID", kwargs)
    Year = getKwargs("Year", kwargs)
    ImageType = getKwargs("ImageType", kwargs, debug=False)
    gain = None
    offset = None
    if SatelliteID in ["GF1"]:
        if SensorID in ["WFV1", "WFV2", "WFV3", "WFV4"]:
            try:
                gain = config["Parameter"][SatelliteID][SensorID][Year]["gain"]
                offset = config["Parameter"][SatelliteID][SensorID][Year]["offset"]
            except KeyError:
                print(f"Cannot find key {SatelliteID}-{SensorID}-{Year}-gain in :{config_file}")
                gain = None
                offset = None
    return gain, offset

2、参数获取

其中考虑到不同卫星所获取的参数不同，采用多参数的方式，对函数进行参数输入，其中要对输入参数进行初步的数据处理，处理方式如下，以防数据无法获取，利用异常处理，设定默认值，对无法获取参数的异常情况做处理，并为调试期设定debug，保证调试信息输出。

def getKwargs(key, kwargs, default=None, debug=True):
    try:
        value = kwargs[key]
    except KeyError as e:
        value = default
        if default is None:
            if debug:
                print(f"\033[31m cant find value by {kwargs}:{e}\033[0m")
    return value

对于影像的信息获取，在此采用遥感影像产品名称解析的方法，对影像信息进行初步提取。

def getInfoFromFilename(filename: str, type=None):
    infoDict = {}
    infoList = filename.split("_")
    if type is None:
        type = infoList[0]
    if type in ["GF1"]:
        infoDict["SatelliteID"] = infoList[0]
        infoDict["SensorID"] = infoList[1]
        infoDict["lon"] = float(infoList[2][1:])
        infoDict["lat"] = float(infoList[3][1:])
        infoDict["Year"] = infoList[4][:4]
        infoDict["month"] = infoList[4][4:6]
        infoDict["day"] = infoList[4][6:8]
    return infoDict

3、空值判断

遥感影像在处理时，存在部分空值、无效值，需要在代码中进行排除。在此采用计算满幅率的方法，对影像进行真实值判断。
注：后续参数为大影像切片计算使用，在此不做处理。

def getFullData(dataset: gdal.Dataset, xoff: int = 0, yoff: int = 0, xsize: int = None, ysize: int = None,
                nodata: float = 0):
    height = dataset.RasterXSize
    width = dataset.RasterYSize
    bandCount = dataset.RasterCount
    if xsize is None:
        xsize = height
    if ysize is None:
        ysize = width
    xoff = max(0, xoff)
    yoff = max(0, yoff)
    if xoff + xsize > height:
        xsize = height - xoff
    if yoff + ysize > width:
        ysize = width - yoff
    data = dataset.GetRasterBand(1).ReadAsArray(xoff, yoff, xsize, ysize)
    for i_bandNo in range(2, bandCount):
        data_i = dataset.GetRasterBand(i_bandNo).ReadAsArray(xoff, yoff, xsize, ysize)
        data = data + data_i
    nodata = bandCount * nodata
    return data != nodata

4、辐射定标

最后进行辐射定标。

def RadiationCalibration(tif_path, out_path, RCP_file, nodata=0):
    """
    辐射定标模块：对影像进行绝对辐射定标。

    :param tif_path: 输入影像数据路径（tif格式文件路径）。
    :param out_path: 输出影像数据路径（tif格式文件路径）。
    :param RCP_file: 辐射定标参数文件路径。
    :param nodata:  空值设置，默认0。
    """
    tif_path = pathlib.Path(tif_path).resolve()
    out_path = pathlib.Path(out_path).resolve()
    RCP_file = pathlib.Path(RCP_file).resolve()
    infoDict = getInfoFromFilename(tif_path.stem)
    if not out_path.parent.exists():
        out_path.parent.mkdir(parents=True, exist_ok=True)
    inDataset = gdal.Open(str(tif_path))
    im_geotrans = inDataset.GetGeoTransform()
    im_proj = inDataset.GetProjection()
    height = inDataset.RasterXSize
    width = inDataset.RasterYSize
    bandCount = inDataset.RasterCount
    print(f"{bandCount} x {height} x {width}")
    driver = gdal.GetDriverByName("GTiff")
    outDataset = driver.Create(str(out_path), height, width, bandCount, GDT_UInt16)
    outDataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
    outDataset.SetProjection(im_proj)  # 写入投影
    Image_block_full = getFullData(inDataset, nodata=nodata)
    gain, offset = RadiometricCalibration(str(RCP_file), **infoDict)
    for i_bandNo in trange(0, bandCount):
        Image_block_band = inDataset.GetRasterBand(i_bandNo + 1).ReadAsArray()
        Image_block_band = Image_block_band * gain[i_bandNo] + offset[i_bandNo]
        Image_block_band = np.where(Image_block_full, Image_block_band, nodata)
        Image_block_band = Image_block_band.astype(np.uint16)
        outband = outDataset.GetRasterBand(i_bandNo + 1)
        outband.WriteArray(Image_block_band)