搞遥感土地利用分析的同行应该都有这种体验手动下载影像、裁剪、掩膜、统计、出图一套流程走下来单一年份还行一旦涉及十几年、几十个年份、十几个分类类别纯手点鼠标能把人累到怀疑人生。我这次把 QGIS 和 PyQGIS 组合起来做了一套针对 MapBiomas 数据的自动化分析脚本从数据下载到分类统计再到专题图布局导出完全用代码驱动。这篇文章就把整套思路、代码细节、踩过的坑全部拆开讲一遍。这套内容适合在 QGIS 里做长时间序列土地覆盖分析的同行参考也适合刚接触 PyQGIS、想把重复操作脚本化的初学者。只要你愿意花一下午把脚本跑通后面再遇到类似的多期栅格分析任务基本都能复用这个框架。1. 需求拆解为什么选 PyQGIS 而不是其他方案1.1 核心需求与选型逻辑这个任务表面上是“分析 MapBiomas 数据”但实际拆开看它是一整套数据工程远程下载数据、按研究区裁剪、统一坐标系、分类统计、时间序列计算、批量出图。全部手工操作不仅慢而且很难保证每一年的处理参数完全一致。不同的重采样方式、不同的裁剪边界、不同的 NoData 处理都会让时间序列分析产生人为偏差。那为什么选 PyQGIS先说结论因为处理链里的核心算法裁剪、重分类、统计、布局出图都是 QGIS 生态里现成能用的而且 PyQGIS 可以直接调用 QGIS 的处理框架不需要自己用 GDAL 从零实现一遍。相比纯 GDAL 脚本PyQGIS 的代码更短出图和 QGIS 项目文件的交互也更方便相比 ArcPy它完全开源免费没有授权问题而且处理结果和你在 QGIS 界面里看到的是一模一样的逻辑。实际做的时候还要注意PyQGIS 有两种运行方式一种是在 QGIS 内置的 Python 控制台里跑另一种是在外部 IDE 里初始化 QGIS 环境跑独立脚本。我这次用的是后者因为要批量处理多年的数据独立脚本更适合自动化还能丢进计划任务定时执行。1.2 先搞清楚 MapBiomas 数据结构MapBiomas 的数据本质上是逐年覆盖的土地利用分类栅格分辨率 30 米每个像元的值代表一个分类代码。比如代码 1 是森林代码 3 是非森林自然植被代码 4 是农业用地代码 5 是非植被区域代码 6 是水体。分析的核心动作就是把研究范围内每个像元的分类代码统计一遍算出各类别面积再看年份之间的变化趋势。这里有个关键认知MapBiomas 的数据集是按瓦片tile组织的不是一次性给一个全国大文件。下载前必须知道研究区落在哪些瓦片范围内。你可以先下载瓦片索引矢量在 QGIS 里加载后用“按位置选择”找出与研究区相交的瓦片编号再去拼 URL 下载。这个步骤看着不起眼但能帮你避开“把整个数据集都下下来”的尴尬。1.3 工作流整体设计整套流程我分成四段第一段是数据获取按研究区边界定位瓦片批量下载对应年份的 GeoTIFF第二段是预处理把矢量边界和栅格统一到同一个坐标系按掩膜批量裁剪顺便处理 NoData第三段是统计计算读取每个年份栅格的像元值用分类代码做唯一值统计换算成面积输出 CSV 表格第四段是可视化一方面用 Matplotlib 画各类别面积的时间序列折线图另一方面用 QgsLayout 排版输出带图例、比例尺、指北针的专题图。这个顺序是固定的因为每一步的输出都是下一步的输入。整体脚本跑完一圈之后你再回看 QGIS 界面会发现它其实就是在替你做那些原本要在菜单里点来点去的事。2. 环境准备与数据下载把地基打牢2.1 初始化无界面 PyQGIS 环境先在外部 IDE 跑 QGIS 脚本有一个重要前提要让 Python 找到 QGIS 的库文件并初始化应用对象。很多初学者第一次跑脚本报错“ModuleNotFoundError: No module named qgis”就是因为没有把 QGIS 的 Python 路径加到 sys.path 里。import sys from qgis.core import QgsApplication # 根据你自己的 QGIS 安装位置调整路径 sys.path.append(/usr/share/qgis/python) QgsApplication.setPrefixPath(/usr/bin/qgis, True) # False 表示不带界面启动 qgs QgsApplication([], False) qgs.initQgis()初始化之后记得在脚本最后调用qgs.exitQgis()释放资源。这个细节容易忘但如果脚本里跑大量循环不释放会导致内存持续占用尤其在 Windows 上常见。提示setPrefixPath的参数因操作系统和 QGIS 版本略有差异。Linux 上一般是/usr/bin/qgisWindows 上通常是 QGIS 安装目录比如C:/Program Files/QGIS 3.28/apps/qgis。可以通过在 QGIS 控制台里执行QgsApplication.prefixPath()来确认当前环境。2.2 用代码按瓦片下载栅格数据下载数据前先准备一份研究区边界矢量。这个矢量可以是任何你手头已有的面图层比如某流域边界、某县界。关键是它的坐标系要明确后续统一坐标参考时要用。下载代码的逻辑是先读取研究区范围计算和哪些瓦片相交然后拼接出下载地址逐块保存。下面是一个简洁可用的版本import requests from qgis.core import QgsVectorLayer, QgsGeometry def get_tiles_for_region(region_path, tile_index_path): region QgsVectorLayer(region_path, region, ogr) tile_index QgsVectorLayer(tile_index_path, tiles, ogr) matched [] for feat in tile_index.getFeatures(): if feat.geometry().intersects(region.extent()): matched.append(feat[tile_id]) return matched def download_tiles(tile_ids, year, output_dir): for tile_id in tile_ids: # 实际地址按 MapBiomas 官方数据服务的规则拼接 url fhttps://example-data-server/{year}/{tile_id}.tif local_path f{output_dir}/{year}_{tile_id}.tif with requests.get(url, streamTrue) as r: r.raise_for_status() with open(local_path, wb) as f: for chunk in r.iter_content(chunk_size8192): f.write(chunk)注意几个细节第一下载地址必须按时效调整MapBiomas 的数据服务地址可能随着 Collection 版本变化第二建议加timeout参数和重试机制以防止网络抖动导致下载中断第三下载完成后顺手检查文件大小如果远远小于正常值大概率是下载到了错误页面。2.3 统一坐标参考系别在这里省时间MapBiomas 原始数据默认是 WGS84EPSG:4326但很多研究区边界是 UTM 投影或者地方坐标系。如果直接在原始坐标系上裁剪统计面积计算会有偏差尤其是高纬度地区。我的做法是在矢量层加载后先用QgsCoordinateReferenceSystem判断它的坐标系统如果和栅格不一致就用QgsCoordinateTransform把矢量转成栅格的坐标系。这样在裁剪时边界能精确对齐后面统计的面积才是真实的地面面积。代码示例from qgis.core import QgsCoordinateReferenceSystem, QgsCoordinateTransform, QgsProject source_crs QgsCoordinateReferenceSystem(EPSG:4326) target_crs QgsCoordinateReferenceSystem(EPSG:32723) # 根据实际需要换成 UTM 带号 transform QgsCoordinateTransform(source_crs, target_crs, QgsProject.instance())如果你的研究区范围跨多个 UTM 分带最稳妥的方案是保留 WGS84 地理坐标做裁剪统计因为栅格像元的面积在 30 米分辨率尺度上经纬度导致的面积误差通常可以接受。但如果做的是精确的碳储量估算建议还是投影到等面积投影上比如 Albers 等积投影。3. 自动分析核心从裁剪到统计全流程3.1 用 QgsProcessing 批量裁剪研究区裁剪是整个流程里最费操作的一步。QGIS 的“按掩膜裁剪栅格”工具在处理 NoData 和边缘像元时做了很多底层优化直接用processing.run()调用它比自己用 GDAL 写裁剪稳得多。import processing from qgis.core import QgsProcessingFeedback def clip_raster(raster_path, mask_layer, output_path): params { INPUT: raster_path, MASK: mask_layer, TARGET_CRS: EPSG:4326, NODATA: -9999, OUTPUT: output_path } result processing.run(native:cliprasterbymasklayer, params, feedbackQgsProcessingFeedback()) return result[OUTPUT]这里TARGET_CRS可以根据研究区需要改成任意投影坐标系。NODATA设置为 -9999是为了在栅格转数组后能统一识别空值。注意MapBiomas 原始数据里 NoData 一般是 0但裁剪后边缘会出现新的 NoData如果你不显式设置统计时就会把边缘空值当成真实类别导致面积虚高。循环批量处理多年份时把每年的输出文件按年份单独建目录存放命名规则统一成{year}_clip.tif。这样后面的统计和出图都能用通配符扫描不需要手工维护文件列表。3.2 读栅格数组按分类代码统计面积裁剪完成后下一步是把每个年份的栅格读入内存做分类统计。这里我用的是 GDAL 的 Python 绑定而不是 QGIS 的栅格接口原因是 NumPy 对多维数组的运算效率明显更高。from osgeo import gdal import numpy as np def calc_area_by_class(raster_path, pixel_size30): ds gdal.Open(raster_path) band ds.GetRasterBand(1) arr band.ReadAsArray() # 去掉 NoData只统计有效像元 arr_valid arr[arr 0] classes, counts np.unique(arr_valid, return_countsTrue) cell_area pixel_size * pixel_size # 平方米 area_km2 counts * cell_area / 1e6 return {int(c): float(a) for c, a in zip(classes, area_km2)}这里有个容易被忽略的坑如果你把裁剪输出设置为某个投影坐标系那像元的实际地面尺寸不再是 30 米乘 30 米而是由投影决定的单位。要正确计算面积最好从 GDAL 数据集里读取GeoTransform或者直接用ds.GetGeoTransform()拿到像元尺寸而不是硬编码 30 米。gt ds.GetGeoTransform() pixel_size abs(gt[1] * gt[5]) # 在 UTM 投影里约等于地面宽度实测下来在 UTM 投影下像素尺寸基本接近 30 米但用代码动态读取比写死更稳。3.3 生成时间序列表格与趋势曲线统计完所有年份后你会得到一张“年份 × 分类”的矩阵。为了方便后续分析和别人复现建议直接导出 CSV编码用utf-8-sig否则 Excel 打开中文表头会乱码。import csv csv_path landcover_trend.csv with open(csv_path, w, newline, encodingutf-8-sig) as f: writer csv.writer(f) writer.writerow([year, class_code, area_km2]) for year, class_dict in yearly_results.items(): for cls, area in sorted(class_dict.items()): writer.writerow([year, cls, area])做趋势图我习惯用 Matplotlib。图不用太花哨重点是显示每一类别面积随年份的增减。画折线图时把森林、农业、水体等主要类别分别用不同颜色区分Y 轴单位统一为平方千米。import matplotlib.pyplot as plt fig, ax plt.subplots(figsize(10, 6)) for cls, color in class_colors.items(): years sorted(yearly_results.keys()) areas [yearly_results[y].get(cls, 0) for y in years] ax.plot(years, areas, labelfClass {cls}, colorcolor, markero) ax.set_xlabel(Year) ax.set_ylabel(Area (km2)) ax.legend() plt.tight_layout() plt.savefig(landcover_trend.png, dpi300)这里可以再扩展一步不仅画绝对面积还可以画各类别占比的堆叠面积图看看研究区在十几年的尺度上土地利用结构怎么演化。很多论文里那种漂亮的年代际变化图其实就是这一步做出来的。4. 可视化落地从配色到出图全细节4.1 分类栅格配色与透明化设置在 QGIS 里手动设置分类配色很简单但脚本里要自动做就得把每个分类代码的符号化规则写清楚。MapBiomas 官方有一套推荐配色比如森林用绿色系农业用黄色系水体用蓝色系。这样出图之后读者一眼就能看懂类别不需要反复看数值。用 PyQGIS 给栅格图层设置样式标准做法是构造一个QgsSingleBandPseudocolorRendererfrom qgis.core import (QgsRasterLayer, QgsSingleBandPseudocolorRenderer, QgsColorRampShader, QgsRasterShader) def style_raster(layer): fcn QgsColorRampShader() fcn.setColorRampType(QgsColorRampShader.Interpolated) # 每个分类代码对应一个颜色 items [ QgsColorRampShader.ColorRampItem(1, QColor(#1a5b2d), Forest), QgsColorRampShader.ColorRampItem(3, QColor(#9ac3a2), Non-forest natural), QgsColorRampShader.ColorRampItem(4, QColor(#f5b45c), Agriculture), QgsColorRampShader.ColorRampItem(5, QColor(#d4271e), Non-vegetated), QgsColorRampShader.ColorRampItem(6, QColor(#2532e0), Water), ] fcn.setColorRampItemList(items) renderer QgsSingleBandPseudocolorRenderer(layer.dataProvider(), 1) renderer.setClassificationMin(1) renderer.setClassificationMax(6) renderer.setShader(fcn) layer.setRenderer(renderer) layer.triggerRepaint()设置好样式后别忘了给 NoData 设置透明色否则黑边会把图面弄得很难看。可以调用layer.dataProvider().setNoDataValue(1, 0)或者直接在渲染器里过滤掉小于等于 0 的像元。4.2 用 QgsLayout 拼装专题图QGIS 的打印布局功能非常强大脚本里可以通过QgsPrintLayout创建布局添加地图、图例、比例尺、指北针然后整体导出 PNG。虽然第一次写代码会觉得 API 冗长但好处是只要写一次后面换年份换区域都能自动复用。from qgis.core import (QgsPrintLayout, QgsLayoutItemMap, QgsLayoutItemLegend, QgsLayoutItemScaleBar, QgsLayoutExporter, QgsProject, QgsLayoutPoint, QgsLayoutSize) project QgsProject.instance() layout QgsPrintLayout(project) layout.initializeDefaults() layout.setName(analysis_map) # 添加地图部件 item_map QgsLayoutItemMap(layout) item_map.setRect(QRectF(20, 20, 160, 120)) layout.addLayoutItem(item_map) # 添加图例 legend QgsLayoutItemLegend(layout) legend.setLinkedMap(item_map) legend.setTitle(Land Cover) layout.addLayoutItem(legend)导出的关键参数是QgsLayoutExporter.exportToImage()。可以设置分辨率 dpi 和背景色。实际测试中300 dpi 的 PNG 已经足够满足论文插图或汇报展示需求。4.3 批量生成多年对比图如果你想把每一年都出一张专题图放进 GIF 里做成动态变化效果可以在循环里逐年份执行上述步骤。每轮循环改一下栅格图层的文件路径刷新项目注册然后重新导出。注意布局部件不需要重新创建只需要更新 map 部件中的图层即可。这一步最值得注意的性能问题是如果一年生成一张图十几个年份就要构建十几次布局。QGIS 的布局对象非常吃内存循环里最好每处理完一年就把布局对象删除或者只保留一个布局实例反复更新地图范围。我在实际操作时是每轮循环重建布局这样内存占用稳定在合理区间。5. 常见问题与排查脚本跑不通的十大坑5.1 下载失败不是代码问题是数据服务问题最常遇到的下载失败有两种一种是响应超时另一种是返回的不是图像而是错误页面。前者用重试机制解决后者要检查拼接的 URL 是否过期。建议在下载循环里加上time.sleep(2)避免高频请求把自己 IP 短时间封掉。还有一种情况是下载过程中断本地生成了不完整的.tif文件。后续 GDAL 打开时会报“不可识别的数据集”。我的经验是写一个文件校验函数在下载后立刻尝试gdal.Open()打不开就自动删掉重新下载。这比事后排查省心太多。5.2 CRS 不一致导致裁剪结果空白裁剪后如果输出栅格是全黑或者面积统计全是 0十有八九是MASK图层和INPUT栅格的坐标系不一致。QGIS 虽然会自动做动态投影转换但某些情况下转换后边界正好落在原始栅格的边缘之外就会产生全空结果。排查方法是裁剪前打印两个图层的crs().authid()确认一致再跑。如果一个是 EPSG:4326一个是 EPSG:32723需要先做重投影。5.3 内存不足大文件读取崩溃处理全国范围的数据时一个裁剪后的栅格文件可能有几十 GB。ReadAsArray()一次性读入内存非常容易崩。可以用ReadAsArray(xoff, yoff, xsize, ysize)分块读取或者先统计每个块的类别再合并结果。def calc_area_by_class_large(raster_path, block_size512): ds gdal.Open(raster_path) band ds.GetRasterBand(1) cols, rows ds.RasterXSize, ds.RasterYSize aggregator {} for x in range(0, cols, block_size): for y in range(0, rows, block_size): block band.ReadAsArray(x, y, block_size, block_size) valid block[block 0] classes, counts np.unique(valid, return_countsTrue) for c, cnt in zip(classes, counts): aggregator[int(c)] aggregator.get(int(c), 0) int(cnt) return aggregator用分块处理后哪怕是 5 GB 的栅格也能在 16 GB 内存的机器上安稳跑完。这是所有脚本里最有价值的小优化之一。5.4 中文字体与乱码问题出图时如果有中文字体需求比如图例标题要写“土地利用类型”Matplotlib 默认字体往往不支持中文会出现方框。解决办法是显式指定中文字体路径。import matplotlib matplotlib.rcParams[font.sans-serif] [Noto Sans CJK SC] matplotlib.rcParams[axes.unicode_minus] False在 Windows 上可以换成[Microsoft YaHei]。CSV 表头乱码则按前面说的用utf-8-sig编码。6. 扩展成真正的自动化工作流6.1 把脚本参数化让同事也能用如果你只是自己用脚本写死研究区路径没问题。但如果想让办公室同事直接拿过去用最好把输入参数抽到脚本顶部或者用配置文件管理。我习惯用一个简单的 JSON 配置{ output_dir: ./analysis/, years: [2010, 2012, 2014, 2016, 2018, 2020], region_shp: ./data/region.shp, tile_index_shp: ./data/mapbiomas_tiles.shp, target_crs: EPSG:4326 }这样换一个研究区只需要改配置文件完全不需要动代码。批量处理多个区域时用 for 循环遍历配置里的 region 列表即可。6.2 无界面运行与定时任务独立脚本初始化 QGIS 环境后可以在不打开 QGIS 的情况下运行。这意味着你可以把整个分析流程加到系统的计划任务里每个月自动跑一次生成最新的土地利用变化报告。Windows 上可以用任务计划程序Linux 上写 crontab。我自己的经验是在 crontab 里跑 PyQGIS 脚本时必须显式设置环境变量比如QGIS_PREFIX_PATH否则有时会找不到 QGIS 插件目录。6.3 再往前走一步的扩展思路这套框架不止能处理 MapBiomas稍微改一改就能套用到其他逐年更新的栅格数据集上。比如连续监测的地表温度产品、植被指数产品或者你自己公司内部的多年分类结果。核心思路不变批量下载、统一坐标系、裁剪统计、时间序列可视化。另外一个可以延伸的方向是空间统计不只是统计面积还可以计算每个类别在研究区内的空间质心变化、边缘密度、破碎化指数。在 PyQGIS 里做这些空间分析也都能找到现成的处理算法完全不用手写底层公式。我在实际项目里用这套脚本帮某机构跑过十几年的土地覆盖变化报告从过去要折腾两周的工作量压缩到半天出结果。最花时间的反而是第一次把代码理顺、把各种下载地址和坐标系验证一遍。等你跑通第一次后续换数据、换区域基本就是复制粘贴改路径的事。这就是做自动化最划算的地方前期投入一小时后面省下几十个小时。