企业官网建设流程全解析

开头

做遥感数据处理的老哥们，估计都有过这样的经历：拿到一景影像，第一件事不是急着做分类或者反演，而是先做预处理。预处理里面最绕不开的两个活儿，一个是按研究区边界把影像“切”出来，另一个就是把大气给影像带来的“滤镜效果”去掉——前者叫裁剪，后者叫大气校正。这两步做好了，后续无论是算植被指数、做地物分类，还是提取地表温度，数据才真正能用。

我之前在几个遥感项目里，被裁剪和大气校正折磨过不少次，走了不少弯路。比如最初用的软件切影像，界面操作慢不说，到了批处理的时候——几十景影像挨个点，点得人怀疑人生。后来发现用Python搭配GDAL来搞，脚本一跑，散散步的功夫全出来了。这篇文章就把我实际工作中用Py和GDAL做遥感影像裁剪和大气校正的整套思路、踩坑记录和可直接抄的代码，复盘一遍。这里给了具体的方案、参数怎么调、为什么这么调，尽量让第一次接触这个的同学也能跟着复现出来。

这篇文章适合哪些人看？就是用Python做遥感数据处理的学生、科研人员和工程师，尤其是到了预处理环节，想通过脚本取代手工操作、把流程自动化、结果可复现的朋友。文章不涉及高深理论推导，重点在于把SHARE出来的路子讲清楚。

1. 整体设计与技术选型：为什么是Python + GDAL

1.1 预处理流程中的两个关键瓶颈

先捋一下遥感影像预处理的常规流程：数据读取 → 辐射定标 → 大气校正 → 几何校正 → 裁剪。其中大气校正和裁剪，是流程中最容易出问题、也最耗时的环节。

裁剪这块，问题不在于“不会切”，而在于“怎么把不规则的矢量边界准确套到影像上”。很多新手直接拿矩形范围去框，结果研究区是个多边形，切出来一大堆无关区域，数据量白白翻了几倍。而且，如果影像很大（比如一景Sentinel-2的10米波段，全幅下来好几个GB），内存不够的话，加载过程卡到人崩溃。

大气校正就更不用说了。大气对地表反射率的影响是叠加在信号里的，不做校正，你算出来的NDVI、地表反照率这些参数，数值整体偏移，后续反演精度直接崩。很多刚入行的同学总会问一个问题：我用ENVI做一次大气校正也就几分钟，为什么要用代码重写一遍？答案很简单——当数据量变成几十景上百景的时候，手动做一遍要几天几夜，而脚本只要跑一次。

1.2 GDAL为什么是首选

在地理空间数据处理圈里，GDAL（Geospatial Data Abstraction Library）的基本地位不用多说。它支持几十种栅格格式，又是命令行和Python绑定的，裁剪、重投影、格式转换、像素级计算都能搞定。我在实践中选择GDAL的原因，总结下来有三条：

• 跨格式能力极强：无论是GeoTIFF、ENVI格式，还是NetCDF，GDAL都读得懂。这意味着前面数据从哪个平台下的，不重要，反正都能转到统一格式来处理。
• 自带命令工具：gdalwarp、gdal_translate这俩命令行工具，配合Python的subprocess调用，或者直接在Python中绑定使用，既能处理大文件，内存控制又很稳。
• 生态成熟，坑少有人填：GDAL二十多年的迭代，边界情况处理得很完善。比如NoData值怎么处理，坐标参考系统怎么转换，它都内置了非常成熟的方案。

1.3 大气校正路径的取舍

大气校正的实现路径有多条：6S模型、MODTRAN模型、LOWTRAN，以及基于影像本身的暗像元法（DOS）。其中6S模型在科研里用得非常多，因为它物理基础相对扎实，输入参数相对容易获取。MODTRAN精度更高，但是参数多且商业授权昂贵，不适合普通用户。

在Python环境里，Py6S库封装了6S模型的调用接口，不用去折腾Fortran编译版本的6S，既能做单波段的大气校正，也能对多光谱影像做逐波段模拟。除了Py6S这类辐射传输模型方法之外，我在处理一些不需要极高精度的场景时，也会用DOS（暗像元减法）方法做一个快速大气校正，虽然精度和物理模型比差了，但胜在简单快速、不需要大气参数数据，作为初筛完全够用。

这一章节总结下来：裁剪用GDAL的gdalwarp，大气校正优先用Py6S + 6S模型，快速预处理场景就用DOS作为替代。下面具体拆解每一步。

2. 数据准备与核心参数解析

2.1 数据格式与组织方式

我习惯在项目文件夹下把数据组织成这样的结构：

project/ ├── data/ │ ├── raw/ # 原始影像 │ │ ├── scene1.tif │ │ ├── scene2.tif │ │ └── ... │ ├── vector/ # 研究区边界 │ │ ├── roi.shp │ │ └── roi.dbf │ └── output/ │ ├── clipped/ # 裁剪结果 │ └── corrected/ # 大气校正结果 ├── scripts/ │ ├── clip.py │ ├── atmospheric_correction.py │ └── utils.py └── config.yaml

为什么要这么组织？因为遥感项目一旦数据量上来，手动翻找文件会疯。数据、脚本、中间结果、最终结果全部隔离，处理流程的每一步都有对应输出，出问题了也好排查是哪个环节挂的。

影像格式我统一建议用GeoTIFF。原因有几个：支持压缩，自带地理参考信息，GDAL读写完全没问题，而且是大多数平台（GEE、Sentinel-2官方数据分发的L2A产品）的标准格式之一。

2.2 坐标参考系统与范围信息核对

很多裁剪失败的案例，最后查来查去，问题出在“坐标系不一致”。边界矢量文件用的是WGS84经纬度（EPSG:4326），影像却是UTM投影坐标系（比如EPSG:32650），直接拿来算交集，出来的结果就是空白或者报错。

所以，在裁剪之前，第一件事是检查矢量边界和影像的坐标系是否一致。代码如下：

from osgeo import gdal, ogr # 打开影像获取投影 ds = gdal.Open('data/raw/scene1.tif') proj = ds.GetProjection() print("影像投影信息:", proj) # 打开矢量边界获取投影 vec_ds = ogr.Open('data/vector/roi.shp') layer = vec_ds.GetLayer() vec_proj = layer.GetSpatialRef() print("矢量边界投影信息:", vec_proj)

如果两边EPSG代码不一致，就用gdal.Warp的dstSRS参数统一到影像的坐标系。这一步很容易被忽略，但却是最关键的坑。EPSG代码可以在读取出投影信息后，用osr模块解析出来，也可以直接在输出信息里看。

2.3 波段范围与位深检查

另外一个容易翻车的地方是影像的位深和波段数。例如Landsat 8 OLI影像有9个波段（包含全色波段），Sentinel-2的L1C产品有12个波段（此处按常规产品组织描述）。做大气校正的时候，一般只处理多光谱波段，全色波段要单独拿出来，不能混在一起。还有像16位整型的影像，直方图展示和后续数值拉伸方式与8位影像完全不同。这些差异如果不处理好，算出来的反射率会差得一塌糊涂。

所以在开工之前，务必跑一下这个命令式检查：

band_count = ds.RasterCount x_size = ds.RasterXSize y_size = ds.RasterYSize print(f"波段数: {band_count}, 尺寸: {x_size} x {y_size}") data_type = ds.GetRasterBand(1).DataType print("数据类型编号:", data_type)

GDAL数据类型编号里，1代表8位无符号整型，2代表16位有符号整型，3代表16位无符号整型，6代表32位浮点型等。这一步确认好了，后面做数值定标、转换浮点型反射率才有底。

3. 基于GDAL的遥感影像矢量裁剪实操

3.1 用gdal.Warp实现精确裁剪

裁剪的核心逻辑其实一句话就能说清楚：让GDAL把影像重采样到矢量边界范围内。而gdal.Warp函数做得非常漂亮——它底层实现了边界对齐、投影转换、重采样一条龙。

基础调用方式是这样的：

from osgeo import gdal gdal.Warp( 'data/output/clipped/scene1_clip.tif', 'data/raw/scene1.tif', cutlineDSName='data/vector/roi.shp', cropToCutline=True, dstNodata=0.0, options=['-overwrite'] )

这四个参数是核心，逐个说清楚：

• cutlineDSName：指定矢量边界文件。这个边界文件可以是Shapefile，也可以是GeoJSON。
• cropToCutline=True：表示输出栅格的范围完全贴合边界。如果不设置这个，输出范围是边界的包围盒（矩形），不是一个多边形，就没达到“裁剪”的效果。
• dstNodata=0.0：设置无效像素值。对于16位影像，NoData通常用0或NaN，具体看原始影像里怎么定义的。
• options=['-overwrite']：允许覆盖已有的输出文件，否则文件存在时会报错。

3.2 重采样方法选择

裁剪必然涉及重采样（因为裁出来的范围大概率不是像素对齐的），重采样方法的选择直接影响裁剪结果质量和后续计算的精度。GDAL支持最近邻、双线性、三次卷积等重采样方法。

我的经验是：

• 做分类、目视解译用最近邻：保持原始像元值，不引入插值造成的光谱信息污染，地物类别边界清晰。
• 做光谱指数反演（NDVI、LAI等）用双线性：平滑效应小，在辐射量上不会太失真。
• 做高精度表面反射率产品、对空间纹理要求高时用三次卷积：更平滑但可能过冲，出现轻微的光谱值越界，应用时要注意截断。

在gdal.Warp中，通过resampleAlg参数控制：

gdal.Warp( 'data/output/clipped/scene1_clip_bilinear.tif', 'data/raw/scene1.tif', cutlineDSName='data/vector/roi.shp', cropToCutline=True, resampleAlg=gdal.GRA_Bilinear, dstNodata=0.0, options=['-overwrite'] )

3.3 批处理加速：一个循环切一百景

在实际生产化流程里，一次性处理几十上百景影像很常见。我的做法是先把所有影像路径放到一个列表里，然后依次调用裁剪函数。这里有一个性能优化点：将矢量边界文件路径和输出目录传参复用，避免重复打开矢量文件。Python的循环调用gdal.Warp，其实效率还不错，因为实际的重采样计算发生在GDAL底层C++代码中，Python这一层只是一个外壳。

一个更具可操作性的批处理模板：

import os from osgeo import gdal raw_dir = 'data/raw' clip_dir = 'data/output/clipped' roi_file = 'data/vector/roi.shp' if not os.path.exists(clip_dir): os.makedirs(clip_dir) for scene in os.listdir(raw_dir): if scene.endswith('.tif'): in_path = os.path.join(raw_dir, scene) out_path = os.path.join(clip_dir, scene.replace('.tif', '_clip.tif')) gdal.Warp( out_path, in_path, cutlineDSName=roi_file, cropToCutline=True, resampleAlg=gdal.GRA_Bilinear, dstNodata=0.0, options=['-overwrite'] ) print(f"完成: {scene}")

有几个细节值得再啰嗦一遍：一是确保目录路径中不要有中文，有些环境下GDAL在处理中文字符路径时会出幺蛾子；二是gdal.Warp本身不会告诉你“裁剪面积比例”，如果想检查裁剪是否过度或不足，可以对比输出文件的大小和像素尺寸是否和预期吻合。

3.4 裁剪结果的边界锯齿与研究区外像素处理

矢量边界如果是高精度多边形，裁剪出的影像边缘往往是斜线切割的，不可避免地会出现锯齿状边缘。这是正常现象，也可以用掩膜来平滑，但高频操作会引入误差。另一个容易被忽略的问题是：研究区外有时会残留黑色填充值或者0值，这个像素在做统计计算（比如计算均值）时会拉低数值。所以在后续计算之前，最好做一个NoData掩膜，把这些背景像素彻底排除。

对于Landsat和Sentinel-2这类影像，原始产品本来就自带QA波段，裁剪之后需要把QA波段也同步裁出来，后续去云、去阴影处理还会用到。我在批处理脚本中，会把所有波段一次性裁出来，比如Landsat 8就把所有波段一起处理，不会单独裁某几个波段再把波段拼回来——因为多次重采样会累积插值误差。

4. 大气校正原理与Py6S实操

4.1 为什么不能跳过大气校正

先举一个非常直观的例子：一景影像中，同一块地物在不同日期的表观反射率相差非常大，但地表真实反射率其实是基本稳定的。差值来源于什么？来源于大气中的分子散射（瑞利散射）、气溶胶散射、水汽吸收等因素，它们会使得传感器接收到的信号既包含了地表反射的贡献，又叠加了一份“大气路径辐射”，还把一部分地表信号衰减掉。

如果不做大气校正，那么做时间序列分析时，变化信息被大气噪声污染，结论完全不可信。这也是很多做长时序生态监测的团队最头疼的地方。大气校正的本质，就是通过模型模拟或统计估计，把这些大气效应从表观反射率里剥离出来，得到地表真实反射率。

4.2 暗像元法（DOS）快速校正

先给出一个简单但能应急的DO方法。暗像元法的假设是：影像上存在很多非常暗的像元（比如清洁水体、深阴影），这些像元的表观反射率几乎全部来自大气路径辐射。通过统计直方图找到1%分位数对应的数值作为暗像元值，然后把它当作大气贡献减掉。

实现起来其实挺快的：

import numpy as np from osgeo import gdal def apply_dos(in_path, out_path, band_nums, dark_percent=0.01): ds = gdal.Open(in_path) # 读取每个波段并统计暗像元 bands_data = [] for b in band_nums: band = ds.GetRasterBand(b) data = band.ReadAsArray() # 只统计非NoData区域 valid_data = data[data > 0] if len(valid_data) == 0: dark_val = 0 else: dark_val = np.percentile(valid_data, dark_percent * 100) bands_data.append((data.astype(np.float32) - dark_val)) print(f"波段 {b} 的暗像元值: {dark_val}") # 写入输出 driver = gdal.GetDriverByName('GTiff') out_ds = driver.Create(out_path, ds.RasterXSize, ds.RasterYSize, len(band_nums), gdal.GDT_Float32) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) for i, arr in enumerate(bands_data): out_ds.GetRasterBand(i + 1).WriteArray(arr) out_ds.FlushCache() print("DOS校正完成")

这个方法最大的局限在于：它把所有波段都用同一个减法，但大气影响在可见光波段强、在近红外波段弱，同一套暗像元减法并不物理自洽。因此DOS只适合快速预览和精度要求不高的场景，真正做科研出结果，还是得用辐射传输模型。

4.3 Py6S + 6S模型的完整流程

在物理模型路线中，6S（Second Simulation of a Satellite Signal in the Solar Spectrum）是被引用最多的模型之一。Py6S是6S模型的Python封装，可以直接在Python里设置参数并获取模拟结果。

我来展示一下完整的单景影像大气校正流程。

首先安装Py6S：

pip install py6s

注意，Py6S内部需要调用Fortran编译的6S可执行文件，安装时一般会自动下载编译，但部分Linux环境下可能需要自行编译。Windows环境直接用pip安装基本没问题，但如果在Linux服务器上装，可能要先把6S编译好。

接下来是调用大气参数模拟：

from py6s import SixS, AtmosProfile, AeroProfile, Geometry, Wavelength def atmospherically_correct(in_path, out_path, bands_meta, geom_info): """ 输入：裁剪后的影像路径、输出路径、波段波长信息、几何信息 """ s = SixS() # -------- 关键参数设置 -------- s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer) s.aero_profile = AeroProfile.PredefinedType(AeroProfile.Continental) s.aot550 = 0.2 # 550nm处气溶胶光学厚度，默认0.2 # 几何信息：太阳天顶角、太阳方位角、卫星天顶角、卫星方位角 s.geometry = Geometry.User() s.geometry.solar_z = geom_info['solar_zenith'] s.geometry.solar_a = geom_info['solar_azimuth'] s.geometry.view_z = geom_info['view_zenith'] s.geometry.view_a = geom_info['view_azimuth'] # 逐波段进行6S模拟，获取大气校正系数 for meta in bands_meta: s.wavelength = Wavelength(meta['wavelength_nm']) s.run() # 6S输出中，关键参数有：x_a、x_b、x_c（用于表面反射率计算） x_a = s.outputs.coef_xa x_b = s.outputs.coef_xb x_c = s.outputs.coef_xc print(f"波段 {meta['name']} 系数: xa={x_a:.6f}, xb={x_b:.6f}, xc={x_c:.6f}")

这里要详细解释三个系数：

• x_a：乘性系数，主要用于修正大气透过率和散射的影响。
• x_b：加性系数，对应大气路径辐射项。
• x_c：常量偏移项。

表面反射率的计算公式是：

[ \rho = \frac{y}{1 + x_c \cdot y} ]

其中：

[ y = x_a \cdot L_{toa} - x_b ]

这里L_toa是传感器接收到的辐亮度值。所以大气校正的过程，本质是对每个波段找到一组系数，然后用上述公式对每一个像元做逐像元的运算。

在拿到逐波段的6S模拟系数后，就可以编写批量处理的代码了。因为6S模拟在波段很多的时候，每跑一次都要几秒钟，如果影像有10个波段，模拟就得几十秒，但每个波段的系数只取决于波段波长和几何条件，跟像元数无关，所以一次性算好系数，后面逐波段校正就只是一个简单的NumPy数组计算了。这样性能上是完全可接受的。

4.4 从DN值到表观反射率再到地表反射率

在跑Py6S之前，还有一个很容易忽略的点：6S模型的输入通常是辐亮度，不是原始DN（数字量化值）。也就是说，你从影像里读出来的DN值，需要先做辐射定标——乘以增益加上偏移，变成辐亮度。

例如某个传感器头文件里已经给出了各波段的增益（radiance_mult）、偏移（radiance_add），那么辐亮度就是：

l_toa = dn * radiance_mult + radiance_add

如果你的影像产品本身就是表观反射率产品（比如部分分发平台直接提供了SR），那可以跳过这步，直接当作6S的输入。但很多L1级产品还是DN值，这一步必须做。

完整的事后校正流程在代码中看起来大致是这样：

import numpy as np from osgeo import gdal def correct_band(data_dn, mult, add, xa, xb, xc): # 辐射定标：DN -> 辐亮度 l_toa = data_dn.astype(np.float32) * mult + add # 6S系数反演 y = xa * l_toa - xb rho = y / (1 + xc * y) # 负值和无效值处理 rho[rho < 0] = 0 rho[rho > 1] = 1 return rho

校正完之后，得到的输出是0到1之间的地表反射率，如果要存储成8位影像方便目视，就乘以10000再转成整型存盘；如果要做定量分析，建议还是直接用Float32存储，保留足够的数值精度。

注意：Py6S模拟的几何参数，太阳天顶角等可以从影像元数据或者头文件里读取。很多L1产品自带的MTL文本里都有SOLAR_ZENITH_ANGLE这类字段（各传感器命名略有差异）。这些参数是逐景变化的，不能复制粘贴，必须实际读取。

4.5 大气校正的常见误区

我见过有人把大气校正数据当成“直接可用的绝对反射率”，这其实是一个普遍误区。6S模型模拟的精度，极大地依赖于你输入的大气参数（气溶胶类型、水汽含量、气溶胶光学厚度等）。如果你设置的参数和实际情况偏差太大，校正结果反而不如直接用表观反射率做相对分析。所以我的习惯是：如果只是做某个单景的相对分类，大气校正优先级没那么高；但如果要做跨时间、跨传感器的定量对比，大气校正是无论如何都不能省的。

5. 将裁剪与大气校正整合成可复用流程

5.1 从原始影像到最终产品的一体化脚本

前面讲的是裁剪和大气校正分开做，但在真实项目中，这两步很少会独立操作。因为它们其实是“流水线”上的两个工位：影像进来，先按边界裁剪，裁完了拿到研究区再跑大气校正。所以我把整个流程写成一套一体化的工作流，保证中间产物能追溯、断了能续跑。

下面是一个整合入口脚本的简化版：

import os from clip import run_clip from atmospheric_correction import run_correction def main(): # 第一步：裁剪 run_clip(raw_dir='data/raw', clip_dir='data/output/clipped', roi_file='data/vector/roi.shp', resample_alg='bilinear') # 第二步：大气校正 run_correction(clip_dir='data/output/clipped', out_dir='data/output/corrected', meta_file='data/metadata.json') if __name__ == '__main__': main()

这样设计的最大好处是模块之间的边界清晰，裁剪如果出问题，不会影响大气校正部分的调试。另外，中间结果跟最终结果分开保存，一方面便于检查裁剪质量，另一方面，如果后续发现大气校正参数需要重调，直接基于裁剪输出重新跑校正步骤就行，不用再把原始大影像拖进来重新裁一遍。

5.2 日志与质量检查

在批处理中，我强烈建议每一步都写日志，记录每个文件的处理时间、波段数、尺寸、是否成功等信息。最简单的实现是用Python的logging模块，在函数关键节点输出INFO。

另一个重要的质量检查点是：裁剪后影像的尺寸和波段数与原始影像是否匹配。有些时候，由于矢量边界和影像只有很小一部分交集，裁剪出来的影像很小，需要在日志里给出“面积比例”提示。还可以计算裁剪前后影像的非零像元总数，对比是否有异常丢失。

5.3 性能优化：多进程批处理

如果数据量大到几十上百景，单线程循环处理每一景需要的时间就很可观了。此时可以用multiprocessing池来并行处理。

from multiprocessing import Pool import os from osgeo import gdal def process_one_scene(scene_name): in_path = os.path.join(raw_dir, scene_name) out_path = os.path.join(clip_dir, scene_name.replace('.tif', '_clip.tif')) gdal.Warp(out_path, in_path, cutlineDSName=roi_file, cropToCutline=True, dstNodata=0.0, options=['-overwrite']) return scene_name if __name__ == '__main__': scenes = [f for f in os.listdir(raw_dir) if f.endswith('.tif')] with Pool(processes=4) as pool: for res in pool.imap_unordered(process_one_scene, scenes): print(f"完成: {res}")

这里有几个要注意的地方：

• 不要在子进程里调用多线程的GDAL操作，否则有GIL或者底层锁冲突风险。每个子进程独立调用GDAL是安全的。
• 输出文件路径如果存在且被其他进程占用，会报错。所以一定要先建好目录，每个进程写不同文件时天然隔离。
• CPU核心数不建议直接拉满，我一般设为物理核心数减1，留出一个核心给IO和系统调度。

6. 常见问题与排查技巧实录

6.1 裁剪结果全黑/空白

这种现象最常见的两个原因：一是矢量边界和影像坐标系不一致，导致求交结果为空；二是cropToCutline设置不当，输出范围跑偏了。

排查顺序：

1. 检查两个文件的投影EPSG编码是否一致。
2. 打印矢量边界的范围（min_x, min_y, max_x, max_y）和影像的范围，看看是否有交集。
3. 在GIS软件或QGIS里叠加上去人工看一眼。

如果坐标系不一致，可以用gdal.VectorTranslate给矢量做重投影，或者直接在gdal.Warp里加cutlineSRS参数指定矢量坐标系的EPSG（较少用，一般直接重投影更保险）。

6.2 GDAL读取影像报'CPLE_OutOfMemory'

大影像一次性读进内存会爆内存。这个问题的根源是ReadAsArray()读入了整个波段。优化方案有两种：

• 使用GDAL的ReadRaster分块读取，一次读一个块，逐块处理再写回。
• 使用gdal.Warp配合-co TILED=YES创建tiled GeoTIFF，并让GDAL内部自己管理块缓存；裁剪时不必一次读完。

对于大气校正环节，因为要做逐像元计算，分块读取几乎是必须的。示例如下：

def block_process(in_path, out_path, block_size=1024): ds = gdal.Open(in_path) driver = gdal.GetDriverByName('GTiff') out_ds = driver.Create(out_path, ds.RasterXSize, ds.RasterYSize, ds.RasterCount, gdal.GDT_Float32) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) band = ds.GetRasterBand(1) out_band = out_ds.GetRasterBand(1) for y in range(0, ds.RasterYSize, block_size): rows = min(block_size, ds.RasterYSize - y) for x in range(0, ds.RasterXSize, block_size): cols = min(block_size, ds.RasterXSize - x) data = band.ReadAsArray(x, y, cols, rows) # 在这里做处理 processed = data * 0.1 out_band.WriteArray(processed, x, y) out_ds.FlushCache()

这个分块读写模式，是我在处理两三个GB级别的影像时最常用的手法，稳定不出错。

6.3 Py6S安装后运行报错

Py6S本身安装容易，但运行时如果提示找不到6S可执行文件，或者Fortran编译失败，绝大多数是环境变量问题。解决路径是：到Py6S的源码目录中找到6S的二进制或源码，手动编译后放到系统能识别的位置，或者设置环境变量指向该可执行文件。

Linux下需要gfortran，Windows下一般直接使用预编译包。如果实在编译不过去，可以尝试用conda环境安装，conda的预编译包在兼容性上通常会好一些。

6.4 大气校正之后数值不合理

校正后的反射率出现大于1或者大量负值，很多人第一反应是“模型参数错了”。我的排查经验是：

• 检查输入辐亮度是否有负值。如果原始DN值里有NoData被当成有效值参与定标，负值就是这么来的。
• 检查大气参数的季节匹配。中纬度夏天的标准大气剖面，如果用在冬天的高纬度数据，结果差得离谱。
• 检查影像是否已经做过辐射定标。如果输入的是DN值直接塞进6S公式，系数完全对不上，结果当然不对。

提示：遇到数值问题时，先在单波段单像元上手动算一遍，用计算器验算公式有没有写错，再排查影像层面。

6.5 裁剪后影像的地理定位偏移

裁剪后的影像在GIS里打开，和其他数据叠加对不齐，最常见的原因是重采样方法引起的像元位置偏移。可以换用最近邻重采样看看是否仍然偏移。如果偏移量为一个整像素，多半是参考网格对齐方式的问题；如果偏移量不规律，那要怀疑原始RPC或几何控制点信息本身有误差。

还有一种可能性是GDAL在裁剪时默认把输出网格对齐到原影像的角点，当你用cropToCutline去贴合矢量边界时，输出网格起点会变化，导致相对原影像产生不超过一个像素的平移。这个在遥感精度要求高的场景下需要额外注意，可以通过在gdal.Warp选项里指定targetAlignedPixels=True和te参数来手动控制输出网格。

6.6 常见问题速查表

7. 实操收尾的一点心得

最后聊几句个人的体会。在长期和数据打交道的过程里，我最大的感受是：遥感预处理这事，代码本身不难，难的是对数据质量的敏感度。比如看到一张裁剪后的影像，多了一个零值带，你会不会在第一眼就觉得“这数据哪里不对劲”？这种敏感度需要靠大量的实操、报错、对比来积累。

另一个体会是，工具链不要贪新鲜。GDAL这个库看起来古老，但它稳定、可靠、性能好，在绝大多数场景下足以应付。很多人一上来就上深度学习框架做影像分割，结果数据没洗干净，预处理环节问题一摞一摞，模型训练效果当然拉胯。把裁剪和大气校正这些“地基活”做扎实，后续才能高枕无忧。

还有一个小技巧：处理完一批数据后，一定挑一景影像，同时导出原始DN影像、表观反射率影像和校正后的反射率影像，放到GIS软件里直观对比一下。肉眼看出来“确实变清楚了”——这比任何日志和数值检查都更让人心里有底。多攒几组这样的对比样本，你对各种传感器数据会建立起非常直观的把握能力。