简介面向遥感数据分析与机器学习建模者的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 ReflectanceMTL里同样给出了另一组系数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, encodingutf-8, errorsignore) 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 大气校正为什么表观反射率还不够用表观反射率已经消除了太阳高度角和日地距离的影响但大气还拦在中间。大气中的分子、气溶胶会把一部分太阳辐射散射回传感器路径辐射同时会吸收和散射地表反射信号导致传感器接收到的信号里混入了大量“非地表”贡献。典型表现是表观反射率影像中蓝色波段明显偏亮因为瑞利散射对短波段的贡献最大水体区域本来应该是暗的但看起来发灰。大气校正就是把路径辐射和大气透过率的影响扣除掉还原出真正的地表反射率。目前主流的做法有三类FLAASHENVI内置。精度较高但需要输入成像时间、大气模型热带/中纬度夏季/中纬度冬季等、气溶胶模型、能见度或气溶胶光学厚度。批量处理时每一景都要单独确认参数对脚本不友好。6S模型。学术精度高但输入参数更复杂计算耗时长通常用于单点或少量像元验证。DOSDark 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. 用PythonGDAL跑通批量预处理脚本设计与主循环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 [fB{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, encodingutf-8, errorsignore) 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, fSR_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[fRADIANCE_MULT_BAND_{band}] add coeffs[fRADIANCE_ADD_BAND_{band}] rad dn * mult add # 暗像元统计直方图1%分位 hist, edges np.histogram(rad, bins10000) cum np.cumsum(hist) / max(1, rad.size) idx np.searchsorted(cum, 0.01) lmin (edges[idx] edges[idx1]) / 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, encodingutf-8) as f: w csv.writer(f) w.writerow([scene, status, message]) w.writerows(rows)主循环里每完成一个景就把status改成DONE并回写。这样中途崩了之后只要重新执行脚本会跳过所有DONE状态的任务。这个习惯帮我省下的重跑时间远超写状态逻辑本身花掉的那十几分钟。4. 几何校正、重采样与批量裁剪的参数设计4.1 Landsat8 L1TP产品要做的几何精校正用RPC还是用GCPLandsat8的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, epsgEPSG:32650): RPC正射重采样输出投影和像元尺寸显式指定 warp gdal.Warp( dst, src, optionsgdal.WarpOptions( rpcTrue, dstSRSepsg, xRes30.0, yRes30.0, resampleAlggdalconst.GRA_Cubic ) ) warp None这里gdal.Warp的返回值不是错误码而是输出栅格对象暴露这个对象能直接读取输出尺寸、波段数、投影很方便在批处理循环里做校验。4.2 重采样方法三类选择对结果的影响重采样方法的选择在批量场景下经常被忽略但它直接改变产品的空间纹理。Landsat8的1~7波段本身是30米分辨率但经过几何校正的投影变换后输出栅格需要重新插值。主流选项有几个最近邻nearest、双线性bilinear、三次卷积cubic。方法优点缺点适用场景最近邻保留原始DN/反射率值不会产生新值边缘锯齿明显分类、土地利用制图双线性平滑计算快略微模糊损坏极值植被指数趋势分析三次卷积纹理锐利视觉效果好可能产生过冲负值影像拼接、目视解译做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_TOUCHEDTRUE的组合效果后者会保留任何与边界有交集的像元适合保证覆盖完整性但如果做面积统计必须去掉这个选项否则面积会被系统性高估。按像元窗口裁剪更轻量用gdal_translate的-srcwin参数指定行列起止适合需要把整景切成固定大小瓦片的场景。两种裁剪方式都建议放到几何校正之后因为重投影会改变原始行列位置先裁剪再校正会让边界处出现大量无效插值。5. 批量预处理避坑五个最容易翻车的现场5.1 MTL文件读取乱码导致定标系数全丢现象脚本跑完输出的表观反射率影像全是0或者黑屏。查日志发现read_mtl_coeffs返回的字典是空的。原因MTL文件虽然以UTF-8编码为主但部分Windows环境下用记事本另存后会出现UTF-8 BOMBOM字符会被Python读成“\ufeff”拼在第一个键名前正则匹配自然失灵。解决读取时用encodingutf-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默认值在任务数为奇数时可能造成一两个进程空转设置chunksizeceil(任务数/进程数)能提升尾部负载均衡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(processes4) as pool: results pool.map(run_one, scene_list, chunksizemath.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预处理的路上少走几趟弯路。本文还有配套的精品资源点击获取