遥感毕设实战：基于Python与开源GIS工具链的端到端处理流程
摘要：许多遥感方向的本科毕设面临数据处理链路不清晰、工具碎片化、结果复现困难等问题。本文以典型土地利用分类任务为例，构建一套基于Python、GDAL、Rasterio和GeoPandas的轻量级遥感处理流水线，涵盖影像预处理、特征提取、模型推理与结果可视化。读者可直接复用该架构，显著降低环境配置成本，提升实验迭代效率，并确保处理流程的可追溯性与可复现性。
1. 毕设常见痛点：数据、坐标、内存三座大山
做遥感毕设，90% 的时间其实耗在“折腾数据”上。我当年踩过的坑，总结下来就三句话：
- 数据格式混乱：导师给的
.img、自己下载的.tif、同学分享的.dat，后缀五花八门，打开才发现波段顺序、位深、无效值全不一样。 - 坐标系不统一：WGS84、UTM、CGCS2000 混着来，QGIS 里看着“叠上了”，其实偏移几百米，一跑分类结果全是“NoData”。
- 内存溢出：笔记本 16 G 内存，ENVI 一打开 8 G 的 Sentinel-2 十波段影像直接卡死，只能降采样、裁切，信息损失心里没底。
传统桌面软件交互友好，但自动化=0；一旦要调参、批处理、写论文附图，鼠标点断手。于是我把目光投向 Python 开源栈：GDAL 负责底层驱动，Rasterio 给 GDAL 套上“Pythonic”外壳，GeoPandas 管矢量，scikit-learn 管分类，整套流程可脚本化、可版本控制、可一键复现。
2. 传统工具 vs Python 开源栈：谁更适合毕设迭代？
| 维度 | ENVI/QGIS 桌面端 | Python 开源栈 |
|---|---|---|
| 交互探索 | 拖图层即可看 | 需写代码 |
| 自动化批处理 | 手动或 ModelBuilder | 脚本/循环/调度 |
| 环境复现 | 换电脑就要重装插件 | 一条requirements.txt |
| 大文件内存策略 | 全盘加载 | 窗口读写、内存映射 |
| 结果可视化 | 出图快 | Matplotlib 可定制 |
结论：桌面端适合“看一眼”，Python 栈适合“跑一百遍”。毕设要调参、要对比实验、要出附录，后者效率翻倍。
3. 端到端代码实战：以“ Landsat-8 土地利用分类”为例
下面给出一条最小可用流水线（Clean Code 版），依赖如下：
rasterio==1.3.8 geopandas==0.13.2 scikit-learn==1.3.0 numpy==1.24.33.1 目录结构
proj/ ├─ data/ │ ├─ LC08_123032_20230809.tif # 原始影像 │ └─ roi.shp # 研究区矢量 ├─ out/ └─ src/ ├─ 01_preprocess.py ├─ 02_feature.py ├─ 03_train.py └─ 04_predict.py3.2 01_preprocess.py：裁剪+重投影+NDVI
import rasterio as rio from rasterio.mask import mask import geopandas as gpd import numpy as np from pathlib import Path SRC = Path('../data/LC08_123032_20220209.tif') ROI = Path('../data/roi.shp') DST = Path('../out/clip.tif') def read_roi(roi_path: Path): """返回与影像相同 crs 的 geojson 列表""" gdf = gpd.read_file(roi_path).to_crs(rio.open(SRC).crs) return [gdf.geometry.__geo_interface__] def clip_and_reproj(src_path, dst_path, shapes): """裁剪 + 重投影到 EPSG:32650，像素 30 m""" with rio.open(src_path) as src: clip, transform = mask(src, shapes, crop=True, indexes=src.indexes) meta = src.meta.copy() meta.update({ 'height': clip.shape[1], 'width': clip.shape[2], 'transform': transform, 'crs': 'EPSG:32650' }) # 写盘 with rio.open(dst_path, 'w', **meta) as dst: dst.write(clip) return dst_path if __name__ == '__main__': shapes = read_roi(ROI) clip_and_reproj(SRC, DST, shapes)3.3 02_feature.py：计算 NDVI 并堆叠波段
import rasterio as rio import numpy as np from pathlib import Path IN = Path('../out/clip.tif') OUT = Path('../out/stack.tif') def calc_ndvi(red, nir): """忽略无效值""" np.seterr(divide='ignore', invalid='ignore') return (nir - red) / (nir + red) def stack_with_ndvi(src_path, dst_path): with rio.open(src_path) as src: prof = src.profile prof.update(count=6, dtype=np.float32) red = src.read(4).astype('float32') nir = src.read(5).astype('float32') ndvi = calc_ndvi(red, nir) with rio.open(dst_path, 'w', **prof) as dst: for b in range(1, 6): dst.write(src.read(b), b) dst.write(ndvi, 6) # NDVI 放第 6 波段 return dst_path if __name__ == '__main__': stack_with_ndvi(IN, OUT)3.4 03_train.py：随机森林训练
import rasterio as rio import numpy as np from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import train_test_split import joblib IMG = '../out/stack.tif' LABEL = '../data/train_samples.shp' # 矢量样本，字段 'class' MODEL = '../out/rf.pkl' def sample_from_shape(img_path, label_path): """矢量→栅格采样，返回 X, y""" with rio.open(img_path) as src: gdf = gpd.read_file(label_path).to_crs(src.crs) X, y = [], [] for cls, geom in zip(gdf['class'], gdf.geometry): # 中心点采样 row, col = src.index(geom.x, geom.y) window = rio.windows.Window(col-1, row-1, 3, 3) patch = src.read(window=window).reshape(-1, 6) X.append(patch.mean(axis=0)) y.append(cls) return np.vstack(X), np.array(y) def train_and_save(X, y, model_path): X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42) clf = RandomForestClassifier(n_estimators=300, max_depth=15, n_jobs=-1, random_state=42) clf.fit(X_train, y_train) print('OA:', clf.score(X_test, y_test)) joblib.dump(clf, model_path) return clf if __name__ == '__main__': X, y = sample_from_shape(IMG, LABEL) train_and_save(X, y, MODEL)3.5 04_predict.py：整景推理 + GeoTIFF 输出
import rasterio as rio import numpy as np import joblib from tqdm import tqdm IMG = '../out/stack.tif' MODEL = '../out/rf.pkl' PRED = '../out/lc_2022.tif' def predict_blockwise(img_path, model_path, dst_path, bs=1024): """分块推理，防止内存爆炸""" clf = joblib.load(model_path) with rio.open(img_path) as src: prof = src.profile prof.update(count=1, dtype=np.uint8, nodata=255) h, w = src.height, src.width with rio.open(dst_path, 'w', **prof) as dst: for i in tqdm(range(0, h, bs)): for j in range(0, w, bs): window = rio.windows.Window(j, i, bs, bs) chunk = src.read(window=window) # reshape 像素×波段 pixels = chunk.reshape(6, -1).T pred = clf.predict(pixels).astype('uint8') # 写回 dst.write(pred.reshape(window.height, window.width), 1, window=window) return dst_path if __name__ == '__main__': predict_blockwise(IMG, MODEL, PRED)运行完四条脚本，即可在out/拿到：
clip.tif：研究区影像stack.tif：叠加 NDVI 的六波段数据rf.pkl：训练好的模型lc_2022.tif：分类结果，类别值 1~5，背景 255
4. 大影像内存管理：窗口、块、缓存
- 窗口读写：Rasterio 的
Window可以把 20 GB 影像切成 1024×1024 的小块，逐块读、逐块写，内存占用恒定在 100 M 左右。 - 块大小选择：不宜太大（爆内存），也不宜太小（频繁 I/O）。经验值：SSD 环境下 2048×2048 像素，机械盘 1024×1024。
- 并发风险：Python 的 GIL 让多线程假并行，真正并行用 multiprocessing 或 joblib.Parallel，但写同一张 GeoTIFF需加锁，否则波段交织。推荐“读-计算-写”分离：先写临时文件，最后
gdalbuildvrt拼回去。 - 数据类型：float64 比 uint16 体积翻倍，中间结果一律
float32足够。
5. 生产环境避坑指南
- CRS 一致性校验：每次
rio.open后断言src.crs == roi.crs，不等即to_crs()。 - 临时文件清理：用
tempfile.TemporaryDirectory包裹，异常退出也能自动清。 - 结果元数据保留：写盘时把
src.profile完整拷过来，再更新必要字段，确保pixel_size、nodata、transform与原始对齐。 - 无效值掩膜：分类前把
nodata区域做成掩膜数组，预测后贴回去，防止“背景被分到水体”。 - 版本锁定：
requirements.txt写死版本，Docker 镜像存一份，明年还能复现。 - 日志：每步输出
shape、dtype、memory_usage，调试用tqdm看进度，出问题能定位。
6. 拓展：把模型换成 XGBoost，把数据换成 Sentinel-2
整套流水线对输入影像的波段数、分辨率、传感器类型零硬编码。只要保持stack.tif的波段顺序与训练样本一致，把RandomForestClassifier换成XGBClassifier，或接入 Sentinel-2 的 10 m 分辨率数据，流程无需改动。想玩深度语义分割？把04_predict.py换成torch.inference，其余预处理、后处理模块可复用。
7. 小结
把毕设从“鼠标流”搬到“脚本流”，最大的收获是实验可迭代：调一个参数，十分钟出图；换一张影像，一条命令搞定。希望这条轻量级开源流水线，能让你把精力花在“算法改进”而不是“格式转换”。祝各位毕业顺利，代码可复现，论文好写，答辩不慌。
