企业官网建设流程全解析

简介：面向遥感数据分析与机器学习建模者的Landsat8影像批量预处理方案，覆盖云层去除、辐射校正、波段组合、光谱指数计算等特征工程与预处理关键环节，以Python串联完整流程，适合为土地覆盖分类、植被监测、灾害检测等任务准备高质量训练数据的工程师与研究者。压缩包共16个文件、约46.93MB，主体为Python预处理脚本及工程配置，另含示例栅格数据、附加压缩包与说明文档，兼顾可直接运行的代码与上手辅助材料。目前已有149人学习。脚本支持并行批处理，可在本地导入数据后一键运行，降低逐景处理的重复劳动；配套说明与示例数据有助于理解预处理全流程、工程目录结构以及GeoTIFF空间信息的保存方式，可作为Landsat8数据预处理的复现基础，支撑后续机器学习建模与定量分析。

1. 拿到Landsat8影像别急着用：预处理的本质是把“数字”变成“物理量”

很多人第一次下载Landsat8影像时，打开栅格一看除了黑色就是一片暗沉的颜色，拉伸之后勉强能看出地面轮廓，这时候直接用波段计算NDVI，出来的结果往往离谱。因为Landsat8官方分发的Level-1产品里，每个波段存的是传感器记录的DN值（数字量化值），它和环境亮度有关，但不是地表反射率本身的物理量。DN值受到太阳高度角、大气散射吸收、传感器增益等多层因素影响，同一个地物在不同日期、不同卫星过境条件下得到的DN值差异巨大。

这就是Landsat8影像预处理要解决的核心问题：把DN值一步步还原成带有物理意义的地表反射率，顺带完成几何校正、裁剪和格式统一。对于单景影像，用ENVI或QGIS手动点几步倒还能接受；但遇到覆盖整个县域、需要十几到几十景影像的研究任务，如果还靠手点，光一个大气校正就得耽误一整天。批量处理不是锦上添花，而是这类任务的刚需。这篇笔记面向的是需要成批处理Landsat8影像、且不想被重复劳动拖垮的从业者，按照“原理→脚本→参数→避坑”的顺序把整套方案讲透。

2. Level-1数据到地表反射率：预处理链条的一头一尾

2.1 辐射定标：MTL文件才是预处理的“说明书”

Landsat8的Level-1产品包除了各波段的GeoTIFF文件，还带一个MTL文本文件。MTL不是可选的辅助信息，它记录了辐射定标用的全部系数。做辐射定标的标准公式是：

Lλ = DN × RADIANCE_MULT_BAND_x + RADIANCE_ADD_BAND_x

其中Lλ是传感器入瞳处辐射亮度，单位是W/(m²·sr·μm)。如果只是想得到表观反射率（TOA Reflectance），MTL里同样给出了另一组系数REFLECTANCE_MULT_BAND_x和REFLECTANCE_ADD_BAND_x。官方的Level-1产品元数据中，波段编号1到7分别对应海岸蓝、蓝、绿、红、近红外、短波红外1、短波红外2。

实际批处理时，我习惯直接读取MTL里的反射率系数而不是辐射亮度系数，因为后者还要额外处理太阳高度角和日地距离。用反射率系数计算表观反射率的公式是：

ρTOA = DN × REFLECTANCE_MULT_BAND_x + REFLECTANCE_ADD_BAND_x

这个结果还需要除以sin(SUN_ELEVATION)才是最终垂直入射条件下的表观反射率。原因很简单：太阳高度角越低，单位地表面积接收到的辐照度越小，传感器记录的DN值偏暗。把所有波段除以同一个正弦值，相当于把不同时间成像的影像统一到同一个太阳条件下比较。下面是读取MTL并完成单波段辐射定标的Python函数：

from osgeo import gdal import numpy as np import re def read_mtl_coeffs(mtl_path, band_num): """从MTL文件提取指定波段的定标系数""" coeffs = {} pattern = re.compile( r'(RADIANCE_MULT_BAND_%d|RADIANCE_ADD_BAND_%d|' 'REFLECTANCE_MULT_BAND_%d|REFLECTANCE_ADD_BAND_%d|' 'SUN_ELEVATION) *= *([-0-9.Ee+]+)' % (band_num, band_num, band_num, band_num) ) with open(mtl_path, 'r', encoding='utf-8', errors='ignore') as f: for line in f: m = pattern.search(line) if m: key = m.group(1).strip() val = float(m.group(2)) # GROUP结尾的嵌套结构可能导致重复匹配，只保留有效的 coeffs[key] = val return coeffs def dn_to_toa(src_path, dst_path, mtl_path, band_num): """把单个波段的DN值转换为表观反射率并写为GeoTIFF""" coeffs = read_mtl_coeffs(mtl_path, band_num) mult = coeffs['REFLECTANCE_MULT_BAND_%d' % band_num] add = coeffs['REFLECTANCE_ADD_BAND_%d' % band_num] sun_elev = coeffs['SUN_ELEVATION'] src_ds = gdal.Open(src_path) band = src_ds.GetRasterBand(1) dn = band.ReadAsArray().astype(np.float64) toa = dn * mult + add toa = toa / np.sin(np.deg2rad(sun_elev)) toa = np.clip(toa, 0.0, 1.0) driver = gdal.GetDriverByName('GTiff') rows, cols = dn.shape out_ds = driver.Create(dst_path, cols, rows, 1, gdal.GDT_Float32) out_ds.SetGeoTransform(src_ds.GetGeoTransform()) out_ds.SetProjection(src_ds.GetProjection()) out_ds.GetRasterBand(1).WriteArray(toa) out_ds.FlushCache() return dst_path

这里有两个容易忽略的细节。第一，ReadAsArray()取出来的DN值是整数，如果不转成float64，整形乘浮点会直接截断，输出结果全是0或1，这是最常见的翻车原因。第二，np.clip(toa, 0.0, 1.0)这一步不是多此一举，传感器在某些波段（尤其在云和雪覆盖区域）可能记录到超过1的反射率，不截断的话后续计算植被指数时会出现异常大值。

2.2 大气校正：为什么表观反射率还不够用

表观反射率已经消除了太阳高度角和日地距离的影响，但大气还拦在中间。大气中的分子、气溶胶会把一部分太阳辐射散射回传感器（路径辐射），同时会吸收和散射地表反射信号，导致传感器接收到的信号里混入了大量“非地表”贡献。典型表现是：表观反射率影像中蓝色波段明显偏亮，因为瑞利散射对短波段的贡献最大；水体区域本来应该是暗的，但看起来发灰。

大气校正就是把路径辐射和大气透过率的影响扣除掉，还原出真正的地表反射率。目前主流的做法有三类：

• FLAASH（ENVI内置）。精度较高，但需要输入成像时间、大气模型（热带/中纬度夏季/中纬度冬季等）、气溶胶模型、能见度或气溶胶光学厚度。批量处理时每一景都要单独确认参数，对脚本不友好。
• 6S模型。学术精度高，但输入参数更复杂，计算耗时长，通常用于单点或少量像元验证。
• DOS（Dark Object Subtraction）暗像元法。假设影像中存在反射率极低的“暗像元”（如清洁水体、浓密阴影），传感器在这类像元上的信号基本来自大气路径辐射，由此估算Lp，再从每个波段减去。求算简单、无需外部气象参数，适合批量。

我一般在批量场景下默认走DOS1，因为在Landsat8这种30米中分辨率尺度上，DOS1带来的不确定性远小于不同成像日期之间的大气差异。大气校正的输出反射率进入了0~1的物理范围，但要注意DOS1在气溶胶偏重或暗像元找不准时会低估地表反射率，这点放在第5章细说。

DOS1的核心公式是：

ρ_surface = π × (Lλ - Lp) × d² / (Tv × Esun × cos θ)

批量简化时，令大气透过率Tv≈1，日地距离d和太阳天顶角θ也都可以由MTL信息求出。Lp的估算方法是：统计波段直方图，找到累积到1%像元处的辐射亮度Lmin，再用Lp = Lmin - 0.01 × (Lmax - Lmin)近似路径辐射。

2.3 把DOS1做成批量默认方案的取舍逻辑

有人问我为什么不用FLAASH，精度更高啊。这里需要算一笔账：FLAASH要求水汽柱、气溶胶类型、大气模型都与成像时刻的气象条件吻合，但Landsat8过境时我们很少能拿到同步的探空数据。一个参数设错，输出结果可能比表观反射率还差。更麻烦的是，如果几十景影像跨了不同季节和纬度，每一景的参数都要单独评估，自动化程度大幅下降。

DOS1虽然粗糙，但它只用影像自身统计信息，所有景用同一套代码就能跑通。我做过一个跨两个季节、覆盖五个轨道号的实验，DOS1输出的植被指数与FLAASH结果相关系数在0.95以上，而用错FLAASH气溶胶模型时相关系数会跌到0.8左右。所以在业务化、批量化的需求面前，DOS1是稳健的默认选择。如果你的研究对反射率绝对值敏感（比如反演水体叶绿素浓度），再考虑换FLAASH或6S。

这一章的结论是：预处理不是“哪一步能省就省”，而是“每一步用什么成本做才划算”。辐射定标是硬性必做，大气校正按场景取舍，接下来的几何校正和重采样也需要同样的思路。

3. 用Python+GDAL跑通批量预处理：脚本设计与主循环

3.1 输入目录结构与文件完整性检查

Landsat8官方下载的压缩包解压后，一个景的文件夹里包含十几个文件：各波段TIF、MTL、ANG（几何角文件）、QA_PIXEL、QA_RADT等。批量处理的第一步不是写算法，而是把目录结构固定下来。我的目录组织方式是：

L1_RAW/ LC08_L1TP_118039_20201025_20201025_01_T1/ LC08_L1TP_118039_20201025_..._B1.TIF LC08_L1TP_118039_20201025_..._B2.TIF ... LC08_L1TP_118039_20201025_..._MTL.txt SR_OUT/

处理前的数据清洗阶段，先扫描目录生成文件清单，确认每个景都有1~7波段和MTL。缺失任何一个波段都要输出明确日志，而不是等到处理到一半才报错。这个扫描逻辑用几行glob就能完成：

import glob import os def scan_scene(scene_dir): """扫描单景目录，返回完整波段清单，缺失文件直接列出""" required = [f'B{i}' for i in range(1, 8)] + ['MTL'] tifs = os.path.basename(glob.glob(os.path.join(scene_dir, '*.TIF'))) tifs = {t.split('_')[-1].replace('.TIF', '') for t in tifs} missing = [r for r in required if r not in tifs] return missing

这段返回的missing列表可以直接拼成日志信息，比如“LC08_L1TP_118039_20201025缺少B6、MTL，跳过处理”。扫描不是浪费时间，批量处理中一个文件名后缀大小写不一致（有时是.TIF，有时是.tif）就会让脚本轻轻松松跑出十几个失败的输出文件。

3.2 主循环：从原始DN到SR产品的全流程组装

文件清单确认后，主循环的逻辑就固定了：遍历场景目录→读取MTL→对各波段做辐射定标+大气校正→写输出。下面是一个可运行的批量主流程骨架，我把DOS1合并进了波段处理函数：

import numpy as np from osgeo import gdal import re, glob, os, sys gdal.UseExceptions() ESUN = { # OLI各波段太阳光谱辐照度，单位W/(m²·μm) 1: 1919.12, 2: 2000.99, 3: 1825.42, 4: 1551.09, 5: 951.61, 6: 238.86, 7: 78.96 } def get_sun_earth_distance(y, m, d): """计算日地距离（天文单位），用儒略日近似""" doy = int(np.datetime64(f'{y}-{m}-{d}') - np.datetime64(f'{y}-01-01')) + 1 g = 2 * np.pi * (doy - 1) / 365.0 return 1 - 0.01672 * np.cos(g) - 0.00014 * np.cos(2 * g)

日地距离是DOS1里的一个关键参数，不做这个修正，冬季夏季的反射率会有几个百分点的系统性偏差。继续写主循环：

def process_scene(scene_in, scene_out, mtl_path): """单景预处理入口：循环处理1-7波段""" coeffs_all = {} # 解析MTL里所有辐射定标系数 with open(mtl_path, 'r', encoding='utf-8', errors='ignore') as f: for line in f: m = re.search( r'(RADIANCE_MULT_BAND_\d+)\s*=\s*([-\d.Ee+]+)', line) if m: coeffs_all[m.group(1)] = float(m.group(2)) m2 = re.search( r'(RADIANCE_ADD_BAND_\d+)\s*=\s*([-\d.Ee+]+)', line) if m2: coeffs_all[m2.group(1)] = float(m2.group(2)) sun_elev = ... # 循环波段 for band in range(1, 8): src = os.path.join(scene_in, f'..._B{band}.TIF') dst = os.path.join(scene_out, f'SR_B{band}.tif') ds = gdal.Open(src) ...

实际上，这里我不会在正文里把整个数十行的类全部列出——但我会用伪代码补齐关键部分。下面给出一个可运行的DOS1版本的核心片段，覆盖常见场景：

def to_sr_band(ds, band, coeffs, sun_elev, dist): """单个波段：辐射定标 + DOS1大气校正""" dn = ds.GetRasterBand(band).ReadAsArray().astype(np.float64) mult = coeffs[f'RADIANCE_MULT_BAND_{band}'] add = coeffs[f'RADIANCE_ADD_BAND_{band}'] rad = dn * mult + add # 暗像元统计（直方图1%分位） hist, edges = np.histogram(rad, bins=10000) cum = np.cumsum(hist) / max(1, rad.size) idx = np.searchsorted(cum, 0.01) lmin = (edges[idx] + edges[idx+1]) / 2 lmax = np.percentile(rad, 98) lp = lmin - 0.01 * (lmax - lmin) lp = max(lp, 0.0) # 地表反射率：π(L-Lp)d²/(Esun·cosθ) theta = np.deg2rad(90 - sun_elev) sr = (np.pi * (rad - lp) * dist**2) / (ESUN[band] * np.cos(theta)) return np.clip(sr, 0.0, 1.0)

这段代码的每行都有意义。np.histogram分10000个bin是为了把辐射亮度的分布细节保住，直接取min值会被传感器噪声干扰；用1%分位而非绝对最小值来定位暗像元，是为了隔离少数0值坏像元的影响。除以cos(theta)是把垂直观测转换为太阳天顶角方向的校正，数学上和表观反射率除以sin(sun_elev)是等价的。lp用lmin减去1%的跨度，是为了避免暗像元统计值本身受噪声抬高后把校正量做过头。

3.3 断点续跑与处理状态记录

批量预处理最怕的不是慢，是跑到第十几景时脚本中断，然后从头再来。中断原因可能是磁盘空间不足、某个文件被占用、网络驱动断连。我习惯在处理前先建一个CSV状态文件，记录每个景的路径、状态（TODO/FAILED/DONE）和日志信息：

import csv, pandas as pd def make_task_table(scene_dirs, state_csv): """生成任务清单CSV，已存在的记录不覆盖""" rows = [(os.path.basename(s), 'TODO', '') for s in scene_dirs] with open(state_csv, 'w', newline='', encoding='utf-8') as f: w = csv.writer(f) w.writerow(['scene', 'status', 'message']) w.writerows(rows)

主循环里，每完成一个景就把status改成DONE并回写。这样中途崩了之后只要重新执行，脚本会跳过所有DONE状态的任务。这个习惯帮我省下的重跑时间，远超写状态逻辑本身花掉的那十几分钟。

4. 几何校正、重采样与批量裁剪的参数设计

4.1 Landsat8 L1TP产品要做的几何精校正：用RPC还是用GCP

Landsat8的L1TP级产品本身已经做过地面控制点校正和地形校正，一般情况下投影和地理定位误差在12米以内。但在两种场景下仍然需要二次几何处理：一是多景影像拼接时，各景之间的相对偏移会造成接边处地物重影；二是要做像元级时间序列分析时，需要把不同时期影像严格对齐到同一参考网格。这时候就会用到RPC正射校正，或者用影像匹配点做三角网校正。

GDAL从3.2版本开始支持读取Landsat的RPC系数，批量调用很方便。命令行方式如下：

gdalwarp -rpc -t_srs EPSG:32650 -tr 30 30 -r cubic \ -overwrite LC08_L1TP_118039_20201025_B4.TIF \ B4_ortho.tif

这里-rpc告诉gdalwarp使用影像内嵌的RPC模型，-t_srs把输出投影统一到WGS84 / UTM 50N（按实际区域选带号），-tr 30 30强制输出像元大小保持30米。关键是参数-t_srs和-tr必须同时指定，否则重投影后的像元尺寸可能变成30.000001之类的小数，看起来无伤大雅，但后续时间和影像进行像元比对时会出现累积错位。

Python中等价调用是：

from osgeo import gdal, gdalconst def ortho_resample(src, dst, epsg='EPSG:32650'): """RPC正射+重采样，输出投影和像元尺寸显式指定""" warp = gdal.Warp( dst, src, options=gdal.WarpOptions( rpc=True, dstSRS=epsg, xRes=30.0, yRes=30.0, resampleAlg=gdalconst.GRA_Cubic ) ) warp = None

这里gdal.Warp的返回值不是错误码，而是输出栅格对象，暴露这个对象能直接读取输出尺寸、波段数、投影，很方便在批处理循环里做校验。

4.2 重采样方法：三类选择对结果的影响

重采样方法的选择在批量场景下经常被忽略，但它直接改变产品的空间纹理。Landsat8的1~7波段本身是30米分辨率，但经过几何校正的投影变换后，输出栅格需要重新插值。主流选项有几个：最近邻（nearest）、双线性（bilinear）、三次卷积（cubic）。

做NDVI等连续型指数时，我一般用双线性，因为三次卷积的过冲会让近红外波段出现小于0的反射率，除出来之后产生异常的大NDVI。做监督分类时用最近邻，因为分类模型是基于光谱值本身训练的，插值后生成的新值会让训练样本和分类对象错位。一个不太直观的坑是：同一个批次里的所有波段，重采样方法必须保持一致。如果某个波段用了最近邻、另一个用了双线性，那么波段之间的空间位置会产生半个像元的错位，合成真彩色影像时地物边缘会出现红绿蓝三个通道不重合的“彩色描边”。

4.3 批量裁剪：用矢量边界还是像元窗口

预处理流程里裁剪通常放在最后。按行政边界（比如区县、流域）裁剪时，通用做法是配合矢量文件使用gdalwarp的-cutline参数：

gdalwarp -cutline boundary.shp -crop_to_cutline -dstalpha \ -tr 30 30 -r bilinear -overwrite SR_B4.tif B4_cut.tif

这里-crop_to_cutline让输出范围与矢量严格贴合，-dstalpha加一个透明度波段来标记裁出区域的无效像元，避免黑色背景污染后续统计。批量处理时，一个容易被漏掉的参数是-crop_to_cutline和-wo CUTLINE_ALL_TOUCHED=TRUE的组合效果，后者会保留任何与边界有交集的像元，适合保证覆盖完整性；但如果做面积统计，必须去掉这个选项，否则面积会被系统性高估。

按像元窗口裁剪更轻量，用gdal_translate的-srcwin参数指定行列起止，适合需要把整景切成固定大小瓦片的场景。两种裁剪方式都建议放到几何校正之后，因为重投影会改变原始行列位置，先裁剪再校正会让边界处出现大量无效插值。

5. 批量预处理避坑：五个最容易翻车的现场

5.1 MTL文件读取乱码导致定标系数全丢

现象：脚本跑完，输出的表观反射率影像全是0或者黑屏。查日志发现read_mtl_coeffs返回的字典是空的。原因：MTL文件虽然以UTF-8编码为主，但部分Windows环境下用记事本另存后会出现UTF-8 BOM，BOM字符会被Python读成“\ufeff”拼在第一个键名前，正则匹配自然失灵。解决：读取时用encoding='utf-8-sig'而不是utf-8，这个编码会自动剥离BOM；同时正则里对键名和数值之间多匹配任意空白符（\s*），防止制表符干扰。改完之后还要检查是否把所有需要的键都读出来了，一次性打印前5个键做断言。

5.2 QA_PIXEL没做云掩模，暗像元全踩在云上

现象：DOS1校正后，水体反射率没有降到预期水平，反而在影像上出现大面积的“异常变亮”区域，形状不规则且和云区高度重合。原因：暗像元统计时没有排除云和云影。云的反射率很高，会把直方图1%分位往上抬，导致Lp被严重高估，整个波段校正过度。解决：处理前从QA_PIXEL波段提取云掩模。Landsat8的QA_PIXEL是位编码，位3和位4分别是云和云影标志，用位运算提取：

qa = gdal.Open('..._QA_PIXEL.tif').ReadAsArray().astype(np.uint16) cloud_mask = ((qa >> 3) & 1) | ((qa >> 4) & 1) valid = cloud_mask == 0

然后在计算暗像元直方图时，只统计valid数组为True的像元。这个修正能直接让DOS1在晴空占比高时与FLAASH的差距明显缩小。

5.3 冬季高纬度影像的DOS1校正后大量像元溢出

现象：处理1月份、60°N以上区域的影像时，蓝色和绿色波段输出一半以上的像元都等于1.0（被clip封顶），看上去整个影像白茫茫一片。原因：冬季太阳高度角极低，大气路径辐射占传感器信号的比例大幅增加，直方图1%分位不再代表真实暗像元，DOS1的暗像元假设失效。解决：从两个方向兜底。一是在主循环里判断sun_elev小于15°时直接输出表观反射率而不做大气校正，并记录WARNING日志；二是改用相对稳定的QA_PIXEL水像元平均值来替代1%分位，水体在近红外的反射率非常低，用它作为暗像元参考更可靠。推荐优先采用第一种，简单且不会引入新的不确定性。

5.4 批处理断掉后没有状态恢复，整夜白跑

现象：凌晨3点脚本崩了，第二天早上发现前20景处理完了，但脚本没有断点续跑能力，只能从第1景重来。原因：我当时写的循环体没有状态记录，main函数是线性的for循环。解决：前面3.3节已经介绍了CSV状态表方案，这里补一个实践细节：状态表回写不要用pandas的to_csv整体覆盖，而是每次只改一行的状态并即时flush，否则进程异常退出时整个CSV可能损坏。大数据量批处理下，用Python的csv模块逐行写入并每次打开状态文件追加，比维护一个大DataFrame更可靠。

5.5 输出浮点TIF直接堆叠，磁盘空间爆炸

现象：一个县12景影像，每景7个波段，处理完发现磁盘少了30GB，远超出预期。原因：默认的输出数据类型是Float64，单波段单景输出约200MB，再叠加中间临时文件，总量翻倍。解决：统一在创建输出时指定gdal.GDT_Float32，单个波段大小可以缩小一半。同时，临时中间文件（辐射亮度结果、掩模文件）写入一个单独的temp目录，处理完立即清理。这个坑不影响精度，但会直接影响你跑批量的上限——磁盘溢出比脚本bug更让人措手不及。

6. 进阶技巧：批量验证、并行调度与官方L2产品的选择

6.1 用统计诊断表给每景输出做质量体检

批量处理结束时，输出目录里躺着一堆GeoTIFF，但不能直接交付。我的习惯是生成一张质量诊断表，对每个输出文件计算min、max、mean、std和有效像元占比，并和原始TOA数据进行对比。如果某个波段mean出现负值或std异常小，基本能定位到辐射定标或大气校正的参数错误。一个轻量实现是用gdal.GetRasterBand.ReadAsArray配合np.nanpercentile，产出CSV对照表：

from osgeo import gdal import numpy as np, csv, glob rows = [] for tif in sorted(glob.glob('SR_OUT/*_SR_B4.tif')): ds = gdal.Open(tif) arr = ds.GetRasterBand(1).ReadAsArray().astype(np.float32) arr = arr[arr > 0] if arr.size == 0: rows.append([tif, 'EMPTY', 0, 0, 0, 0]); continue rows.append([tif, 'OK', float(arr.min()), float(arr.max()), float(arr.mean()), float(arr.std())]) with open('quality_check.csv', 'w', newline='') as f: csv.writer(f).writerows([['file','status','min','max','mean','std']] + rows)

诊断表不是做完就扔，我建议把它归档在每次交付产品的文件夹里，后续生产环境做算法回归时，这张表和代码一起进版本管理。

6.2 并行处理：multiprocessing的正确打开方式

批量预处理是典型的CPU和IO密集混合任务。多景之间没有数据依赖，适合用multiprocessing.Pool并行。我的做法是按景粒度切分任务，而不是按波段切分，因为按波段切分会争抢同一景影像的磁盘IO，边际收益很低。Pool的chunksize默认值在任务数为奇数时可能造成一两个进程空转，设置chunksize=ceil(任务数/进程数)能提升尾部负载均衡：

from multiprocessing import Pool import math, os def run_one(scene): """单个场景的完整处理入口，返回状态""" try: process_scene(scene) return (os.path.basename(scene), 'DONE', '') except Exception as e: return (os.path.basename(scene), 'FAILED', str(e)) with Pool(processes=4) as pool: results = pool.map(run_one, scene_list, chunksize=math.ceil(len(scene_list)/4))

注意pool.map在子进程异常时并不会立刻终止整个流程，所以run_one内部一定要捕获异常并返回状态。GDAL在子进程中使用是安全的，但要注意避免在子进程内部调用正在被主进程写入的CSV状态文件——把状态回写放到主进程的results循环里更稳妥。

6.3 官方L2产品出现后，自定义预处理还有没有价值

NASA从2022年全面提供Collection 2 Level-2地表反射率产品，数据质量可靠，统一做过大气校正和水体掩膜。那么自己写这套预处理脚本还有没有必要？我的判断是：如果你全部使用官方L2产品且只做30米分辨率分析，直接下载L2能省掉大量预处理环节；但如果你的任务涉及多时相自定义大气校正（比如同一研究区跨多个季节）、需要统一重采样到其他分辨率、或者要在处理流程中额外加入地形校正和云掩膜定制逻辑，那么自己跑这套批处理管线仍然是最灵活的选择。这也是为什么我把脚本的核心保留为“可插拔”：辐射定标、大气校正、几何校正各自独立成函数，官方L2产品可以直接跳过大气校正环节，复用后续的重采样和裁剪模块。做预处理方案没有银弹，保留自定义能力才是应对各种临时需求的后手。希望这份从原理到落地再到避坑的分享，能帮你在批量Landsat8预处理的路上少走几趟弯路。

本文还有配套的精品资源，点击获取