简介面向遥感图像处理初学者与需要落地绝对辐射定标的开发者这份压缩包围绕Python与GDAL展开聚焦把影像原始数字值转换为具有物理意义的辐射值这一关键环节。包内共9个文件以xml、json和py为主压缩后约16KBPython脚本承载波段读取与定标计算JSON用于组织辐射校正参数XML配合元数据与工程配置另有少量工程辅助文件整体轻量便于查阅与复用。已有903人学习说明该主题在遥感数据处理学习中有一定关注度。内容覆盖定标参数读取、逐波段线性校正、结果写回与元数据更新等流程并提及大气校正的衔接思路可作为GDAL读写、增益与暗电流处理及排错思路的参考模板适合希望快速理解绝对辐射定标实现脉络的读者。1. 一景影像对不上往往不是大气的问题同一片农田两个时相的影像摆在一起植被指数差了 0.15第一反应是大气影响或者物候变化把 DN 值拉出来一看才发现两景的量纲根本不在一个水平线上。绝对辐射定标要处理的就是这件事把传感器记录的原始 DN 值还原成有物理单位的辐亮度再进一步换算成大气顶表观反射率让不同时间、不同传感器、不同成像条件下的数据具备可比性。这份资源给的是 python gdal 遥感场景下的一套完整实现从元数据解析、系数管理到分块读写、命令行调用都串了起来。适合刚接触遥感影像预处理、需要把 L1 产品推到 L2 的分析人员也适合手里已经有一堆脚本、想把定标环节做成稳定模块的开发者。它不解决大气校正但对定标这一步该踩的坑交代得比较全。2. 项目骨架与数据流GDAL 负责 IOnumpy 负责算元数据负责参数2.1 目录结构与模块职责这套代码没有做成大包大揽的框架按职责切成四层入口层、元数据层、算法层、IO 层。入口层只做参数解析和任务调度元数据层负责把产品说明书翻译成程序能用的系数对象算法层是纯函数、只接 numpy 数组IO 层封装 GDAL 的打开、建文件、分块读和分块写。这么切的原因很实际——算法层要能被单元测试直接调用不需要真的准备一景影像IO 层要能替换比如以后换成别的库读写算法一行不用动。路径职责关键点calib/cli.py命令行入口参数校验、任务分发只依赖 argparse不放业务逻辑calib/meta.py解析元数据文件产出SceneMeta兼容文本键值对与 XML 两种格式calib/spec.py数据类定义BandCalib/SceneMeta系数与场景元数据分离calib/algo.pydn_to_radiance、radiance_to_reflectance纯 numpy无 IO 依赖calib/raster.py分块读写、输出创建、地理参考复制全部 GDAL 调用集中在此calib/runner.py编排读元数据 → 遍历波段 → 分块处理唯一知道全流程的地方tests/系数解析与算法的小样本测试用手工构造的 3×3 数组这种切法的直接好处是排错路径清晰。结果量级不对先看meta.py解析出的系数结果数组形状错乱去看raster.py的窗口计算中间某个值算错单独调algo.py里的函数传两个数进去就能复现。血泪经验是定标这类数值处理最怕把所有东西揉在一个三百行的脚本里改一个系数要重跑半小时才知道对不对。2.2 定标系数的来源与解析系数不是拍脑袋来的公开的 L1 级产品一般在元数据里带两类字段一类是辐亮度定标系数命名形如RADIANCE_MULT_BAND_n和RADIANCE_ADD_BAND_n另一类是反射率定标系数REFLECTANCE_MULT_BAND_n和REFLECTANCE_ADD_BAND_n。太阳高度角和日地距离也在同一份文件里。解析器要做的事情就是把这些散落的键值对收拢成一个结构体顺手做单位归一。# calib/meta.py import re import xml.etree.ElementTree as ET from pathlib import Path # 元数据里常见的几种键名统一映射到内部字段 KEY_PATTERNS { radiance_mult: re.compile(rRADIANCE_MULT_BAND_(\d)$, re.I), radiance_add: re.compile(rRADIANCE_ADD_BAND_(\d)$, re.I), refl_mult: re.compile(rREFLECTANCE_MULT_BAND_(\d)$, re.I), refl_add: re.compile(rREFLECTANCE_ADD_BAND_(\d)$, re.I), } def parse_kv_text(path: Path) - dict: 解析 MTL 风格的 KEY VALUE 文本忽略注释行和空行 raw {} for line in path.read_text(encodingutf-8, errorsignore).splitlines(): line line.strip() if not line or line.startswith(#) or not in line: continue k, _, v line.partition() raw[k.strip().upper()] v.strip().strip() return raw def build_scene_meta(raw: dict) - SceneMeta: meta SceneMeta( sun_elevationfloat(raw.get(SUN_ELEVATION, 90.0)), earth_sun_distancefloat(raw.get(EARTH_SUN_DISTANCE, 1.0)), quantize_valuefloat(raw.get(QUANTIZE_VALUE, 1.0)), ) for field_name, pat in KEY_PATTERNS.items(): for key, val in raw.items(): m pat.match(key) if not m: continue band int(m.group(1)) calib meta.bands.setdefault(band, BandCalib()) # 系数按字段名落到对应属性缺哪项补默认值 attr { radiance_mult: radiance_mult, radiance_add: radiance_add, refl_mult: reflectance_mult, refl_add: reflectance_add, }[field_name] setattr(calib, attr, float(val)) return metaparse_kv_text只做字符串层面的清洗把等号两边拆开、去掉引号、键名统一大写不做任何数值假设。build_scene_meta用一组正则把带波段号的键归到对应波段对象上这样 7 个波段也不用手写 7 段代码。quantize_value的默认值给 1.0是因为不少产品的 DN 就是量化后的整数不需要再缩放但确实有产品要乘一个比例留个口子比事后改算法省事。注意第 5 章会讲这个字段用错会直接导致结果差三个数量级。XML 格式的元数据用ElementTree走一遍就行难点不在解析在于不同产品的标签层级不一样有时是平铺的band number1有时嵌在radiance下面。我的做法是先iter()全树按 tag 名匹配不写死路径:def parse_xml(path: Path) - dict: tree ET.parse(path) raw {} for node in tree.iter(): tag node.tag.split(})[-1].upper() # 去掉可能的命名空间 if node.text and node.text.strip(): raw[tag] node.text.strip() # 波段级字段常见写法属性 number 决定波段号 num node.attrib.get(number) if num and node.text: raw[f{tag}_BAND_{num}] node.text.strip() return raw这段的关键是tag.split(})[-1]去掉命名空间前缀否则带ns0:的标签一个都匹配不上这类问题在解析行业标准格式的 XML 时特别常见。2.3 数据流与内存策略为什么必须分块完整一景多光谱影像单波段动辄上万个像素的边长uint16 存储时几百 MB转换成 float32 就是翻倍。如果走整幅 ReadAsArray → 算 → 整幅 WriteArray的路子峰值内存是原文件的 3 到 4 倍机器上同时跑几个任务直接会被系统干掉。分块不是优化项是能不能跑完的前提。实现上的取舍是把块大小做成参数默认 1024×1024。太小了 GDAL 的调用开销占主导太大了省不下内存。还有一个容易被忽略的点分块窗口要和输出文件的分块设置对齐否则写入时 GDAL 内部还要做一次重排速度会掉得很明显。所以建输出文件时统一设TILEDYES块大小和读取块保持一致。参数默认值影响block_size1024单次读写的像素边长内存占用近似block² × 4 × 波段数TILEDYES输出文件按块组织与分块写入对齐BIGTIFFIF_NEEDED超过 4GB 自动切大文件格式COMPRESSDEFLATE浮点数据压缩比一般但能省一半磁盘分块之后还要处理边界影像宽高不一定是块大小的整数倍最后一块的窗口要min(block, width - xoff)截一下否则 GDAL 会抛越界错误。这个逻辑在raster.py里统一封装上层的runner.py只管传窗口坐标不用每次自己算。3. 核心算法拆解从 DN 到辐亮度再到反射率的完整链路3.1 辐亮度定标乘加两个系数就够了吗辐亮度的公式本身很朴素L M_L × Q_cal A_LM_L是乘系数A_L是加系数Q_cal是量化后的 DN。看起来一行就能写完实际落地时有三个细节决定结果对不对。第一是Q_cal。有些产品的 DN 本身就是量化值直接代进去有些产品的 DN 需要先乘QUANTIZE_VALUE常见是 0.01 或 1.0这一步漏了整体就差 100 倍。第二是加系数的量级A_L一般是负值或者零用来修正暗电流别看它小在低反射地物上占比不小。第三是数据类型dn * mult时如果 dn 还是 uint16numpy 会按整型运算小数部分直接被截掉结果全错。必须先把数组转成 float32 再算。# calib/algo.py import numpy as np def dn_to_radiance(dn: np.ndarray, mult: float, add: float, quantize: float 1.0) - np.ndarray: dn: 原始数字量化值数组 mult: 辐亮度乘系数单位 W/(m2*sr*um) add: 辐亮度加系数同单位 quantize: DN 的量化比例多数产品为 1.0 返回 float32 辐亮度数组 # 关键先转 float32避免整型运算截断小数 arr dn.astype(np.float32) if quantize ! 1.0: arr arr * np.float32(quantize) out arr * np.float32(mult) np.float32(add) return out.astype(np.float32, copyFalse)astype(np.float32)放第一行是刻意的这条语句是整个函数最容易被忽略也最容易翻车的地方。后面所有常量都显式包一层np.float32是为了避免 numpy 把 float32 数组和 python float 运算时悄悄提升成 float64内存翻倍不说大数据量下速度也掉。copyFalse让返回时不再复制一次能省一块内存峰值。3.2 表观反射率太阳高度角和日地距离两个修正项反射率比辐亮度多用两个量太阳高度角和日地距离。公式是ρ (π × L × d²) / (ESUN × sin(h))L是上一步算出来的辐亮度d是日地距离天文单位ESUN是该波段的大气顶太阳辐照度h是太阳高度角。注意分母里是sin(h)不是cos(h)——因为sin(高度角) cos(天顶角)很多实现把天顶角代进去用 cos写成cos(h)就等价于把太阳当成了在天顶低纬度冬季成像的影像会整体偏亮。另一个坑是角度单位。元数据里的太阳高度角通常是度numpy 的三角函数要弧度转换这步必须显式写np.deg2rad一行的事漏了就差出几十倍的反射率。日地距离如果元数据里没给可以用日期算个近似值误差在千分之一量级对反射率影响不大。def radiance_to_reflectance(radiance: np.ndarray, esun: float, sun_elevation_deg: float, earth_sun_distance: float 1.0) - np.ndarray: esun: 大气顶太阳辐照度W/(m2*um)按波段取值 sun_elevation_deg: 太阳高度角单位度来自元数据 earth_sun_distance: 日地距离天文单位默认 1.0 # 角度转弧度漏了这一步结果会差几十倍 elev_rad np.deg2rad(np.float32(sun_elevation_deg)) sin_elev np.sin(elev_rad) # 极区冬季太阳高度角可能为负钳一下避免除零或负值 sin_elev np.maximum(sin_elev, np.float32(0.01)) d2 np.float32(earth_sun_distance) ** 2 coef np.float32(np.pi) * d2 / (np.float32(esun) * sin_elev) return (radiance * coef).astype(np.float32, copyFalse)sin_elev的钳位不是洁癖。高纬度地区冬季成像时太阳高度角可能接近 0 甚至为负直接除会得到 inf 或者符号翻转的反射率后续做指数计算全是 NaN。给一个 0.01 的下限结果虽然不可信但至少不会污染整幅影像的统计。ESUN按波段查表一般在 1500 到 2000 之间不同波段差异不小不能全波段共用一个值。如果产品元数据直接给了反射率系数那更省事ρ M_ρ × Q_cal A_ρ不用自己算太阳角度和日地距离这两条路径在代码里做成两个策略由--level参数切换。3.3 分块读写的落地写法算法是纯函数真正把它接到文件上的地方在raster.py。这里要做四件事打开源文件、按窗口读数组、调算法、写进输出文件的对应窗口位置。输出文件的地理变换、投影、NoData 必须在创建时就设好否则后面补很麻烦。# calib/raster.py from osgeo import gdal import numpy as np def open_src(path: str) - gdal.Dataset: gdal.UseExceptions() # 出错直接抛异常比返回 None 好排查 return gdal.Open(path, gdal.GA_ReadOnly) def create_dst(path: str, src: gdal.Dataset, dtypegdal.GDT_Float32) - gdal.Dataset: driver gdal.GetDriverByName(GTiff) dst driver.Create( path, src.RasterXSize, src.RasterYSize, src.RasterCount, dtype, options[ TILEDYES, BLOCKXSIZE1024, BLOCKYSIZE1024, COMPRESSDEFLATE, PREDICTOR3, # PREDICTOR3 对浮点有效 BIGTIFFIF_NEEDED, ], ) # 地理参考必须复制否则输出是一张没有坐标的图 dst.SetGeoTransform(src.GetGeoTransform()) dst.SetProjection(src.GetProjection()) return dst def iter_windows(width: int, height: int, block: int 1024): 生成不越界的窗口 (xoff, yoff, xsize, ysize) for y in range(0, height, block): h min(block, height - y) for x in range(0, width, block): w min(block, width - x) yield x, y, w, hgdal.UseExceptions()是这几行里最值钱的一句。默认情况下 GDAL 出错只往 stderr 打一行警告然后返回 None脚本继续往下跑最后得到一个空文件或者半截结果。开启异常之后文件打不开、磁盘满了、权限不对都会当场抛出来。PREDICTOR3专门给浮点数据用压缩比能比默认高不少配DEFLATE是最稳的组合。iter_windows返回的是生成器大影像上万个块也不会一次性把坐标列表铺在内存里。min那两行处理的就是最后一行、最后一列的残缺块。把这几块串起来的主循环在runner.pydef run(src_path: str, dst_path: str, meta, levelreflectance, block1024): src open_src(src_path) dst create_dst(dst_path, src) for xoff, yoff, xsize, ysize in iter_windows(src.RasterXSize, src.RasterYSize, block): for b in range(1, src.RasterCount 1): band_calib meta.bands.get(b) if band_calib is None: raise KeyError(f波段 {b} 缺少定标系数检查元数据波段号是否对齐) dn src.GetRasterBand(b).ReadAsArray(xoff, yoff, xsize, ysize) rad dn_to_radiance(dn, band_calib.radiance_mult, band_calib.radiance_add, meta.quantize_value) if level reflectance: out radiance_to_reflectance( rad, band_calib.esun, meta.sun_elevation, meta.earth_sun_distance) else: out rad dst.GetRasterBand(b).WriteArray(out, xoff, yoff) # 每处理完一块刷一次缓存崩了不至于全丢 dst.FlushCache() dst None # 显式释放确保 TIFF 目录写入磁盘band_calib is None的检查看着啰嗦实际很有用产品元数据里的波段编号和影像的波段顺序不一定一致某个波段缺系数时静默跳过会得到一张部分波段有值、部分波段全零的影像肉眼看不出来统计时才发现。直接抛出来跑一次就知道哪对不上。dst None这行必须写。GDAL 的 Python 绑定依赖引用计数来落盘函数返回前不释放文件句柄挂着后续再读这个文件可能读到不完整的头信息。我见过有人因为少了这一句在 Notebook 里反复重跑每次都读到上一次的残留结果。4. 跑起来环境、参数与一次完整执行4.1 环境与依赖依赖很少GDAL 和 numpy 两个测试用 pytest。GDAL 的安装是老生常谈的翻车点pip 直接装gdal经常编译失败原因是它要和系统里的 GDAL 库版本严格匹配。稳妥的做法是走 condaconda install -c conda-forge gdal numpy版本由 conda 统一解决。如果只能用 pip先装GDAL对应的 wheel或者用rasterio自带的 GDAL但这套代码直接用 osgeo不适用。# 推荐conda 一把装齐避免 GDAL 版本地狱 conda create -n radcalib python3.10 -y conda activate radcalib conda install -c conda-forge gdal numpy pytest -y # 验证 GDAL 是否可用以及版本 python -c from osgeo import gdal; print(gdal.__version__)Python 选 3.10 是因为数据类的写法在 3.7 之后才稳定3.10 是 conda-forge 上 GDAL wheel 覆盖比较全的版本。gdal.__version__打出来之后顺手确认一下和 conda 里gdalinfo --version是否一致不一致说明环境里混进了两套 GDAL后面会出各种莫名其妙的符号错误。4.2 命令行参数与调用方式入口设计成两种用法命令行跑单景Python API 跑批。命令行参数不多但每个都有默认值兜底。参数必填默认说明--src是无源影像路径GeoTIFF--meta是无元数据文件文本或 XML--dst是无输出路径--level否reflectance可选radiance/reflectance--block否1024分块边长内存紧张时调小--dtype否Float32输出数据类型别改成整型# 输出表观反射率 python -m calib.cli \ --src ./data/scene_001.tif \ --meta ./data/scene_001_MTL.txt \ --dst ./out/scene_001_refl.tif \ --level reflectance \ --block 1024 # 只做辐亮度给后面的大气校正用 python -m calib.cli \ --src ./data/scene_001.tif \ --meta ./data/scene_001_MTL.txt \ --dst ./out/scene_001_rad.tif \ --level radiance批处理走 API 更顺手写个循环遍历目录把元数据文件按同名规则匹配上去from pathlib import Path from calib.meta import parse_kv_text, build_scene_meta from calib.runner import run for tif in Path(./data).glob(*.tif): meta_file tif.with_name(tif.stem _MTL.txt) if not meta_file.exists(): print(f跳过 {tif.name}找不到元数据) continue meta build_scene_meta(parse_kv_text(meta_file)) run(str(tif), f./out/{tif.stem}_refl.tif, meta, levelreflectance)批处理里最容易出事的是元数据缺失continue之前打一行日志跑完看一眼就知道漏了哪几景比默默跳过强。4.3 结果自检不看统计值等于没跑跑完之后必须验验的方法不是打开看图是看统计数据。反射率的合理范围在 0 到 1 之间水体、浓密植被偏低裸土、云偏高极值偶尔能到 1.2 左右大气顶反射率允许超 1。如果最小值是 -32768 或者最大值是 65535说明 NoData 或者数据类型出了问题。# 逐波段统计输出 min / max / mean / stddev gdalinfo -stats ./out/scene_001_refl.tif | grep -A 6 Band 1-stats会生成一个.aux.xml缓存文件第二次读就快了。看几个关键信号min 接近 0 而 max 在 0.8 左右正常min 是负数且绝对值很大多半是 NoData 没处理max 只到 0.02大概率系数单位错了乘系数本应是 1e-2 量级却当成了 1e-4 又或者重复缩放了。标准差特别小比如 0.001也要警惕可能整幅影像被同一个常数覆盖那是读窗口的偏移写错了。5. 避坑与排查定标环节最常见的五个坑5.1 数据类型被截断结果全是整数现象输出影像的反射率只有 0 和 1 两个值或者辐亮度全是整数小数部分消失。原因源数据是 uint16dn * mult时 numpy 按整型规则运算乘出来再赋值给整型数组小数被截掉。更隐蔽的一种是输出文件用GDT_UInt16创建算法算的是浮点写进去还是被截断。解决算法入口第一行就astype(np.float32)输出统一GDT_Float32。改代码之前先gdalinfo看一眼输出的数据类型是UInt16就基本确定是这个原因。这类问题属于玄学范畴不算就发现不了一看就恍然大悟。5.2 NoData 被当成有效值参与运算现象影像边缘一圈或者整幅出现极端的负反射率统计最小值是 -0.3 甚至更低。原因源文件的填充值是 0 或某个哨兵值读进来直接乘加0 × mult add得到的是加系数本身如果加系数是负的出来就是负数。填充值本身还会把整幅的均值拉偏。解决读块之后先做掩膜mask dn nodata算完把掩膜位置的输出设成输出文件的 NoData 值并在create_dst里SetNoDataValue。别指望用np.nan蒙混过去float32 的 NaN 写进 TIFF 后很多软件读出来是乱值。5.3 地理参考丢失输出成了没有坐标的图现象输出文件能在 QGIS 里打开但拖进去不定位或者提示图层没有坐标参考。原因driver.Create新建的数据集默认没有地理变换和投影必须从源文件复制。常见写法是先Create再处理忘了SetGeoTransform和SetProjection这两句。解决在create_dst里统一补上别分散到调用方。验证方法很简单gdalinfo看输出里有没有Origin和Coordinate System is两段缺一段就是没设。这个坑不痛但烦等下游做裁剪时才发现整批结果都要重跑。5.4 系数单位理解错整体差 100 倍现象反射率最大值只有 0.01 左右或者辐亮度比正常值大两个数量级。原因元数据里的乘系数通常以1e-2量级给出比如2.0000E-02有人看到科学计数法以为是需要额外除以 100或者把QUANTIZE_VALUE又乘了一遍。反过来把已经含单位的系数当成了无量纲比例也是常见误读。解决拿一个已知地物做锚点。深水体的反射率在 0.02 到 0.05健康植被近红外在 0.3 到 0.5。算完抽几个像素对一下量级不对立刻回去查系数单位和quantize。我一般的习惯是先把系数打印出来看一眼数量级比起跑完再回头查省时间。5.5 太阳高度角用错低太阳高度影像整体偏亮现象冬季或者高纬度影像的反射率明显偏高同一地物比夏季影像高出 30% 以上。原因两种写法都会导致这个结果——直接用cos(高度角)而不是sin(高度角)等于是把太阳放到了天顶或者角度没转弧度就丢进np.sin得到的是一个毫无意义的小数。第三种是元数据里给的是天顶角字段名看岔了又转了一次。解决把sun_elevation的取值打印出来正常范围在 10 到 90 度之间如果出现 0.3 这种数就是弧度如果出现 60 多但结果是偏亮就检查是不是用了 cos。radiance_to_reflectance里的np.deg2rad不要省宁可多写一行。6. 二次开发把定标做成可插拔的处理器6.1 抽出处理器接口接新传感器不用改主流程现在的runner.py里还硬编码着如果是 reflectance 就走这条分支传感器一多就会变成一堆 if。改造的方向是定义一个处理器协议每个传感器或者每种定标级别实现一个类主流程只负责调度。from typing import Protocol import numpy as np class BandProcessor(Protocol): name: str def prepare(self, meta) - None: ... def process(self, dn: np.ndarray, band: int) - np.ndarray: ... class ReflectanceProcessor: 表观反射率处理器直接套用产品给的反射率系数 name reflectance def prepare(self, meta): self.meta meta def process(self, dn, band): c self.meta.bands[band] arr dn.astype(np.float32) * np.float32(c.reflectance_mult) return (arr np.float32(c.reflectance_add)).astype(np.float32) class RadianceProcessor: 辐亮度处理器先定标再按太阳角度换算 name radiance_to_refl def prepare(self, meta): self.meta meta self.sin_elev max(np.sin(np.deg2rad(meta.sun_elevation)), 0.01) def process(self, dn, band): c self.meta.bands[band] rad dn.astype(np.float32) * np.float32(c.radiance_mult) np.float32(c.radiance_add) coef np.pi * self.meta.earth_sun_distance ** 2 / (c.esun * self.sin_elev) return (rad * coef).astype(np.float32)prepare里做一次性计算比如太阳角度只算一次process只接数组和波段号这样主循环里就是out processor.process(dn, b)一行。加一个新传感器实现一个类注册进去就行runner.py不用动。6.2 验证改造有没有改坏三次对比法改完之后不能只看能跑通要有可比对的基准。我的做法是保留一份改之前的输出改完在同一景上重跑做三次对比逐波段统计量min/max/mean/stddev差异在 1e-4 以内随机抽 1000 个像素做逐点差值最大值应为 0输出文件的地理变换、投影、NoData 三项与基准完全一致。前两项验数值第三项验元信息缺一项都可能埋雷。对比项方法通过标准统计量gdalinfo -stats逐波段比对均值差异 1e-4逐点值抽样 1000 像素求差最大绝对差为 0元信息GetGeoTransform/GetProjection完全一致重构之后我踩过一次重构时把quantize的乘法挪到了外面单波段测试全过整景跑出来偏了 1%原因是测试用例里quantize都是 1.0根本没覆盖到。从那以后我每次动算法层都会强制走一遍带非 1.0 量化系数的合成用例再跑三次对比宁可多花十分钟也别等交付前一天才发现基准全废。希望这些对你有用。本文还有配套的精品资源点击获取