Python与Landsat遥感影像的作物估产实战：从NDVI到随机森林

做农业遥感这几年，我越来越觉得，作物估产是一条完整的数据链路，而不是某个算法的独角戏。用Python玩转农业大数据，对我来说不是一句口号，而是每天都在做的事：从Landsat遥感影像中读取波段，预处理成干净的反射率数据，算出NDVI这类植被指数，再结合田间产量记录完成建模。这篇教程就把我最近半年反复跑通的一条工作流完整复盘出来，从Python环境搭建、Landsat数据下载，一直讲到作物产量预测模型怎么构建、怎么评估、怎么避开我在实操中踩过的坑。

如果你正准备做农业遥感，或者手里已经有一批田间产量数据但不知道怎么和卫星影像结合，这篇文章应该能帮你省下不少摸索时间。后面所有代码我都尽量给出能直接改参数运行的版本，处理对象是Landsat 8/9影像，用到的是Python生态里常见的rasterio、geopandas、rasterstats和scikit-learn。有一定Python基础最好，没有的话照着装环境、按步骤跑也能走通，只是遇到报错时需要多点耐心。

1. 项目概述与整体设计思路

1.1 为什么是Landsat影像和Python的组合

先回答一个很多人问过我的问题：做作物产量预测，为什么优先用Landsat，而不是Sentinel-2或者MODIS？

Landsat项目从1972年一直延续到现在，Landsat 8和Landsat 9两颗卫星目前在轨运行，每16天重访同一地点一次，两颗卫星交叉覆盖后大概8天就能拿到一景影像。它提供30米空间分辨率的红、绿、蓝、近红外、短波红外等波段，这个分辨率对地块尺度的农业监测来说非常合适：比MODIS的250米精细得多，能看清田块内部差异；又比Sentinel-2的10米数据更适合做长时序分析，因为Landsat的历史存档从80年代就开始了，可以反推过去几十年的植被变化。

更重要的是，Landsat数据是完全免费开放的。之前做项目时，动不动要买高分辨率商业影像，一景几万块钱，小课题组根本扛不住。Landsat的出现直接让农业遥感研究门槛降了一大截。配合Python，整个链条都能自动化：rasterio和GDAL处理栅格，geopandas处理矢量地块，numpy做数值计算，rasterstats做分区统计，scikit-learn做回归建模。这些库全部免费，而且社区资料极其丰富，遇到问题基本一搜就有答案。

我自己的体会是，30米分辨率的Landsat虽然不能精确到单株作物，但对于县域级、农场级的产量估测已经够用。尤其是当你的产量数据是地块尺度的、面积足够大时，Landsat影像的混合像元影响可以被有效摊薄。这也是我后来反复跟学生强调的：选数据不是越高分辨率越好，而是要和你的样本单元匹配。

1.2 完整技术路线：从原始影像到一亩地的估算产量

整个项目的技术路线可以拆成六个环节，每一步都有明确的产品输出：

1. 确定研究区域和时间窗口，比如“华北某县冬小麦，关键生育期是3月中旬到6月上旬”。
2. 下载研究区内的Landsat影像，优先选云量低于10%的景。
3. 做预处理：辐射定标、大气校正、云掩膜、影像裁剪。如果直接用Level-2表面反射率产品，前两步可以省掉。
4. 计算植被指数，最常用的是NDVI，再按生育期提取峰值、累积值等特征。
5. 把产量数据和遥感特征表对齐，整理成模型输入。
6. 训练回归模型，评估预测精度，输出产量图或区域总产量。

我见过不少新手一上来就想用深度学习做端到端预测，把整张影像直接扔给卷积神经网络。精神可嘉，但实际项目里通常走不通，因为农业产量数据样本量太少，一年只有几十个到几百个地块，撑不起一个大模型。更合理的思路是：先用遥感影像提取可解释的特征，比如NDVI、EVI这些和作物生长状态直接相关的量，再用随机森林这类对中小样本友好的模型做回归。这样既稳定，又能在出问题时追查到是哪个特征出了问题。

后面几节我会按这条路线逐步展开，每到一个环节都会给出能直接运行的Python代码，以及我在实践中踩过的坑。接下来先从环境搭建说起，环境装不对，后面所有代码都会在import那一步卡住。

2. 环境搭建与数据准备

2.1 Python开发环境配置：从源头减少折腾

很多人在环境配置阶段就放弃了，其实大部分遥感库安装报错，根源都是用了系统自带的Python和pip硬扛。我现在的标准做法是：装Miniconda，建独立环境，用conda-forge渠道装地理相关库。

步骤很简单，Windows、Linux、macOS都类似。先去Miniconda官网下载安装包，装完打开Anaconda Prompt或终端，执行：

conda create -n agro python=3.9 -y conda activate agro

这里建议用Python 3.9或3.10，目前rasterio、geopandas这些库在这些版本上兼容性最好。接下来安装核心库：

conda install -c conda-forge gdal rasterio shapely pyproj -y pip install geopandas pandas numpy scikit-learn matplotlib jupyter

为什么要先用conda装GDAL相关库？因为GDAL是很多遥感库的底层依赖，如果用pip直接装，经常会出现“gdal-config: command not found”或者版本冲突。conda-forge仓库会把gdal、rasterio、shapely这些打包成一套兼容版本，安装体验顺滑很多。

编辑器方面，我推荐VSCode，免费、插件丰富、启动快。装好Python扩展后，按Ctrl+Shift+P，选择“Python: Select Interpreter”，找到名为agro的conda环境，再新建终端，确认命令行左端显示(agro)，就说明环境激活了。这一步经常有人忽略，结果明明装好了库，运行时还是提示ModuleNotFoundError，其实所有库都在conda环境里，但编辑器默认用的却是系统Python。

2.2 Landsat数据下载与ArcGIS Pro加载查看

数据下载是第一步，也是最容易被低估的一步。我常用的渠道是USGS EarthExplorer，注册账号后按路径选择下载。注意几个筛选条件：云量阈值设在10%以下，优先选Landsat 8/9的Collection 2 Level-2产品，时间范围覆盖作物关键生育期。

Landsat的命名规则乍一看很唬人，比如“LC08_L2SP_119039_20230515_20230518_02_T1_SR_B4.TIF”，拆开看其实很清楚：LC08表示Landsat 8，L2SP表示Level-2表面反射率产品，119039是行列号，20230515是采集日期，SR表示表面反射率，B4是波段编号。Level-2产品的好处是已经完成辐射定标和大气校正，拿到手就能直接用，能省掉预处理里最麻烦的一段。

下载后，我习惯先拿到ArcGIS Pro里快速看一眼。新建地图，在“地图”选项卡下拉菜单里选择“添加数据”，选中下载的GeoTIFF文件。单波段影像默认显示成灰度，想要看得舒服，用“栅格图层”的“符号系统”切换波段组合。Landsat 8/9的假彩色影像波段组合是红=5、绿=4、蓝=3，这个组合下植被呈红色，水体和建筑一目了然，非常适合检查影像质量和云覆盖情况。

如果你的项目需要按地块统计，那就得先有地块边界。在ArcGIS Pro里创建新的地类图斑矢量图层，步骤是：在“目录”窗格右键数据库或文件夹，新建面要素类；坐标系一定要选和影像一致，比如UTM投影带或者WGS 1984地理坐标系，后面用Python合并时才能对得上。建好后切换到“编辑”选项卡，用“创建要素”工具沿着田块边界绘制多边形，避免让多边形跨过道路、水塘这类非农用地。属性表里添加几个字段：field_id、crop_type、plant_date、yield_kg_ha，这些字段后面会在建模环节派上大用场。画完后导出成Shapefile或GeoPackage，放到项目文件夹下。

2.3 产量样本数据准备：决定上限的一环

遥感影像再干净，没有对应的产量数据，模型就无从训练。产量数据的来源可以是测产报表、联合收割机产量图、合作社记录或者统计年鉴，但我强烈建议至少有一年内地块尺度的实测数据，否则模型只能学到平均值。

产量数据整理有两个关键点。第一，单位必须统一，我遇到过好几次地块面积用的是亩，产量单位却是吨，一算就差了15倍，模型输出完全没法解释。建议全部换算成kg/ha，即每公顷公斤数。第二，坐标必须准确，尽量记录每个地块的中心经纬度或者有对应边界，这样和遥感像元匹配时才能落在正确的位置。样本表大概长这样：

还有一个空间尺度匹配问题。Landsat像元是30米×30米，如果地块小于一个像元，或者形状特别狭长，提取到的像元值会掺入大量非植被信息。我一般会把面积小于0.2公顷的地块直接剔除，或者在统计时只取地块边界内部的核心像元，避开边缘混合区。宁可样本少一点，也好过噪声大一点。

3. 遥感影像预处理全流程

3.1 辐射定标与大气校正：能省则省，但有前提

在Landsat早期产品里，影像记录的是DN值（无量纲的数字量化值），它受太阳高度角、传感器增益、大气散射等多个因素影响，不能直接用来比较不同日期、不同区域的植被状态。辐射定标是把DN值转换成大气顶反射率，大气校正再进一步消除大气分子和气溶胶的影响，得到地表反射率。

我们为什么在乎地表反射率？因为NDVI这类植被指数，本质上就是红光和近红外反射率的比值运算。如果不做大气校正，同一个地块在不同日期算出来的NDVI会忽高忽低，完全没法用于生育期对比。

不过，使用Collection 2 Level-2产品时，数据已经完成了表面反射率转换，这一步可以跳过。我现在的项目基本全部使用L2产品，省下的时间用来处理更让人头疼的云掩膜。你只需要处理QA_PIXEL波段，按位运算提取云和云阴影。

import rasterio import numpy as np qa_path = 'LC08_L2SP_119039_20230515_20230518_02_T1_QA_PIXEL.TIF' with rasterio.open(qa_path) as qa: qa_arr = qa.read(1) profile = qa.profile # Landsat Collection 2 QA波段：bit 4是云，bit 5是云阴影 cloud_mask = ((qa_arr & (1 << 4)) > 0) | ((qa_arr & (1 << 5)) > 0)

这里用了按位与运算。QA波段每个像元的二进制位记录了不同质量标志，比如第4位是云，第5位是云阴影。理解了这个逻辑，就能灵活提取各种掩膜。提取后，把云和云阴影像元在后续的NDVI计算中设为无效值，能显著减少异常峰值的干扰。

3.2 影像裁剪、拼接与NDVI计算

拿到若干景Landsat影像后，先检查投影坐标系和范围。如果研究区横跨多景影像，可以用rasterio.merge或gdal.Warp做镶嵌；如果只是关心某个县域，建议用研究区矢量直接裁剪，这样数据量小、处理快。

裁剪我习惯用rioxarray，代码逻辑很清晰：

import rioxarray from shapely.geometry import box import geopandas as gpd gdf = gpd.read_file('study_area.shp') # 确保矢量与栅格坐标系统一 gdf = gdf.to_crs(epsg=32649) xds = rioxarray.open_rasterio('LC08_L2SP_119039_20230515_20230518_02_T1_SR_B4.TIF') xds = xds.rio.clip(gdf.geometry, gdf.crs, drop=True)

然后就是NDVI计算。Landsat 8/9波段编号中，B4是红光波段，B5是近红外波段，NDVI公式为：

NDVI = (NIR - RED) / (NIR + RED)

import rasterio import numpy as np red_path = 'LC08_L2SP_119039_20230515_20230518_02_T1_SR_B4.TIF' nir_path = 'LC08_L2SP_119039_20230515_20230518_02_T1_SR_B5.TIF' with rasterio.open(red_path) as src_r: red = src_r.read(1).astype('float32') profile = src_r.profile with rasterio.open(nir_path) as src_n: nir = src_n.read(1).astype('float32') ndvi = (nir - red) / (nir + red + 1e-10) ndvi = np.clip(ndvi, -1, 1) profile.update(dtype=rasterio.float32, count=1, compress='lzw') with rasterio.open('ndvi_20230515.tif', 'w', **profile) as dst: dst.write(ndvi, 1)

代码里加1e-10是为了防止分母为0，clip是把数值限定在-1到1之间。这一步虽小，但如果不处理，个别像元会出现几百这种离谱值，后面统计特征时会把均值带歪。

有了单期NDVI，还要做多时间的特征合成。比如提取研究区整个生育期的NDVI时间序列，按月排序后计算每个像元的峰值、累计值、平均值。我通常在生成完各期NDVI后，用numpy把它们堆叠成三维数组：

ndvi_stack = np.stack([ndvi_20230315, ndvi_20230410, ndvi_20230505, ndvi_20230525], axis=0) peak_ndvi = np.nanmax(ndvi_stack, axis=0) cum_ndvi = np.nansum(ndvi_stack, axis=0)

峰值NDVI对应作物生长最旺盛期的叶面积指数，累积NDVI则反映整个生育期的光合累积量，这两个特征在产量预测中通常贡献最大。

3.3 从栅格提取地块级特征：把影像变成表格

模型吃的是表格，不是图片。所以下一步要把NDVI栅格和地块矢量结合起来，计算每个地块内的统计值。这个环节我强烈推荐rasterstats库，一行代码就能完成分区统计：

import geopandas as gpd import pandas as pd from rasterstats import zonal_stats gdf = gpd.read_file('sample_fields.shp') # 确保矢量和栅格都在同一投影坐标系 gdf = gdf.to_crs(epsg=32649) stats = zonal_stats(gdf, 'ndvi_20230515.tif', stats=['mean', 'max', 'min', 'median', 'count']) gdf = gdf.join(pd.DataFrame(stats))

zonal_stats返回的count是每个地块内有效像元数量，这个字段很重要。如果count太小，说明地块面积和30米分辨率不匹配，统计结果大概率受混合像元影响，考虑剔除。聚合时，我习惯同时保留median和mean，因为NDVI分布通常偏态，像元中位数比均值更稳健，能减少灌溉沟渠、田埂造成的高值像元干扰。

把多个时期的NDVI分区统计后，再和产量表left join，就得到了一个完整的建模输入表。每一行是一个地块，列是field_id、year、peak_ndvi、cum_ndvi、median_ndvi、yield_kg_ha这些字段。到这一步，你的数据已经从“一堆卫星影像”变成了“一张可以交给机器学习模型训练的表”，接下来才是建模的活。

4. 产量预测模型构建与评估

4.1 特征工程：先想清楚产量被什么决定

很多新手拿到那张表后第一件事就是打开sklearn跑线性回归，结果往往不理想。问题通常出在特征上：只用单期NDVI来预测最终产量，信息量太单薄。

作物的产量形成是一个累积过程，取决于播种后的出苗率、分蘖数、最大叶面积指数、灌浆持续期、收获指数等多个因素。因此我建议至少提取以下几类特征：

如果数据条件允许，还可以加入播种日期、品种类型、土壤有机质含量这些非遥感信息。实测经验告诉我，把气候特征加进模型后，R²在大多数年份能提升0.05到0.1，尤其当研究区内降雨年份间差异明显时。

提取完特征后，先做一次简单的相关性检查：

import pandas as pd df = pd.read_csv('sample_data.csv') print(df.corr()['yield_kg_ha'].sort_values(ascending=False))

如果发现peak_ndvi和cum_ndvi相关性超过0.85，没必要同时保留高相关特征，模型会变得不稳定。我一般会先用相关性矩阵筛一遍，再保留5到8个相对独立的特征，样本量小的时候尤其如此。

4.2 模型选型：线性回归、随机森林与XGBoost的取舍

农业产量和遥感特征之间通常是弱线性关系，更准确地说存在显著的非线性：长势太差产量低，长势中等不错，长势过高反而可能因为倒伏或病害导致减产。所以纯线性回归做基线可以，但上限有限。

我在项目中对比过三种常用模型：

产量数据经常只有几十到几百个样本，随机森林在这类场景里表现稳定，不需要做太多数据标准化，对异常值也不敏感，所以我把它作为默认选择。完整建模代码如下：

import pandas as pd from sklearn.model_selection import train_test_split, GridSearchCV from sklearn.ensemble import RandomForestRegressor from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error from sklearn.pipeline import Pipeline from sklearn.preprocessing import StandardScaler df = pd.read_csv('sample_data.csv') X = df.drop(columns=['field_id', 'year', 'yield_kg_ha']) y = df['yield_kg_ha'] X_train, X_test, y_train, y_test = train_test_split( X, y, test_size=0.2, random_state=42 ) param_grid = { 'n_estimators': [100, 200, 300], 'max_depth': [5, 8, 10, None], 'min_samples_leaf': [1, 2, 4] } rf = RandomForestRegressor(random_state=42) grid = GridSearchCV(rf, param_grid, cv=5, scoring='r2', n_jobs=-1) grid.fit(X_train, y_train) best_rf = grid.best_estimator_ y_pred = best_rf.predict(X_test) print('R2:', r2_score(y_test, y_pred)) print('RMSE:', mean_squared_error(y_test, y_pred, squared=False)) print('MAE:', mean_absolute_error(y_test, y_pred)) print('Best params:', grid.best_params_)

这里的GridSearchCV是在训练集内部做五折交叉验证选择参数，再用选出的模型在完全没有参与训练的测试集上做最终评估。如果R²在你自己的测试集上明显低于训练集，说明过拟合，需要减少max_depth或者把min_samples_leaf调大。

有一点必须提醒：如果你按地块随机划分train_test_split，同一年的多个相邻地块很可能同时出现在训练集和测试集，由于它们共享相似气候条件，验证指标会偏乐观。更严格的做法是按年份划分：

train_df = df[df['year'].isin([2019, 2020, 2021])] test_df = df[df['year'] == 2022] X_train = train_df.drop(columns=['field_id', 'year', 'yield_kg_ha']) y_train = train_df['yield_kg_ha'] X_test = test_df.drop(columns=['field_id', 'year', 'yield_kg_ha']) y_test = test_df['yield_kg_ha']

这种时间验证策略才更接近真实业务：我们永远是用过去几年的数据训练模型，去预测还没到来的一季产量。按年份划分后R²通常会比随机划分低一点，但这才是模型实际部署时的真实水平。

4.3 结果解读：真正决定产量精度的往往不是算法

模型训练完成后的结果解读，往往比调参更能说明问题。以我常跑的冬小麦数据为例，随机森林在测试集上的R²能做到0.70到0.80，RMSE在300到400 kg/ha之间。对于区域尺度估产来说，这个精度已经具备业务参考价值。

特征重要性输出非常简单：

importance = pd.Series(best_rf.feature_importances_, index=X.columns) importance.sort_values(ascending=False, inplace=True) print(importance)

在大多数年份里，peak_ndvi和cum_ndvi的重要性排在最前，这很符合直觉：峰值NDVI对应扬花灌浆期的最大叶面积，直接决定光合产物总量；累积NDVI则把全生育期的长势累积下来，和最终生物量的关系更密切。如果发现某个年份是最晚一期的平均值占主导，就要警惕是不是影像时间没有覆盖到生长高峰期，导致特征没有提取到位。

误差从哪里来？我总结下来主要有五个来源：第一，物候期不一致，同一张影像里有的地块刚返青，有的已经抽穗，NDVI对比失去意义；第二，气象灾害，比如灌浆期突然来了高温或倒伏，遥感特征还没反映出来，产量已经掉下去了；第三，田间管理差异，同一个地块不同区域施肥量、灌溉量不同，30米像元内混在一起；第四，产量数据本身的记录误差；第五，云掩膜造成的某些日期NDVI缺失。

理解这些误差来源的意义在于：模型预测结果不是终点，交叉检验之后还要回到地块上去看异常样本。我会把预测误差最大的前10个地块单独提出来，叠加到影像上看是否存在边界与道路吻合、是否遭遇过药害或者涝渍。这一步虽然费时间，但能让你对这套方法到底在哪些条件下可靠形成非常清楚的认识。

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

5.1 Python环境与库安装问题

这条流程里最容易劝退新手的，就是import失败。我整理几个高频问题：

特别强调一点：不要混着用conda和pip反复安装同一个库。我在项目早期有一回先用conda装了rasterio，后来又用pip装了一个更高版本，把底层libtiff搞冲突了，连读影像都会崩溃。后来我统一原则：基础地理库用conda-forge装，纯Python库用pip装，装完就不动版本。

5.2 影像处理与数据准备中的坑

先从最伤脑筋的波段编号说起。Landsat 5/7的红波段是B3，近红外是B4；而Landsat 8/9的红波段是B4，近红外是B5。如果方案文档沿用旧数据集的编号，你的NDVI很可能算反，植被区域直接变成负值。我在代码里会写一个变量保存传感器类型和波段映射，所有用波段的地方统一走映射表。

投影不一致是我遇到第二多的坑。ArcGIS Pro里新创建图斑时默认可能用地理坐标系，而Landsat影像通常是UTM投影；用rasterstats做分区统计时，除非栅格和矢量坐标系一致，否则系统会重新投影，轻则慢，重则错位。我的习惯是一进项目就把所有数据统一转成影像的UTM坐标参考，后续所有处理都在这个坐标参考下进行。

NoData的处理也容易翻车。如果影像中有NoData值，直接参与NDVI计算会产生NaN，统计特征时如果没跳过NaN，一个地块的均值直接变成空值。建议在读取波段时就检查并设置有效值范围，计算NDVI前把无效像元设成np.nan，后续统计时统一用nanmean和nanmax。

5.3 产量建模中的数据陷阱

遥感估产模型精度不高，经常不是算法选错，而是数据构造出了问题。我复盘自己的历史项目，发现最典型的几个陷阱：

• 样本量太少。几十个地块做train_test_split，测试集可能只有十来个样本，R²波动极大。这时改用留一法交叉验证或按年份划分会更诚实。
• 空间自相关。相邻地块长相相似，天然属于相近的产量水平；随机划分时模型相当于“偷看”了邻居答案。空间分块交叉验证能缓解，但在小范围研究区很难完全避免，所以结果外推要保守。
• 产量记录误差。田间测产通常会有5%到10%的误差，如果个别地块记录了异常值（比如把亩产850公斤误写成8500公斤），模型会被拉偏。建模前一定要画箱线图看y分布，超过3倍四分位距的样本直接剔除或回查记录。
• 时间错位。有次一个学生拿到的产量数据是前一年的，但影像下载的是当年，他还没发现就直接建模，结果R²为负。产出数据的时间字段必须和影像时间严格对应，这个错误看起来低级，但在数据多线流转时很容易发生。

这些坑都踩过一遍之后，我对模型结果的判断会多一分谨慎。如果某一次预测结果显示全县产量比上年增长20%，我不会直接出报告，而是先找产量最高的几个地块和高低产年份的NDVI分布做对比，确认不是特征计算错误导致的虚假信号。

最后说一点个人体会。跑通这套从Landsat影像到作物产量预测的流程之后，我最大的感受是：遥感估产工程里，80%的时间花在数据获取、预处理和样本整理上，真正调模型只占很小一块。如果你手头正好有某个地区的地块产量数据，建议先挑一个小范围、两三年的样本，把从影像下载到建模评估的链路完整跑一遍，再慢慢扩大覆盖范围。这样即使中途报错，你也能更快定位是数据问题、坐标问题还是模型参数问题，而不是对着几百个地块的数据一头雾水。希望这篇教程能帮你少走点弯路，把更多精力花在实际问题上。