遥感+python 1.2 辐射定标

目录

• 遥感+python 1.2 辐射定标
• • 一、辐射定标概念
  • • 绝对定标
    • 相对定标
  • 二、辐射定标原理
  • 三、代码实现
  • • 1、数据准备
    • 2、参数获取
    • 3、空值判断
    • 4、辐射定标
    • 5、源码下载

  本章节，笔者主要讲述辐射定标的概念，原理，即代码实现。

  如同尺子上有度量标准、秤有计量标准一样，卫星观测需要有个精确的数据标准。为了让卫星观测到的数据更加真实地反映实际物理量，需要对卫星观测的数据进行定标。

一、辐射定标概念

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

绝对定标

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

相对定标

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

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

二、辐射定标原理

  将记录的原始DN值转换为表观反射率²；将图像的亮度灰度值转换为绝对辐射亮度。
L λ = G a i n ∗ D N + O f f s e t . L_{\lambda} = Gain*DN+Offset. Lλ​=Gain∗DN+Offset.
式中，L为辐射亮度，Gain为辐射增益，Offset为辐射偏置。

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

三、代码实现

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 = Noneoffset = Noneif 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 = Noneoffset = Nonereturn gain, offset

2、参数获取

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

def getKwargs(key, kwargs, default=None, debug=True):try:value = kwargs[key]except KeyError as e:value = defaultif 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.RasterXSizewidth = dataset.RasterYSizebandCount = dataset.RasterCountif xsize is None:xsize = heightif ysize is None:ysize = widthxoff = max(0, xoff)yoff = max(0, yoff)if xoff + xsize > height:xsize = height - xoffif yoff + ysize > width:ysize = width - yoffdata = 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_inodata = bandCount * nodatareturn 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.RasterXSizewidth = inDataset.RasterYSizebandCount = inDataset.RasterCountprint(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)

5、源码下载

以上是部分代码展示，详情请下载。

<返回目录

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！