简介这套Landsat8影像批量预处理方案面向人工智能与机器学习从业者重点解决遥感数据清洗、云遮挡去除、辐射与大气校正、波段合成、光谱指数计算等特征工程问题。压缩包约46.93MB共含16个文件核心为2个Python批处理脚本配合GeoTIFF示例影像如GMTED2km.tif、README说明文档、JSON/XML配置及.idea工程文件与备份便于对照理解项目工程结构。借助脚本与文档可系统梳理图像加载、云处理、波段校正、特征创建、数据尺度统一与结果保存的完整流程并能以matplotlib/seaborn做结果验证适合在土地覆盖分类、植被监测、灾害检测等场景中进一步建模。目前已有149人学习面向具备Python基础、希望将rasterio、numpy、scikit-image等库用于实际遥感预处理工作的入门与进阶用户。资源来源于网络分享仅供学习交流使用。1. Landsat8影像数据预处理为什么批量预处理才是遥感落地的真正分水岭做遥感的人都有这种经验从USGS拿到一景Landsat8 L1级数据打开一看是16位DN值直方图挤成一团云像棉被一样盖着半幅影像——这样的数据直接丢给机器学习模型出来的结果基本是玄学。Landsat8影像数据预处理从来不是可选项而是决定下游分类、植被监测、变化检测能否成立的前置条件。这个项目给出了一套完整的批量预处理方案核心是PreprocessL8.py脚本覆盖了辐射定标、大气校正、云检测、波段组合与特征输出并且天然支持文件夹级别的批量跑批。适合三类人一是刚接触遥感数据、被L1级产品搞得焦头烂额的学生二是需要给模型喂大量干净样本的算法工程师三是做农业、环境、城市规划项目但不想在预处理上反复返工的从业者。下面我把这套脚本的工程结构、参数逻辑和实际操作一步步拆开讲清楚。2. 读懂工程骨架项目结构、依赖关系与数据组织方式2.1 项目目录结构与数据组织逻辑先把整个压缩包解开目录结构比较清晰。.vscode和.idea分别是VS Code和PyCharm的工程配置目录settings.json里通常写好了Python解释器路径和代码检查规则vsc.xml和misc.xml属于IDE辅助文件跑代码的时候用不上但能帮你省去配置环境的时间。核心文件是PreprocessL8.py它负责预处理主流程base.py提供公共函数和常量GMTED2km.tif是全球地形高程数据分辨率为2公里主要用来做地形校正或辅助大气校正中的高程插值README.md.zbak是被改过后缀的说明文档把.zbak改回.md就能看。数据组织是整个批量处理的地基。我一般建议把Landsat8原始数据按照LC08_景号_日期的格式分文件夹存放每个文件夹里放着该景的全部波段文件B1到B11和MTL元数据文件。脚本默认遍历指定目录下的所有子目录自动识别每个子目录里的Landsat8产品这个设计比手动指定单景文件路径高效得多——你只需要把原始数据按一景一个文件夹丢进去剩下的交给脚本。data/ ├── LC08_20180601/ │ ├── LC08_L1TP_123034_20180601_20180601_01_RT_B1.TIF │ ├── LC08_L1TP_123034_20180601_20180601_01_RT_B2.TIF │ ├── ...B3-B11 │ └── LC08_L1TP_123034_20180601_20180601_01_RT_MTL.txt └── LC08_20180701/ ├── ... └── LC08_L1TP_123034_20180701_20180701_01_RT_MTL.txt2.2 依赖库选型与基础函数解读base.py里封装的依赖和工具函数是整个预处理流程的地基。这个项目主要依赖rasterio读写GeoTIFF、numpy做矩阵运算、pyproj处理投影转换加上concurrent.futures做并行处理。用rasterio而不是GDAL命令行原因在于rasterio的transform和reproject接口可以直接和numpy数组互操作在Python生态里做波段运算时不用频繁落盘效率和代码可读性都好很多。# base.py 核心功能片段简化版 import rasterio import numpy as np from rasterio.warp import calculate_default_transform, reproject from concurrent.futures import ThreadPoolExecutor def load_band(band_path): 加载单个波段返回 (数组, 元数据) with rasterio.open(band_path) as src: return src.read(1), src.profile def stack_bands(band_paths): 把多个波段堆叠成一个多光谱数组 bands [] for path in band_paths: data, _ load_band(path) bands.append(data) return np.stack(bands, axis0)load_band返回数组和profile两个对象数组用于计算profile用于后续写出时保留地理空间信息。stack_bands用np.stack沿第一维堆叠波段生成(波段数, 行, 列)的三维数组这是后续所有波段运算的标准输入格式。用ThreadPoolExecutor而不是ProcessPoolExecutor是因为rasterio底层GDAL在读取文件时IO密集多线程切换已经能榨干磁盘吞吐多进程反而会增加内存复制开销——一景Landsat8影像解压后全波段接近1.5GB进程拷贝会直接打爆内存。3. 批量预处理核心链路辐射定标、大气校正与云检测的参数设置3.1 辐射定标从DN值到反射率的两个必经步骤Landsat8 L1级产品提供的是量化后的DN值范围在0到65535之间。辐射定标分两步走先把DN值乘上MTL文件里的RADIANCE_MULT_BAND_x加RADIANCE_ADD_BAND_x得到辐射亮度值再通过太阳高度角和日地距离换算成表观反射率。PreprocessL8.py里是直接从MTL解析这些参数不需要手工查表。# PreprocessL8.py 辐射定标核心段 import re from datetime import datetime def parse_mtl(mtl_path): 解析MTL文件提取定标参数 params {} with open(mtl_path, r, encodingutf-8) as f: for line in f: m re.match(r\s*(\w)\s*\s*(.), line.strip()) if m: params[m.group(1)] m.group(2).strip().strip() return params def dn_to_reflectance(dn, mult, add, sun_elev, d): DN转表观反射率 radiance dn * mult add reflectance (np.pi * radiance * d ** 2) / (1320 * np.sin(np.deg2rad(sun_elev))) return reflectanceparse_mtl用正则逐行解析MTL把所有GROUP里的键值对拍平成字典好处是不需要关心MTL的具体层级结构取参数时直接用params[RADIANCE_MULT_BAND_4]就能拿到。日地距离d用DN转反射率时乘的那个平方项是当天日地距离和平均日地距离的比值常见做法是用儒略日查表或插值项目里直接用天文学公式计算。这里的1320是Landsat8的太阳等效辐照度ESUN值不同波段这个值不同实际处理时应该从MTL或传感器常数表逐波段取值。3.2 大气校正的三种路线及本方案的选择大气校正是预处理里最影响结果的一步也是坑最多的一步。业界常用三条路线一是用ENVI的FLAASH模块效果最好但需要输入很多大气参数没法批量自动化二是用6S模型自己写接口灵活但代码量大三是用简化方法——基于暗像元DOS做相对校正不需要外部大气参数适合大范围批量的场景。这个项目的定位是批量预处理所以走的是DOS路线。核心假设是影像里存在反射率近似为零的暗像元比如清洁水体、阴影区这些像元在表观反射率上的抬升量近似等于大气程辐射。脚本自动扫描影像直方图的经验百分位把该值作为大气校正的偏移量扣除。# DOS大气校正 def dos_correction(reflectance, percentile1.0): 暗像元法大气校正 # 取直方图指定百分位的反射率作为haze值 haze np.percentile(reflectance[reflectance 0], percentile) corrected (reflectance - haze) / (1 - haze) return np.clip(corrected, 0, 1)percentile参数取1.0是常见做法表示取所有正值像元反射率分布中最暗的1%作为大气影响估计。如果影像里水体面积大这个值会偏小如果全是山地森林可能要调大到2.5。np.clip把结果限制在0到1之间避免出现负反射率或大于1的异常值。这个方法的缺点是不适合定量遥感分析但对机器学习特征输入来说它保留了波段间的相对光谱关系已经够用。3.3 云检测与云掩膜避免污染样本的关键参数Landsat8影像的云检测有现成的FMask算法可以引用但那个算法实现复杂、依赖多。这个项目用的是简化方案——利用Landsat8的卷云波段B9和热红外波段B10来识别厚云和薄云。B9波段对高空卷云敏感B10热红外对云顶低温敏感两者结合能识别大部分云区。def cloud_mask(b9, b10, b9_threshold0.03, b10_threshold270): 生成云掩膜 # 卷云波段阈值分割 cirrus_mask b9 b9_threshold # 热红外亮温阈值 temp_mask b10 b10_threshold # 合并并做膨胀消除边缘效应 cloud np.logical_or(cirrus_mask, temp_mask) return cloudb9_threshold设为 0.03 是经验值如果影像里薄云较多可以降到 0.02但会有误判雪地的风险。b10_threshold单位是开尔文270K约-3℃以下判为云。生成的二值掩膜最终会写入输出文件模型训练时可以直接将掩膜区域的像元权重置为零。这里有个细节云检测要在辐射定标之后、地形校正之前做因为地形阴影也会降低亮温贸然校正会抳掉部分阴影信息。4. 特征工程在影像预处理中的落地波段组合、光谱指数与PCA降维4.1 波段组合的输出配置与典型应用场景Landsat8的OLI传感器有9个波段加上热红外共11个但不是每个波段对下游模型都有用。常规做法是把波段分成三组自然色组合B4-B3-B2用于人工目视检查地表真实感强假彩色组合B5-B4-B3突出植被适合农业监测地质组合B7-B6-B4对矿物和土壤敏感适合环境调查。PreprocessL8.py允许你同时输出多组波段组合而不只是输出全波段堆叠。# 构建输出波段组合 def select_band_combo(stacked, combo_name): 按名称选择波段组合返回RGB三波段数组 combos { natural: [3, 2, 1], # B4-B3-B2索引从0开始 false_color: [4, 3, 2], # B5-B4-B3 geology: [6, 5, 3], # B7-B6-B4 all_bands: list(range(11)) } idx combos[combo_name] return stacked[idx, :, :]索引从0开始所以natural组合取原波段的第4、3、2波段对应存储顺序里下标3、2、1。all_bands输出全部11个波段适合训练多光谱模型时用。我一般建议分类任务用all_bands目视解译输出一张false_color同时存着方便每次跑完快速核对预处理效果。4.2 光谱指数计算NDVI、NDWI、EVI的实现细节光谱指数是遥感特征工程里性价比最高的一类特征。NDVI是归一化植被指数公式是(NIR - Red) / (NIR Red)对Landsat8来说就是(B5 - B4) / (B5 B4)。NDWI用(B3 - B5) / (B3 B5)提取水体。EVI增强型植被指数加入了蓝色波段B2和调整参数能在高植被覆盖区减少饱和。def spectral_indices(stacked): 计算常用光谱指数 # 波段索引: 0B1, 1B2, 2B3, 3B4, 4B5, 5B6, 6B7... red, blue, green, nir, swir1 stacked[3], stacked[1], stacked[2], stacked[4], stacked[5] ndvi (nir - red) / (nir red 1e-10) ndwi (green - nir) / (green nir 1e-10) evi 2.5 * (nir - red) / (nir 6 * red - 7.5 * blue 1) # 将指数矩阵stack成特征通道 return np.stack([ndvi, ndwi, evi], axis0)分母里加的1e-10是防零除的常规操作虽然遥感数据反射率不可能严格为零但浮点计算和掩膜处理后偶尔会出现0值加上这个保底能避免出现nan或inf。np.stack将三个指数堆叠成新增的3个特征通道和原始波段一起输出。EVI公式里的系数6和7.5是行业标准值不同传感器这套系数几乎一样不用改。做植被相关机器学习任务时我建议把EVI做成必选特征它在高密度植被区的抗饱和能力比NDVI强很多。4.3 数据尺度统一标准差归一化与分位数裁剪做完波段堆叠和指数计算后不同波段的数值范围差异很大直接喂给模型会出问题。Landsat8的短波红外波段和热红外波段量纲不同热红外是开尔文温度动辄300K而反射率波段只有0到1。特征工程里一般用Z-score和分位数裁剪两步解决。def normalize_features(feature_array, clip_percentile(2, 98)): 分位数裁剪 标准差归一化 lo, hi np.percentile(feature_array, clip_percentile) clipped np.clip(feature_array, lo, hi) mean clipped.mean(axis(1, 2), keepdimsTrue) std clipped.std(axis(1, 2), keepdimsTrue) normalized (clipped - mean) / (std 1e-8) return normalizedclip_percentile的(2, 98)是经验参数它先把两侧极端离群值裁掉避免个别亮像元拉高整个波段的方差。std 1e-8防止恒定波段比如全零掩膜区域导致除零。这里有个容易翻车的点mean和std必须按每个波段单独计算不能在整个三维数组上统一计算否则数值范围大的波段会主导整个特征空间模型学到的全是热红外波段的分布差异而丢失了光谱细节。5. 批量落地与常见问题排查跑批失败的血泪经验5.1 批量处理脚本的调度与进度管理批量处理的核心约束不是脚本逻辑而是计算资源与异常隔离。一个文件夹里可能有几十景影像其中一景出问题不应该中断整个队列。PreprocessL8.py用concurrent.futures做并行调度每景影像作为一个任务提交任务内捕获异常并记录日志跑完一景再跑下一景。# 批量处理入口 from concurrent.futures import ThreadPoolExecutor, as_completed import logging logging.basicConfig(filenamepreprocess.log, levellogging.INFO) def process_scene(scene_dir, output_root): 处理单景影像返回状态 try: # 1. 解析MTL # 2. 辐射定标 # 3. 大气校正 # 4. 云检测 # 5. 特征计算与输出 logging.info(f{scene_dir} processed.) return f{scene_dir}: OK except Exception as e: logging.error(f{scene_dir} failed: {str(e)}) return f{scene_dir}: FAIL scene_dirs [d for d in glob.glob(data/*) if os.path.isdir(d)] with ThreadPoolExecutor(max_workers4) as executor: futures [executor.submit(process_scene, d, output) for d in scene_dirs] for future in as_completed(futures): print(future.result())max_workers4是CPU核心数和IO负载的折中值如果机器是8核但磁盘是普通机械硬盘开到4就行如果是NVMe固态开到8也没问题。glob.glob(data/*)列出所有子目录作为单景输入所以你的目录里千万不要混入非Landsat8数据的文件夹否则解析MTL那一步会直接报错。日志文件preprocess.log是排查问题的第一入口每失败一景都会记录具体异常的堆栈方便定位是数据问题还是参数问题。5.2 频繁翻车的六个坑及处理方法坑一MTL文件解析失败。现象是脚本一进parse_mtl就报KeyError: RADIANCE_MULT_BAND_4。原因是下载的数据产品类型不是标准L1TP或者MTL版本过旧缺少部分键值。解决方法是先检查MTL文件里是否有RADIANCE_MULT_BAND_4这个键如果没有说明是Collection 1早期格式需要改用REFLECTANCE_MULT_BAND_4字段或者直接打印MTL的前50行确认结构再调整正则。坑二云掩膜把大片植被误判为云。现象是输出影像上植被区域出现大面积空洞。原因是卷云波段B9阈值0.03太低高海拔地区地表反射率低B9值本身偏低但还没到云的程度。解决方法是把b9_threshold调高到0.04或0.05同时增加一个限制条件——只有B9和B10同时超阈值的像元才判为云避免单波段误判。坑三输出GeoTIFF的地理信息丢失。现象是合成的RGB影像用GIS软件打开时没有坐标或者位置偏差了几公里。原因是写出时只用了rasterio.open的默认参数没有把原始波段的transform和crs传进去。解决方法是写文件时显式带上参考波段的profile。def save_geotiff(output_path, data, profile): 带空间信息写出 profile.update(dtypefloat32, countdata.shape[0], compresslzw) with rasterio.open(output_path, w, **profile) as dst: dst.write(data.astype(float32))坑四内存爆掉。现象是处理到第3景时程序被系统杀掉日志里没有报错。原因是多线程同时处理多景影像每景的栈数组和中间变量占满内存。解决方法是检查max_workers是否开得过大同时注意在每景处理完save_geotiff后调用gc.collect()释放大数组或者直接把max_workers1串行处理慢一点但稳定。坑五影像拉伸后颜色发灰。现象是输出的真彩色组合图片看起来对比度极低。原因是反射率数据本身在0到0.3之间分布直接用matplotlib的imshow默认拉伸范围不对。解决方法是显示前先做2%线性拉伸或者把数据缩放到0-255整型后再显示。预处理流程里常见做法是对输出单独生成一张_stretched.tif用于目视检查。坑六投影坐标系不一致。现象是多景影像拼接后边界处地物错位。原因是不同景的UTM分带可能不同或者个别景是WGS84地理坐标。解决方法是处理完单景后统一用rasterio.warp.reproject转到目标CRS常见做法是先挑一景影像的CRS作为基准其他影像全部重投影到该基准上。6. 进阶实践把预处理输出直接接入机器学习训练管线预处理不是终点产出高质量特征才是目标。这套脚本运行后会在输出目录里生成每景对应的_features.tif和_cloudmask.tif前者是波段加指数的多通道特征文件后者是云区域掩膜。把这个输出接入机器学习训练管线关键一步是做好样本采样——遥感影像一个像元就是一个样本但相邻像元高度相关直接全量训练会导致模型过拟合到空间自相关上。我之前做土地覆盖分类时踩过这个坑把整景影像所有像元灌进随机森林训练精度99%验证精度却摔到70%多。原因是训练集和验证集来自同一景影像空间相邻像元的光谱几乎一致模型记住的是位置而不是地物规律。从那以后我每次做遥感训练都强制走一遍空间解耦采样每景影像先按网格切块从不同空间位置抽样本确保训练和验证样本来自不同地理区域这样精度才有参考价值。# 从特征影像中采样训练样本 def sample_pixels(feature_path, mask_path, labels_path, samples_per_class2000): 空间解耦采样 with rasterio.open(feature_path) as src: features src.read() # (C, H, W) with rasterio.open(mask_path) as src: cloud src.read(1) with rasterio.open(labels_path) as src: labels src.read(1) valid (cloud 0) (labels 0) valid_pixels np.argwhere(valid) # 每隔5个像元采一个降低空间自相关 sampled_idx valid_pixels[::5] np.random.shuffle(sampled_idx) X features[:, sampled_idx[:, 0], sampled_idx[:, 1]].T y labels[sampled_idx[:, 0], sampled_idx[:, 1]] # 按类别数量均衡 unique_classes np.unique(y) result_x, result_y [], [] for cls in unique_classes: cls_idx np.where(y cls)[0] if len(cls_idx) samples_per_class: cls_idx np.random.choice(cls_idx, samples_per_class, replaceFalse) result_x.append(X[cls_idx]) result_y.append(y[cls_idx]) return np.vstack(result_x), np.concatenate(result_y)valid_pixels[::5]是关键步把候选像元隔5个采样等效于间距约150米30米分辨率乘5这个间隔大幅削弱了相邻像元的光谱相关性。cloud 0保证样本不会取自云区域。按类别均衡那里用了np.random.choice对样本量不足的类别做无放回抽样超过阈值的类别截断到samples_per_class防止优势类别碾压模型。这套流程我在做植被监测项目时反复跑过Landsat8的预处理脚本处理完约20景影像采样得到约8万个像元样本用随机森林和XGBoost各跑了一遍整体分类精度比不做空间解耦高了近7个百分点其中差异主要在山地和阴影区的类别上体现。预处理到特征工程再到模型训练整套链路是贯通的任何一个环节的参数选择都会传导到最终结果。希望这套批量预处理的拆解和踩坑记录能帮到你少走几趟我当年走过的弯路。本文还有配套的精品资源点击获取
