版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
从官网下载的哨兵2 L1C数据是大气表观反射率产品，而大气校正输入数据必须为辐射亮度数据，大气表观反射率和辐射亮度值存在以下数学关系：

转换公式

辐射亮度转换为大气表观反射率公式：

ρλ=πLλd2/ESUNλsinθ\rho_{\lambda}={\pi L_{\lambda}d^2}/{ESUN_{\lambda}sin\theta} ρλ=πLλd2/ESUNλsinθ
大气表观反射率转换为辐射亮度公式
Lλ=ρλESUNλsinθ/πd2\ L_{\lambda}={\rho_{\lambda}ESUN_{\lambda}sin\theta}/{\pi d^2} Lλ=ρλESUNλsinθ/πd2
参数说明：
ρλ：大气表观反射率；Lλ：辐射亮度；d：日地距离；ESUNλ：太阳辐照度；θ:太阳高度角\rho_{\lambda}：大气表观反射率；L_{\lambda}：辐射亮度；d：日地距离；ESUN_{\lambda}：太阳辐照度；\theta:太阳高度角 ρλ：大气表观反射率；Lλ：辐射亮度；d：日地距离；ESUNλ：太阳辐照度；θ:太阳高度角

参数查找

日地距离、太阳辐照度

xml文件

将哨兵2MTD_MSIL1C.xml以notepad++打开，日地距离可从文件中节点为U得知,太阳辐照度—Solar_Irradiance就在下面几行。

IDL读取

栅格元信息中 ‘EARTH SUN DISTANCE’即为日地距离，Solar_Irradiance即为太阳辐照度，这里注意IDL读哨兵2xml数据读出来是3个栅格对象——10米（B2，B3，B4，B8波段）；20米（B5，B6，B7，B8A，B11，B12波段）；60米（B1，B9，B10波段），我这里只对10米进行数据转换，只读取了10米波段的栅格数据，所以这里太阳辐照度顺序对应为—B2,B3,B4,B8.

file = 'G:\L1C_T49SCC_A037624_20220905T033145\S2A_MSIL1C_20220905T032531_N0400_R018_T49SCC_20220905T052134.SAFE\MTD_MSIL1C.xml'
raster = (e.openraster(file))[0]
print,raster.metadata

太阳高度角

xml文件

在相对路径如GRANULE/L1C_T51RTQ_A015101_20180514T024123下包含另一个名为MTD_TL的xml文件，在节点Mean Sun Angle下有两个值，太阳天顶角（ZENITH ANGLE）；太阳方位角（AZIMUTH_ANGLE），太阳高度角和太阳天顶角互余，所以就为90-太阳天顶角。

IDL读取

同样读取10米栅格的元信息，SUN ELEVATION即为太阳高度角。

代码实现

file = 'G:\L1C_T49SCC_A037624_20220905T033145\S2A_MSIL1C_20220905T032531_N0400_R018_T49SCC_20220905T052134.SAFE\MTD_MSIL1C.xml'inopenraster = e.openraster(file)metadata = inopenraster[0].metadatasun_eleva = metadata['SUN ELEVATION']sun_distance = metadata['EARTH SUN DISTANCE'];创建对象存储转换之后的每个波段数据band_arr = objarr(4)for i=1,4 do begin;进行表观反射率转换为辐射亮度转换;公式构建band_SOLAR = (metadata['SOLAR IRRADIANCE'])[i-1]expression1 = 'b'+strtrim(string(i),2) + '*'expression2 = band_SOLAR * sin(sun_eleva/180 * !PI)expression3 = !DPI * sun_distance*sun_distance * 100000expression = expression1+strcompress(string(expression2))+'/'+strcompress(string(expression3))print,expression;进行band math计算band_Task=ENVITask('PixelwiseBandMathRaster')band_Task.INPUT_RASTER = inopenraster[0]band_Task.EXPRESSION = expression
;    ndvi_Task.output_raster_uri = 'G:\Sentinel2_Radio\' + 'b'+strtrim(string(i),2) + '_final5.dat'band_Task.Executeband_arr[i-1] =  band_Task.output_rasterendfor;波段组合gridTask = ENVITask('BuildGridDefinitionFromRaster')gridTask.INPUT_RASTER = inopenraster[0]gridTask.ExecutestackTask = ENVITask('BuildLayerStack')stackTask.INPUT_RASTERS = band_arrstackTask.GRID_DEFINITION = gridTask.OUTPUT_GRIDDEFINITIONstackTask.Executeprint,stackTask.output_raster_uri

对比该结果和ENVI官方哨兵2辐射定标插件结果，出来是一样的，那就没得问题。
参考资料：
https://blog.csdn.net/u013471015/article/details/88663209
https://www.cnblogs.com/enviidl/p/16386359.html

IDL学习——哨兵2 L1C数据辐射定标相关推荐

Sentinel-2数据辐射定标及大气校正
使用工具: 1.SNAP平台及 Sen2cor插件处理过程: 1.下载SNAP平台和Sen2cor插件,建议安装最新版本 2.SNAP平台直接安装: 3.Sen2cor解压之后(文件夹名为:Sen ...
landsat TM数据辐射定标和flaash大气校正
才发现我写了那么多都是废话,想知道landsat TM 辐射定标和大气校正的,请移步这里:http://blog.sina.com.cn/s/blog_764b1e9d0102v59e.html 有视 ...
关于哨兵2号数据辐射定标
最近在做sentinel2的TOA对比,下载到了L1C级别的数据但是DN值不是TOA值查找了相关文档找到了一个计算公式然后就去找相关参数 RADIO_ADD_OFFSET 和 QUANTIFICA ...
遥感IDL二次开发（辐射定标）
1.程序功能: 在IDL环境下进行辐射定标,并将结果返回ENVI. 2.运行步骤: 2.1 在ENVI中打开TM数据,并将数据传回到IDL工作空间: 首先打开ENVI+IDL8.5,导入TM数据,然后 ...
辐射定标与大气校正(ENVI和6s模型对比)
辐射定标与大气校正实验目的与任务实验目的遥感图像通常是用无量纲的数字量化值(DN) 记录信息的,进行遥感定量化分析时,常用到辐射亮度值.反射率值和温度值等物理量.通过辐射定标可以实现DN值与这此 ...
tm影像辐射定标_ENVI+遥感图像的辐射定标
实验:遥感图像的辐射定标 1. 实验目的与任务: ( 1 )了解辐射定标的原理: ( 2 )使用 ENVI 软件自带的定标工具定标 ( 3 )学习使用波段运算进行辐射定标. 2. 实验设备与数据: 设 ...
Sentinel-2（哨兵-2）L1C数据辐亮度（辐射定标）和TOA反射率的获取说明
哨兵-2 L1C级数据是我们常用的一种数据,但是如和从L1C数据中获取需要的大气顶表观反射率(TOA)以及大气顶的辐亮度()呢? 获得反射率TOA一般就是用DN值除以10000(一万),这个值是固定的 ...
Linux系统下Sen2Cor对Sentinel哨兵2号遥感数据预处理(辐射定标和大气校正）Sen2Cor下载，使用
Sen2Cor插件的作用是对哨兵2号数据进行大气校正和辐射定标,将L1C数据,处理成为L2A数据.网上用window的教程有很多,这是linux的下载安装使用Sen2Sor的教程,我的系统是ubunt ...
IDL学习——处理自带经纬度文件的遥感影像，以哨兵5P数据为例
IDL学习--处理自带经纬度文件的遥感影像,以哨兵5P数据为例自带经纬度文件处理流程 1.寻找经纬度文件 2.构建GLT文件 3.使用GLT文件进行几何校正 4.批处理过程中遇到的问题最近一直在做 ...

IDL学习——哨兵2 L1C数据辐射定标

转换公式

参数查找

日地距离、太阳辐照度

xml文件

IDL读取

太阳高度角

xml文件

IDL读取

代码实现

IDL学习——哨兵2 L1C数据辐射定标相关推荐

最新文章

热门文章