简介这是一份“遥感Python绝对辐射定标”配套工程资源面向遥感、GIS及Python学习者解决如何利用GDAL库将遥感图像原始数字值转换为有物理意义的辐射值这一核心问题。压缩包共9个文件、约16KB包含一个Python定标脚本、一份JSON辐射定标参数文件、若干XML工程配置及Git忽略规则等结构精简便于对照教程逐行理解定标流程。内容以绝对辐射定标为主线覆盖数据准备、元数据与定标参数读取、各波段Band对象遍历、辐射校正值计算以及按需结合大气校正模型的处理思路代码框架清晰可快速改写成多波段批量处理工具或作为独立模块嵌入其他遥感流程。已有903人学习下载适合作为课程设计、科研预处理或GDAL二次开发的入门实操样板尤其有助于理解从DN值到反射率/辐射值的换算链条。1. 绝对辐射定标为什么DN转成物理量才是定量分析的前提做遥感定量分析的人多半被同一件事坑过两期影像在同一块农田上算NDVI趋势一会儿升一会儿降换了传感器之后所有结论都反了。原因在于DN不是一个物理量它只是传感器量化输出的一串整数增益、偏置、太阳高度角不同都会让这串数跑偏。遥感Python绝对辐射定标要做的就是把DN转换成辐亮度或表观反射率让不同时期、不同传感器、不同光照条件下的影像能在同一物理尺度上互相比较。这个操作看起来只是一次乘加和一次余弦改正但系数从哪里读、单位怎么统一、输出怎么保存每个环节都有细节。这篇笔记会沿着一条完整链路走读TIF、解析元数据、批量定标写回再把验证方法和踩坑记录一次讲清新老手都能直接照着改。2. 定标链路增益、偏置与太阳参数如何把DN变成表观反射率2.1 绝对辐射定标与相对辐射定标的分工遥感里说的“定标”其实有两层。相对辐射定标只是把不同影像的DN分布对齐比如做直方图匹配让两期影像看起来亮度一致但输出仍然是DN不是物理量。绝对辐射定标则必须把DN换算成传感器入瞳处的辐亮度或者进一步换算成表观反射率。这个区别在做时序分析和定量反演时非常关键NDVI差值分析用相对定标勉强能过但要反演叶绿素浓度、水体悬浮物含量输入必须是辐亮度或反射率这类有单位的物理量否则经验模型里的系数全部错位。严格说起来绝对辐射定标只是整个辐射校正链路的前半段。DN经乘加运算得到辐亮度再结合太阳角度和日地距离得到表观反射率这一层一般不涉及大气。真正的地表反射率还需要大气校正把大气散射和吸收的影响去掉。我在处理业务数据时会把DN到表观反射率这一段统称为绝对辐射定标把大气校正作为单独一步处理这样职责清楚也方便排查问题。下面所有代码默认输出到表观反射率这一层。2.2 直读式定标系数MTL与XML里的Mult和Add现在主流光学遥感卫星的L1级产品官方元数据文件里会直接给出定标系数用户不需要自己去查传感器的绝对定标参数表。以常见的30米多光谱数据为例它的元数据文件通常叫MTL.txt里面有一组字段长这样RADIANCE_MULT_BAND_1 0.0211 RADIANCE_ADD_BAND_1 -5.2311 REFLECTANCE_MULT_BAND_1 0.000023 REFLECTANCE_ADD_BAND_1 -0.099 SUN_ELEVATION 56.52细看会发现这一份文件里同时给了两套系数一套是辐亮度定标系数公式是L ML * DN AL另一套是表观反射率定标系数公式是ρ (Mρ * DN Aρ) / sin(θe)其中 θe 是太阳高度角。两套系数不能混用ML与Mρ相差上千倍代入错的系数结果立刻乱套。另有一类10米级多光谱数据的L1C产品处理方式更简单量化值DN除以10000就是表观反射率不需要逐波段查找增益和偏置。但这条规则只适用于该产品的L1C级别如果拿到的是L2A地表反射率产品那已经是大气校正后的结果再做一次“定标”反而画蛇添足。所以动手之前先确认产品级别这比选Python库还重要。对于一些科学卫星数据定标系数不是写在文本里而是存在HDF文件属性中。比如MODIS的L1B辐亮度产品每个波段有两个属性一个叫radiance_scales一个叫radiance_offsets读取方式变成查属性表而不是查文本。HDF属性读取用Python的h5py或pyhdf也很快几十行能搞定但这里先不展开后面提到的代码逻辑是一样的。2.3 传统公式里的日地距离与ESUN两种反射率算法对比如果拿到的是旧一代卫星数据元数据里没有直读式反射率系数就需要走经典公式计算表观反射率ρ π * L * d² / (ESUN * cos θz)其中L是辐亮度d是日地距离修正因子ESUN是波段平均太阳辐照度θz是太阳天顶角。d不是一个天文单位距离的绝对值而是一个以年为周期的比例系数常见近似公式是d 1 - 0.01673 * cos(0.0172 * (DOY - 2))DOY是影像获取日期的年积日1月1日为第1天。这个近似公式的精度已经足够误差在0.1%以内完全不影响后续反演。ESUN则是各波段的太阳等效辐照度在数据产品的官方文档里能查到不同卫星、不同波段的ESUN都不一样不能跨数据源混用。把直读式和经典式放在一起对比各自的适用情况就很清楚了。对比项直读式反射率系数经典式系数来源元数据文件直接给出官方辐射定标参数表需要的太阳参数SUN_ELEVATION太阳天顶角或高度角需要的辅助参数无日地距离、ESUN适用产品近十年主流光学卫星老一代卫星或无反射率系数的数据出错概率较低但要小心波段对应较高参数多了容易张冠李戴我一般优先用直读式系数因为它少传两个参数少两个出错点。经典式保留在代码里做兜底处理老数据时再用。3. 用Python跑通最小定标流程读DN、算辐亮度、写回GeoTIFF3.1 最小命令rasterio读DN并完成辐射定标先给一个能直接跑的最小示例。假设你已经从影像产品中拿到了一个波段TIF和对应的定标系数第一步就是把DN读进来做乘加运算。import rasterio import numpy as np with rasterio.open(scene_band1.tif) as src: dn src.read(1).astype(float64) profile src.profile gain 0.0211 # 示例值必须替换成元数据里的RADIANCE_MULT_BAND_x offset -5.2311 # 示例值必须替换成元数据里的RADIANCE_ADD_BAND_x radiance gain * dn offset print(radiance range:, radiance.min(), radiance.max())这段代码取的是单波段。点开TIF时注意波段序号的起始位置rasterio的read(1)读的是第一个波段对应元数据里的BAND_1。如果影像的波段序号从0开始的产品就必须先确认波段映射关系否则定标完的波段和实际物理波段对不上。这是最容易在批量处理时翻车的细节。把DN转成float64是刻意的。很多光学卫星影像是uint16存储DN范围能到0到65535增益乘进去之后数值可能跨好几个数量级用uint16中间存储会溢出用float32虽然精度够但在极端高值和低值像元混合的场景下有舍入风险。定标是逐像元线性变换我用float64跑一遍不心疼反正后面还会转回float32写出。3.2 自动解析MTL文件别再手动抄系数手动把gain和offset填进代码能跑通一个波段但处理一景影像十几个波段时效率太低而且很容易把BAND_5的系数贴到BAND_6上。更好的办法是直接从MTL文件里读取系数。import re def parse_mtl(mtl_path): params {} with open(mtl_path, r, encodingutf-8) as f: for line in f: line line.strip() if not in line: continue key, value line.split(, 1) params[key.strip()] value.strip() return params mtl parse_mtl(scene_MTL.txt) band 1 refl_mult float(mtl[fREFLECTANCE_MULT_BAND_{band}]) refl_add float(mtl[fREFLECTANCE_ADD_BAND_{band}]) sun_elev float(mtl[SUN_ELEVATION]) # 直读式反射率公式 toa (refl_mult * dn refl_add) / np.sin(np.radians(sun_elev))这个parse_mtl函数看起来很简陋但对MTL这类“缩进由空格构成、分组用GROUP/END_GROUP包裹”的文本逐行按等号切分反而比正则更稳。正则匹配如果要处理GROUP嵌套、键名重复容易写出体积膨胀且难维护的模式。逐行split的唯一注意点是某些字段值可能带引号如果引号存在需要再处理一次def clean_value(v): v v.strip() if v.startswith() and v.endswith(): return v[1:-1] return v解析完成后mtl就是一个以字段名为键的字典接下来无论你是做单波段测试还是批量处理都可以从同一个字典取系数。取值时用f-string拼波段号循环里非常方便而且波段号写错会在第一次运行时直接报KeyError比静默读到错误系数好查得多。3.3 批量定标多波段并写回GeoTIFF保留地理信息还能压缩有了解析函数下一步就是批量处理一景影像的所有波段并把结果写成新的GeoTIFF。如果还是用老办法——先转成numpy数组计算再用rasterio新建文件写CRS和transform很容易把地理信息写丢。正确姿势是直接从输入文件的profile复制再修改必要字段。import rasterio import numpy as np def calibrate_toa_band(src_path, dst_path, refl_mult, refl_add, sun_elev): with rasterio.open(src_path) as src: dn src.read(1).astype(float32) profile src.profile.copy() toa (refl_mult * dn refl_add) / np.sin(np.radians(sun_elev)) profile.update( dtypefloat32, nodata-9999.0, compressdeflate, predictor3, ) with rasterio.open(dst_path, w, **profile) as dst: dst.write(toa, 1)逐行说三个参数。nodata设置成-9999是因为定标后有效反射率理论上在-0.1到1.2之间-9999足够远离有效范围后续统计时一眼能识别。compressdeflate是GeoTIFF里通用性最好的压缩方式不要用LZWdefrate在大多数平台兼容性更好而且压缩比相当。predictor3是浮点数据专用的预测器它按浮点数的指数部分做差分预测对反射率这种取值相对连续的浮点数组文件体积能比不设预测器再小20%到40%。批量处理时我会建一个波段映射表把输入文件名列表和波段序号对齐band_files [ scene_B2.TIF, scene_B3.TIF, scene_B4.TIF, scene_B5.TIF, ] for idx, src_path in enumerate(band_files, start1): band_idx idx 1 # 根据产品实际波段号调整 mult float(mtl[fREFLECTANCE_MULT_BAND_{band_idx}]) add float(mtl[fREFLECTANCE_ADD_BAND_{band_idx}]) dst_path fTOA_B{band_idx}.TIF calibrate_toa_band(src_path, dst_path, mult, add, sun_elev) print(done:, dst_path)这里最值得注意的就是波段映射。不同产品的文件名后缀与实际波段号的关系并不固定有的用B2、B3、B4有的用Band1、Band2。批量循环中如果不对应程序不会报错只会默默产出错误的反射率图。所以我常在循环里加一个断言assert refl_mult 0 and abs(refl_add) 1.0, 系数疑似读错多光谱反射率定标系数的增益通常是个很小的正值偏置通常在零附近。如果读出来的是负数或上百十有八九是波段号对应错了。4. 定标路上翻车最多的地方5个常见问题排查4.1 定标结果全黑是失败还是显示设置问题现象定标完的浮点影像在GIS软件里打开整幅图是全黑的看不到任何地物轮廓。原因反射率浮点影像的数值集中在0到0.5之间而大多数软件默认按全范围线性拉伸显示。如果影像里有一个异常的极大值比如云顶反射率接近1或者有一个负的无效值-9999显示系统会把有效值压缩到极窄的灰度区间看起来就是全黑。解决先别怀疑定标失败用统计量确认数值分布。with rasterio.open(TOA_B4.TIF) as src: arr src.read(1) mask arr -1000 # 过滤nodata valid arr[~mask] print(np.nanpercentile(valid, [2, 50, 98]))如果2%分位数和98%分位数落在0到0.6之间定标结果基本正常。显示时手动拉伸到这两个分位数或者用软件的2%线性拉伸就能正常看。我习惯把定标代码里顺手输出一组统计量处理完第一眼就能判断有没有低级错误。4.2 辐亮度大了1000倍先检查单位是μm还是nm现象定标出来的辐亮度数值比经验估算大了1000倍或者某个波段的辐射值明显不合理。原因辐亮度的单位在不同数据产品里不完全一致。光学遥感常用的单位是W/(m²·sr·μm)但部分传感器产品在文档里给出的是W/(m²·sr·nm)两者差了1000倍。如果按μm的单位读入数据却套用nm系数整幅图所有像元都会放大1000倍。这个错误不报错、不影响显示拉伸但下游反演模型会彻底崩溃。解决拿到数据先看波段中心波长的量级。中心波长在几百nm量级的传感器辐亮度如果写成nm单位数值会非常小接近0.001级别μm单位的数值则在1到10之间。根据量级判断该用哪套系数。更稳妥的做法是定标后输出一个中位数水体在红波段辐亮度应该在个位数水平如果中位数变成几千几乎可以断定单位错位。4.3 天顶角与高度角混用反射率直接超过1现象计算表观反射率后大量像元值超过1甚至出现几十倍的极端值。原因太阳参数有两种表示方式太阳天顶角是垂直方向与太阳光线的夹角太阳高度角是水平面与太阳光线的夹角两者互余差90度。很多元数据文件里写的是SUN_ELEVATION也就是高度角但计算反射率时有人会把它当成天顶角用余弦去算。天顶角50度对应的余弦是0.64高度角50度对应的正弦是0.77混用之后反射率会被放大或缩小约20%。更极端的案例是把高度角直接代入余弦天顶角代入正弦结果反射率超过2。解决明确记住一句口诀——元数据给的是高度角就用正弦给的是天顶角就用余弦。下面的代码同时支持两种模式def to_toa_reflectance(dn, refl_mult, refl_add, sun_elevationNone, sun_zenithNone): if sun_elevation is not None: denominator np.sin(np.radians(sun_elevation)) elif sun_zenith is not None: denominator np.cos(np.radians(sun_zenith)) else: raise ValueError(需要提供太阳高度角或天顶角) return (refl_mult * dn refl_add) / denominator额外提醒角度转弧度这步不能省numpy的三角函数默认按弧度计算直接传入50度结果会非常荒谬。np.radians就是用来做这个转的很多人漏掉它。4.4 多景影像共用一份定标系数时序分析数据全部作废现象批量处理一个月内的多景影像时个别景的NDVI时序出现突变但原始影像直方图看起来正常。原因光学卫星的定标系数不是永久不变的。卫星发射后会定期更新绝对定标参数同一天的L1产品可能因为元数据版本更新增益系数发生微小变化。如果整个批次处理时只从某一景的MTL里读了一次系数全批次共用了那其他几景影像都会带上系统偏差。偏差小的时候看不出来一旦做NDVI阈值分类或变化检测边界就出来了。解决每一景影像单独解析它自己的MTL文件哪怕来自同一目录。批量处理时脚本里的循环应该是“一景一读”而不是“一批共读一份”import pathlib scene_list pathlib.Path(scenes).glob(*_MTL.txt) for mtl_path in scene_list: mtl parse_mtl(str(mtl_path)) sun_elev float(mtl[SUN_ELEVATION]) # 再按波段处理该景的TIF文件前缀与MTL同名我一般还会把处理日期和MTL文件路径记录到日志文件里这样复查数据的时候能快速定位是哪一份MTL产生了问题。4.5 写入float影像时把nodata当成有效值现象定标后的影像在查询统计量时最小值总是-9999或者融合分析时边缘出现大片异常低值。原因原始影像的边界外像元本身是无效值DN可能是0或者某个填充值。定标公式只做乘加和除法不会自动识别这些无效值结果边缘外像元也参与了运算变成负值或极大值。解决把掩膜信息从输入影像带出来在输出前标记无效值。with rasterio.open(src_path) as src: dn src.read(1).astype(float32) mask src.read_masks(1) # 有效像元为255无效为0 profile src.profile.copy() toa (mult * dn add) / np.sin(np.radians(sun_elev)) toa np.where(mask 0, -9999.0, toa) profile.update(dtypefloat32, nodata-9999.0, compressdeflate, predictor3) with rasterio.open(dst_path, w, **profile) as dst: dst.write(toa, 1)read_masks(1)读取的是波段的有效区域掩膜返回0到255的数组只有像元有效时是255。np.where把无效位置替换成nodata。这是最容易忽略的一步尤其当原始影像四周有黑边时不定标看不出来定标后黑边变成巨大的负值会让后续所有统计和可视化都受到污染。5. 定标精度验证没有地面真值也能做的三个检查5.1 用官方示例像元回算验证公式与系数无误大多数公开遥感产品的文档里会给出至少一个示例像元的DN值以及对应的辐亮度或反射率。这是一种非常廉价的验证方式。先把你的定标函数写成一个纯函数然后把官方示例DN代入def dn_to_toa(dn, refl_mult, refl_add, sun_elev): return (refl_mult * dn refl_add) / np.sin(np.radians(sun_elev)) official_dn 8200 official_toa 0.2513 # 假设文档给出的参考值 my_toa dn_to_toa(official_dn, refl_mult, refl_add, sun_elev) print(relative error:, abs(my_toa - official_toa) / official_toa * 100)相对误差在0.1%以内说明公式和系数都没问题。如果误差超过1%先查有没有把Mρ和ML写反再查太阳角度单位。这个验证步骤最好写进批处理脚本的开头每次换数据源先跑一遍能省掉后面排查的半天时间。5.2 用典型地物的反射率范围做定性校验没有官方示例也没关系遥感图像灰度直方图里隐含了地物信息。定标后的表观反射率有个经验范围清洁深水体在近红外波段通常低于0.05植被在近红外波段通常高于0.3裸土在红波段大约0.15到0.3云和雪可以冲到0.8以上。拿到定标影像后我习惯直接统计几个已知地物点的反射率比如一片湖泊、一片农田、一个裸土区域import rasterio points { water: (420000, 3260000), crop: (421000, 3261000), bare: (422000, 3262000), } with rasterio.open(TOA_B4.TIF) as src: for name, (x, y) in points.items(): row, col src.index(x, y) value src.read(1, window((row, row 1), (col, col 1)))[0, 0] print(name, value)近红外波段水体反射率如果超过0.1或者裸土接近0说明定标有问题。用这种方式做不了精确验证但能在十分钟内筛掉80%的低级错误。注意这些阈值适用于表观反射率如果做了大气校正得到地表反射率水体值会更低裸土上限也有变化。5.3 多期影像一致性检验定位定标异常时序分析场景下还有一个更实用的验证逻辑同一地物在多期影像里的表观反射率虽然会有波动但不会出现成倍变化。如果同一块农田在A影像里近红外反射率是0.35B影像里变成0.8排除极端天气和作物收割因素后第一个怀疑对象就是定标系数。def check_consistency(toa_a, toa_b, thresh0.5): ratio toa_b / (toa_a 1e-6) bad np.abs(ratio - 1.0) thresh return np.mean(bad) * 100 # 返回异常像元比例这个检验函数适合在两期影像覆盖重叠区域时跑异常比例超过5%就要回去查太阳角度或者系数是否取自正确的波段。我自己处理几十年土地变化数据时这套方法抓出过两期MTL的SUN_ELEVATION字段读取顺序错误的问题比肉眼盯直方图快得多。6. 把这套流程封装成命令行工具让下一次定标不再重新踩坑我习惯把解析MTL、定标、写GeoTIFF封装进一个独立脚本用命令行参数驱动。这样每次处理新数据不用再打开Jupyter手动改路径批处理和多景并行也能复用同一套逻辑。一个精简版的结构如下# toa_calib_cli.py import argparse import json import pathlib import rasterio import numpy as np def parse_mtl(mtl_path): params {} with open(mtl_path, r, encodingutf-8) as f: for line in f: line line.strip() if not in line: continue key, value line.split(, 1) params[key.strip()] value.strip() return params def calibrate_file(src_path, dst_path, mult, add, sun_elev): with rasterio.open(src_path) as src: dn src.read(1).astype(float32) mask src.read_masks(1) profile src.profile.copy() toa (mult * dn add) / np.sin(np.radians(sun_elev)) toa np.where(mask 0, -9999.0, toa) profile.update(dtypefloat32, nodata-9999.0, compressdeflate, predictor3) with rasterio.open(dst_path, w, **profile) as dst: dst.write(toa, 1) def main(): parser argparse.ArgumentParser(descriptionTOA绝对辐射定标工具) parser.add_argument(--mtl, requiredTrue, helpMTL/XML元数据文件路径) parser.add_argument(--band, typeint, requiredTrue, help波段号) parser.add_argument(--src, requiredTrue, help输入TIF路径) parser.add_argument(--dst, requiredTrue, help输出TIF路径) args parser.parse_args() mtl parse_mtl(args.mtl) mult float(mtl[fREFLECTANCE_MULT_BAND_{args.band}]) add float(mtl[fREFLECTANCE_ADD_BAND_{args.band}]) sun_elev float(mtl[SUN_ELEVATION]) calibrate_file(args.src, args.dst, mult, add, sun_elev) with open(pathlib.Path(args.dst).with_suffix(.json), w, encodingutf-8) as f: json.dump({ src: args.src, band: args.band, mult: mult, add: add, sun_elev: sun_elev, }, f, indent2, ensure_asciiFalse) if __name__ __main__: main()命令行运行方式python toa_calib_cli.py --mtl scene_MTL.txt --band 4 --src scene_B4.TIF --dst TOA_B4.TIF脚本会在输出TIF旁边生成一个同名JSON记录这次定标用的系数和太阳高度角。此举不是为了好看是为了下次查数据时不用再翻原始MTL——看一眼JSON就知道某幅影像的数值是怎么来的。我在这套小工具上修修补补了两年翻车最多的地方从来不是numpy而是系数有没有取自对应波段、太阳角度有没有搞错。现在把关键参数落在JSON里每次处理完顺手核对一眼基本没再出过低级错误。希望帮到你。本文还有配套的精品资源点击获取