本文分享一套经过多次迭代的Python 遥感批处理代码，一键搞定 HDF/Tiff 处理、矢量裁剪、多时相求均值、分区统计以及带地名标注的精美出图。
做遥感分析时，最耗时的往往不是核心算法，而是繁琐的数据预处理：
今天，我将分享一套“保姆级”的 Python 遥感批处理脚本。这套代码经过了多次实战调试，专门解决了上述所有痛点。
一、代码核心功能一览
这套代码集成了以下 5 大核心功能：
自动坐标对齐：不管你的影像和矢量是投影坐标还是经纬度，代码会自动检测并强制对齐，杜绝“图是图，线是线”的分离现象。
批量矢量裁剪：利用 Geopandas 和 Rasterio，自动按行政边界裁剪，并去除黑边。
多时相求均值：自动堆叠多期影像，忽略无效值（NoData/云），计算多年平均 NDVI。
分区统计 (Zonal Statistics)：一键统计每个区县的 Mean/Max/Min 值，并导出为 Excel/CSV 表格。
精美可视化：
○白底修正：自动将背景 0 值设为白色透明，不再是一片红。
○边界叠加：完美叠加行政区划黑框线。
○智能标注：自动读取 SHP 字段，在地图上标注地名。
二、核心技术点解析
在使用 Python (Matplotlib + Rasterio) 绘图时，最容易踩的坑就是坐标范围不匹配。
❌ 错误做法：
直接使用 plt.imshow(data)。这样画出来的图，坐标轴是像素行列号（例如 0~2000），而矢量文件是经纬度（例如 30°N~40°N）。两者不在一个量级，导致矢量边界线“隐身”。
✅ 正确做法 (本代码方案)：
我们引入了 rasterio.plot.plotting_extent，获取栅格数据的真实地理坐标范围，并传给 imshow：
from rasterio.plot import plotting_extent# 获取图片的地理坐标范围 [左, 右, 下, 上]raster_extent = plotting_extent(src)# 画图时指定 extent，强行对齐经纬度ax.imshow(data, extent=raster_extent, ...)
此外，为了让出图更美观，我们对地名标注做了“描边”处理，即使在深色植被背景下，文字依然清晰可见：
# 文字描边效果text.set_path_effects([PathEffects.withStroke(linewidth=2, foreground='white')])
如何使用这套代码？
使用非常简单，你只需要修改代码最上方的 CONFIG 配置区：
CONFIG = { # 1. 基础路径 "root_dir": r"D:\20260122", "input_folder": "GEE_MODIS_NDVI", "shp_folder": "shp", "shp_filename": "henan_city.shp", # 2. 地名标注 "label_field_name": "name", "label_font_size": 12, "show_labels": True, "show_boundaries": True, # 3. 数据参数 "nodata_value": -9999, "scale_factor": 1, # 4. 功能开关 (建议只开出图，省时间) "run_hdf_to_tif": False, "run_clip": False, # 裁剪过的不用再裁了 "run_mean": False, # 均值算过的不用再算了 "run_zonal": False, # 统计过的不用再统了 "run_plot": True, # <--- 只运行这个！ # 5. 绘图参数 "plot_cmap": "RdYlGn", "plot_vmin": 0.0, "plot_vmax": 0.9, "plot_dpi": 300}# 设置中文字体 plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False
三、成果展示
1. 完美的白底 NDVI 分布图（叠加市界和地名）：
2. 自动生成的统计表格：
def check_paths(): print("------ 正在检查文件路径 ------") if not os.path.exists(CONFIG["root_dir"]): print(f"[错误] 找不到根目录: {CONFIG['root_dir']}") return False shp_path = os.path.join(CONFIG["root_dir"], CONFIG["shp_folder"], CONFIG["shp_filename"]) if not os.path.exists(shp_path): print(f"[错误] 找不到shp文件: {shp_path}") return False print(f"[通过] 路径检查正常。") return Truedef make_output_dirs(base_path): dirs = { "tif": os.path.join(base_path, "out", "1_converted_tif"), "clip": os.path.join(base_path, "out", "2_clipped"), "mean": os.path.join(base_path, "out", "3_mean_result"), "stats": os.path.join(base_path, "out", "4_zonal_stats"), "plot": os.path.join(base_path, "out", "5_plots") } for d in dirs.values(): os.makedirs(d, exist_ok=True)
四、结语
由于代码较长（含详细注释和异常处理），微信公众号回复：C003，获取代码与数据。
如果你觉得这个脚本对你有帮助，欢迎点赞、在看、分享给更多被 GIS 手动操作折磨的朋友！
