不认识 MapBiomas 的人第一次看到这种栅格数据大概率会手足无措一年一期、30米分辨率、一个像素一个分类代码覆盖整片大陆动辄就是几十个GB的TIFF。手动在 QGIS 里加载、裁剪、查属性表、配色、出图一来二去半天就没了。更别提要做时间序列分析的时候1985 到 2023 年近四十期数据靠手工去点基本就是灾难。这篇文章就是把“用 PyQGIS 自动分析 MapBiomas 并可视化”这整条流水线拆开揉碎从数据下载、投影统一、裁剪掩膜到像素级面积统计、多年批量处理再到自动出图导出一套完整流程。适合被重复性工作折磨的 GIS 从业者、遥感数据分析师也适合准备在毕业论文里快速产出一批专题图的学生。我会把代码、参数和踩坑记录都写在下面。1. 项目拆解与整体思路1.1 我们要解决什么问题MapBiomas 数据的一个核心特点是每个像素的灰度值就是地物分类代码比如 3 代表非森林自然植被、4 代表农业、5 代表牧场、9 代表城市、11 代表水体。它不像普通遥感影像那样需要你做分类而是别人已经把分类结果做好了你直接拿去统计就行。但问题恰恰出在“直接拿去用”这几个字上。原始数据按瓦片发布一景 Tile 覆盖范围很大分辨率 30 米单年文件动辄几个 GB。你手上只有一个流域、一个县或者一个保护区边界直接拿原始栅格跑统计不仅慢还会把大量无关区域的像素算进去。再加上不同年份、不同瓦片之间的投影坐标系可能有细微差异NoData 值处理不当最后统计出来的面积能偏差得让你怀疑人生。这篇文章要解决的就是这套脏活批量加载、统一投影、按研究区裁剪、分年统计各类别面积、生成变化曲线、自动套配色出图。目标是让这些原本要在 QGIS 界面里点一小时的操作变成一个脚本跑三分钟的事情。1.2 方案设计与技术选型为什么用 PyQGIS 而不是直接用 GDAL 或者纯 NumPy我的判断是PyQGIS 站在 QGIS 这个生态上既能拿到算法库的处理能力又能直接复用 QGIS 的渲染、布局和导出能力。如果你用纯 GDAL读完波段之后要自己写重投影、裁剪、渲染配色、出图这套东西代码量会翻好几倍。而 PyQGIS 里一个processing.run()就能调用 QGIS 自带的裁剪、重分类等算法样式和布局这些和我们在界面上操作完全对应学习成本低很多。整条流水线我把它拆成四段数据准备确认研究区范围、下载对应瓦片、检查坐标系和分辨率。预处理统一投影、按矢量边界裁剪、处理 NoData。统计与建模像素级分类统计、多年面积时间序列、变化检测。可视化与输出设置分类配色、自动创建布局、批量导出专题图。这个顺序千万不要乱。很多人一上来就写统计代码结果投影没统一面积算得有问题返工成本极高。规格化流程的意义在于每一步的输出都被下一步消费前面错了后面白做。1.3 环境准备PyQGIS 怎么才能跑起来QGIS 自带的 Python 环境里已经预装了 PyQGIS 和相关依赖理论上打开内置终端Python Console就能直接开始写。但我建议你别在控制台里写长脚本而是在本地用文件编写再丢到 QGIS 里执行。我的环境版本组合供参考QGIS 3.28 LTRPython 3.9操作系统是 64 位 Linux。如果你用 Windows路径写法注意区分反斜杠和正斜杠尤其是栅格路径里有空格的时候建议一律用正斜杠。验证环境是否可用第一段代码就这么写from qgis.core import QgsApplication, QgsProject, QgsRasterLayer print(QGIS 版本:, QgsApplication.version()) print(PyQGIS 导入成功)如果你在 QGIS 内置终端跑QgsApplication已经被初始化了直接能用。如果你要在独立的 Python 脚本里跑需要手动初始化。建议刚开始不要折腾独立脚本就先在 QGIS 内置终端里写等整条链路跑通了再考虑脱离界面批量执行。提示处理大栅格的时候QGIS 内置终端是单线程的脚本里如果跑了很重的算法界面会暂时卡住这是正常的。想改善体验可以在循环里加QgsApplication.processEvents()让界面有时间刷新。2. 数据获取与预处理2.1 MapBiomas 数据下载与图层选择MapBiomas 的数据集有若干 Collection 版本年份从 1985 年到最近年份逐年发布每个年份都是一期独立的栅格数据。下载时你需要选三样东西Collection 版本、目标年份、瓦片编号。瓦片选择这一步很多人想偷懒直接把地图范围里的瓦片全下下来。我试过一个 6 行 4 列的瓦片范围下载量轻松超过 40GB后面处理起来又慢又占硬盘。正确做法是先用矢量边界生成最小外接矩形然后人工确认覆盖这个矩形的瓦片编号只下载这些瓦片。以某研究区为例跨两个瓦片我就只下载对应编号的两三个文件而不是把整个数据目录端回来。这一步能省掉 80% 的磁盘空间和后续 IO 时间。2.2 影像重投影与范围统一MapBiomas 原始数据的投影一般是地理坐标系例如 EPSG:4326单位是度。问题在于如果用经纬度的栅格直接算面积那每个像素的真实面积会随纬度变化根本不能用“像素数 × 固定面积”来统计。我的做法是统一重投影到Albers 等积投影。名字里带“等积”两个字意思是它能在较大区域内保证面积不变形对于这种大陆尺度的土地覆盖统计是最合适的。实际操作中你可以根据研究区纬度选择合适的中央经线和标准纬线如果嫌参数麻烦直接用数据集官方推荐的投影参数也行。重投影这一步推荐用 QGIS 内置的“重投影”算法在 PyQGIS 里通过processing.run()调用import processing params { INPUT: raster_path, TARGET_CRS: EPSG:XXXX, # 替换为你的等积投影 RESAMPLING: 1, # 双线性分类栅格推荐最近邻是0 OUTPUT: projected_path } processing.run(gdal:warpreproject, params)注意一个细节分类栅格在重投影时重采样方法务必用最近邻NEAREST不要用双线性或三次卷积。因为分类代码是类别编号不是连续数值插值会产生 2.7 这种根本不存在的类别污染后续统计。这个坑我第一次跑就踩了处理完的数据里冒出来一堆“新类别”当时我还以为是数据本身有问题。2.3 裁剪与研究区掩膜投影统一后下一步是把大栅格裁到研究区范围内。这里用矢量边界对栅格进行掩膜裁剪params { INPUT: projected_path, MASK: boundary_vector_path, NODATA: 0, ALPHA_BAND: False, CROP_TO_CUTLINE: True, KEEP_EXTRAA_BANDS: False, OUTPUT: clipped_path } processing.run(qgis:cliprasterbymasklayer, params)CROP_TO_CUTLINETrue的意思是不仅掩膜还把范围裁掉这样输出栅格只有研究区外接矩形那么大后续numpy读取的时候内存压力小很多。还有NODATA参数。裁剪后研究区之外的部分会变成 NoData之内的像素保留原分类值。我建议把 NoData 显式设为 0 或者 255这样在统计的时候好过滤。注意掩膜裁剪本身不会改变分类值但它会在边界处产生非常细碎的像素。如果你的矢量边界和栅格边界不完全对齐边界会出现半个像素。对 30 米分辨率的 MapBiomas 来说单个像素面积是 900 平方米边界区域几十个像素的偏差在统计上影响不大可以忽略。3. PyQGIS 自动化统计从栅格到数据表3.1 核心统计思路拆解MapBiomas 的分类栅格统计面积的核心逻辑只有一句话统计每个分类代码出现的像素个数再乘以单个像素的真实面积。单个像素面积的计算有两种方式。如果你用的是等积投影像素尺寸在投影单位下基本不变那单像素面积就是分辨率乘分辨率比如 30 米分辨率就是 30 × 30 900 平方米。如果你用的还是经纬度的栅格那就得逐像素计算纬度方向上的面积麻烦得多。这也再次说明投影统一是前置条件。整个统计流程的伪代码是这样的for each 年份: 加载裁剪后的栅格 读取唯一值和像素计数 过滤掉 NoData 值 计算每个类别的面积 像素数 × 单体面积 存入数据表有了逐年数据表你就能画时间序列、计算增长量和变化率最后再把这些数字和图表配套进专题图里。3.2 基础代码框架加载图层与常用对象PyQGIS 处理栅格绕不开三个类QgsRasterLayer表示栅格图层QgsRasterDataProvider用来实际读取波段数据QgsProject用来管理图层并映射到 QGIS 界面。加载一个裁剪后的栅格到 QGIS 项目里代码是这样的from qgis.core import QgsProject, QgsRasterLayer raster_path /path/to/clipped_2020.tif layer QgsRasterLayer(raster_path, MapBiomas 2020) if not layer.isValid(): raise Exception(图层加载失败请检查路径和文件格式) QgsProject.instance().addMapLayer(layer)isValid()这个检查一定要做。栅格文件格式损坏、路径写错、空间参考缺失都会导致图层无效而无效图层上的任何操作都会在后边以奇怪的方式报错。我通常先加载一个图层进项目看看再用脚本处理基本上能避免一半的“为什么结果为空”问题。3.3 重分类把代码映射成可读类别MapBiomas 的分类代码是一堆数字3、4、5、9、11……数字本身没有可读性。在出图前我们要把数字映射成中文类别名和配色信息。我维护一张映射表把它放在代码最前面MAPBIOMAS_CLASSES { 1: {name: 森林, color: #1f7003}, 3: {name: 非森林自然植被, color: #27a02c}, 4: {name: 农业, color: #d6bc74}, 5: {name: 牧场, color: #ecec8e}, 9: {name: 城市, color: #d4271e}, 11: {name: 水体, color: #2c55a1}, 12: {name: 其它非植被, color: #af2a2a}, 15: {name: 牧场农业混合, color: #b8af4e}, }为什么要单独维护一张映射表而不是在统计代码里写死因为后续的可视化和面积统计都要引用这里的信息。如果你把类别名散落在各个函数里改一次类别对应关系要改三个地方很容易改漏。一张表统一管理统计和配色都从这里取后续扩展类别时只需改一处。如果研究区的分类集合不是全类别你需要先看一眼栅格的实际值分布再定映射表可以用下面的方法快速侦察import numpy as np def unique_values(path): layer QgsRasterLayer(path, tmp) provider layer.dataProvider() extent layer.extent() block provider.block(1, extent, layer.width(), layer.height()) data block.numpy() values, counts np.unique(data[~data.mask], return_countsTrue) return dict(zip(values.tolist(), counts.tolist()))这里我把 NoData 用 numpy 的 masked array 处理了。provider.block()返回的就是一个np.ma.MaskedArrayNoData 位置自动标记为 mask。统计唯一值的时候一定用data[~data.mask]把 NoData 剔除掉。3.4 面积统计像素级计算核心统计我不用 QGIS 的统计面板直接用 NumPy 跑速度会快得多。一个 30 米分辨率的裁剪栅格像素量一般在千万级别NumPy 处理也就是几秒的事。def compute_area_by_class(raster_path, pixel_size900.0): layer QgsRasterLayer(raster_path, tmp) provider layer.dataProvider() extent layer.extent() block provider.block(1, extent, layer.width(), layer.height()) data block.numpy() valid data[~data.mask] values, counts np.unique(valid, return_countsTrue) area_dict {} for code, cnt in zip(values, counts): if code in MAPBIOMAS_CLASSES: area_dict[code] cnt * pixel_size return area_dict这里有个我经常犯迷糊的点provider.block()的第四、五个参数是“宽高”而不是“在栅格中的行列范围”。第一次写的时候我把extent和width, height传错了读出来的数组维度完全对不上。所以每次写这行代码我都提醒自己看一眼layer.width()和layer.height()。再提醒一个细节pixel_size这个变量默认是 900也就是 30×30但前提是你已经做了等积投影。如果你用的是别的投影这个值必须重算否则面积就失真了。如果你觉得 NumPy 太底层也可以用 QGIS 的 “rasterlayeruniquevaluesreport” 算法它能直接输出唯一值统计表包括计数。但对于批量处理来说NumPy 不需要生成中间文件内存里直接算完效率更高。两者我都试过批量处理一千期小图斑时NumPy 方案大概快五倍。3.5 多年批量处理从手点到一键前面的统计函数只处理单年现在用循环把它变成时间序列。MapBiomas 从 1985 年开始到最近年份结束逐年遍历import csv from os.path import join years range(1985, 2023) results [] for year in years: raster_file fmapbiomas_clip_{year}.tif raster_path join(clipped_dir, raster_file) if not os.path.exists(raster_path): print(f跳过缺失年份: {year}) continue areas compute_area_by_class(raster_path) row {year: year} for code, info in MAPBIOMAS_CLASSES.items(): row[info[name]] round(areas.get(code, 0) / 10000, 2) # 平方米转公顷 results.append(row)为什么要转成公顷MapBiomas 的面积尺度太大用平方米动辄就是几亿设置图例的时候数位过长不方便阅读。公顷是土地覆盖统计里比较常用的单位比平方公里精细一些比平方米友好得多。输出 CSV 也很直接with open(class_area_by_year.csv, w, newline, encodingutf-8) as f: writer csv.DictWriter(f, fieldnames[year] [c[name] for c in MAPBIOMAS_CLASSES.values()]) writer.writeheader() writer.writerows(results)有一年我因为没有加encodingutf-8用 Excel 打开 CSV 中文全是乱码后来在csv模块里指定编码才解决。Windows 上的道友尤其注意utf-8加上Excel 可能默认按 GBK 读取但至少保存的数据本身是干净的。3.6 变化检测快速找出面积变化最大的类别时间序列统计完成后变化检测就是纯表格操作了。我拿“森林”类别举例。假设上面的results已经是一个列表的 dict我要提取每年的森林面积并计算“相对于上一年”的变化量import pandas as pd df pd.DataFrame(results) df[forest_change] df[森林].diff() df[forest_change_rate] df[forest_change] / df[森林].shift(1) * 100这样一下就能看到哪一年的森林面积骤减哪一年有明显恢复。如果研究区历史上发生过大规模事件比如干旱、火灾或农业扩张这个表就会在对应年份出现明显的异常点也就是变化检测的入口。再进一步如果你想定位“空间上哪里有变化”那就需要做两期栅格的逐像素比较。用 NumPy 读取两个年份的数组按类别布尔比较forest_2020 (data2020 1) forest_2019 (data2019 1) converted np.logical_and(~forest_2020, forest_2019) # 从森林变成其他得到的布尔矩阵可以统计像素数、换算面积也可以进一步转为矢量多边形。注意转矢量时要用gdal:polygonize算法并且先对布尔矩阵做一个中值滤波去掉单像素噪声不然会产生海量碎片多边形后续根本没法看。4. 可视化出图自动生成专题地图4.1 分类配色设计统计完之后大多数人的需求并不是一张 CSV而是一张拿得出手的专题图。MapBiomas 官方有一个标准的分类配色方案表现为不同类别用不同颜色区分森林是绿色系农业是浅黄/棕色系水体是蓝色系城市是红色系。你完全可以照着这个逻辑来自定义配色。我在代码里维护的MAPBIOMAS_CLASSES字典中就已经包含color键了就是为了渲染时直接取用。配色有两点建议类别之间色差要足够大至少能让人一眼分辨出林地和农田。水体和城市这两个类别虽然面积可能不大但要用高对比颜色蓝色和红色因为人工检查时往往最先注意这两类。4.2 图层样式自动设置QGIS 里设置栅格分类样式对应的是QgsPalettedRasterRenderer。它接收一组 (值, 颜色, 标签) 作为分类条目。from qgis.core import QgsPalettedRasterRenderer, QgsColorRampShader, QgsRasterShader def apply_mapbiomas_style(layer): colors [] for code, info in MAPBIOMAS_CLASSES.items(): colors.append(QgsPalettedRasterRenderer.Class(code, QColor(info[color]), info[name])) renderer QgsPalettedRasterRenderer(layer.dataProvider(), 1, colors) layer.setRenderer(renderer) layer.triggerRepaint()这个代码把每个分类值对应到一种颜色并显示类别名。如果不设置你会发现栅格加载出来全是灰的或者假彩色因为 QGIS 默认的栅格渲染器不会自动识别分类数据的语义。还有一个坑有时候图层设置了样式但出图时却看不到颜色原因是 QGIS 的布局项默认跟随图层的当前渲染器而设置渲染器后没有调用triggerRepaint()。我每次设置完样式都会主动刷新一下确保界面和布局都拉到最新状态。4.3 布局与出图样式设置好后出图环节是个体力活。QGIS 布局脚本化有点绕但掌握套路后就很顺手。核心逻辑分四步创建布局、添加地图项、添加图例和比例尺、导出图片。我先给一个简化的出图代码from qgis.core import ( QgsPrintLayout, QgsLayoutItemMap, QgsLayoutItemLegend, QgsLayoutItemScaleBar, QgsLayoutPoint, QgsLayoutSize, QgsUnitTypes ) project QgsProject.instance() def generate_layout(layer, output_path): layout QgsPrintLayout(project) layout.initializeDefaults() layout.setName(专题图) # 地图项 map_item QgsLayoutItemMap(layout) map_item.setRect(20, 20, 260, 170) # x, y, w, h单位 mm map_item.setLayers([layer]) layout.addLayoutItem(map_item) map_item.attemptResize(QgsLayoutSize(260, 170)) map_item.setExtent(layer.extent()) # 图例 legend QgsLayoutItemLegend(layout) legend.setLinkedMap(map_item) layout.addLayoutItem(legend) legend.attemptResize(QgsLayoutSize(60, 90)) # 比例尺 scalebar QgsLayoutItemScaleBar(layout) scalebar.setLinkedMap(map_item) scalebar.setUnits(QgsUnitTypes.DistanceMeters) scalebar.setNumberOfSegments(3) scalebar.update() layout.addLayoutItem(scalebar) scalebar.attemptMove(QgsLayoutPoint(20, 210)) # 导出 exporter QgsLayoutExporter(layout) exporter.exportToImage(output_path, QgsLayoutExporter.ImageExportSettings())这套流程里最需要调的是坐标和尺寸。attemptResize和attemptMove的单位不是像素而是布局页面的单位通常是毫米我刚上手的时候把长宽写成像素结果图在整个布局里只剩一个小角落。图例摆放也要注意如果你加了多个图层图例会默认把项目里所有图层都列出来。我只想显示当前分类栅格的图例就调用legend.setLinkedMap(map_item)并手动移除多余图层。这一步可以在布局初始化循环里过滤一遍for item in layout.items(): if isinstance(item, QgsLayoutItemLegend): root item.modelRootGroup() for child in root.children(): if child.name() ! layer.name(): root.removeChildNode(child)4.4 组合一键输出多年图集当你有 1985 到 2023 年的数据想做一套逐年专题图时手工出图的成本完全不可接受。脚本化就是解决最后一公里问题的核心。我的做法是用for year in years循环每一期生成自己的图层路径调用apply_mapbiomas_style和generate_layout输出文件名带年份for year in years: layer_path join(clipped_dir, fmapbiomas_clip_{year}.tif) layer QgsRasterLayer(layer_path, fMapBiomas {year}) if not layer.isValid(): continue QgsProject.instance().addMapLayer(layer) apply_mapbiomas_style(layer) output_path join(output_dir, fmapbiomas_{year}.png) generate_layout(layer, output_path)这样跑完几十张图全部生成而且配色统一、范围一致、图例格式统一。你要是做 PPT 汇报或者写报告这套图集直接能用。也可以生成 GIF 动态图来展示时间序列变化把所有 PNG 用 pillow 合成帧动画会非常直观。注意导出 PNG 的分辨率要保持一致不然合成 GIF 时会拉伸。5. 实战中容易踩的坑与排查记录5.1 CRS 不统一导致面积错乱症状统计出来的面积和官方报告对不上甚至不同年份之间同一地区的森林面积差了好几倍。原因某个年份的栅格重投影参数写错了或者下载到的瓦片本身 CRS 不是预期的地理坐标被强行按等积投影读取。排查流程用QgsRasterLayer.crs().authid()打印每个输入图层的 CRS和TARGET_CRS对比。如果发现某个年份的输出 CRS 不对就回看gdal:warpreproject的日志多半是重投影目标参数被覆盖了。我在一次批量处理中因为循环里复用了同一个params字典前一年的TARGET_CRS被后一年意外修改导致后几年全部投影错误。解决方案是每个循环迭代里重新创建params字典不要复用同一个对象。5.2 NoData 像素污染统计结果症状分类统计里多出一个极小的类别比如“255”或者“0”并且面积忽大忽小。原因裁剪后的栅格周围一圈是 NoData但 NumPy 读取时如果 mask 没正确识别NoData 值会变成真实像素参与统计。更常见的是初始下载的原始栅格本身就有边界 NoData重投影时没有正确赋值。解决办法裁剪时显式设置NODATA并在 NumPy 统计里用 masked array 过滤。在compute_area_by_class里我特意用了data[~data.mask]这一行就是防这个的。经验如果数据里出现未知代码不要直接忽略先打印出来看看数值。有些 MapBiomas 版本会把“烧毁区域”标记为特定代码如果研究区发生过火灾这些代码是真有含义的。忽略它们会低估某些类别面积。5.3 QGIS 卡死与内存溢出症状处理某一期大栅格时QGIS 界面卡死甚至直接崩溃。原因大栅格文件加载进 QGIS 时也会占用大量内存。如果循环里每个年份的图层都被addMapLayer加进去了项目会越来越臃肿几十期图层全堆在 Layers 面板里内存很容易爆。解决策略在循环里处理完一个图层后用QgsProject.instance().removeMapLayer(layer)把它移出项目或者根本不要addMapLayer只做统计和渲染后直接删掉。如果只是想出图用QgsRasterLayer但不要添加到项目样式设置和布局引用都是可以的。在重循环里适当调用QgsApplication.processEvents()让 QT 事件循环喘息。我后来稳定的做法是统计阶段不往项目里加任何图层全部用独立变量操作只有最终要出图的图层才加入项目并在出完图后一并清理。5.4 图例颜色和地图对不上症状地图配色是正确的但图例里某一类显示的颜色和地图里的颜色不一样。原因我遇到过两种。一是自定义了 QgsPalettedRasterRenderer 后没有清理旧的渲染器缓存二是图例关联到了错误的图层比如关联到了另一个没有套用样式的临时图层。排查步骤先确认地图项里的图层和渲染器一致再确认图例的 model 中引用的图层名匹配。大部分情况都是“图例没有关联到地图项”造成的。5.5 脚本运行报错 QuickFixModuleNotFoundError: numpyQGIS 自带 Python 环境一般带 numpy。如果是在独立终端里跑可能需要pip install numpy。QgsApplication not initialized独立脚本里忘了建 QgsApplication 实例。参考官方文档按顺序QgsApplication([], False)→initQgis()。Algorithm not found: qgis:xxx算法 ID 拼写错误或版本不同。在代码里输出processing.algorithmHelp(算法ID)查看参数列表。汇总成一张速查表问题现象可能原因解决动作面积偏大或偏小CRS 没统一或像素尺寸用错layer.crs()打印确认重投影后再统计面积统计出现奇怪类别NoData 参与统计numpy 里用~data.mask过滤界面卡死项目图层过多内存不足统计时不加图层只处理不显示图例与地图颜色不一致图例关联错图层/缓存未刷新重新关联地图项triggerRepaint()中文输出乱码CSV 编码问题encodingutf-8-sig可兼容 Excel我的实操体会与扩展建议整套流水线跑通之后我最大的感受是这类“半自动”脚本对日常分析帮助远大于重写一套完整的处理平台。你不要指望第一次写出来的代码就完美先从单年份手动跑一次确认统计结果和 QGIS 界面手动操作结果一致再把循环套进去。如果直接上手循环出了问题你根本不知道是数据问题、参数问题还是裁剪问题。还有一个经验把中间产物重投影后栅格、裁剪后栅格都保留下来别贪图省空间直接覆盖原文件。因为没准你后面想换个研究区边界重新裁剪或者想对比不同 Collection 版本到那时候原始数据和中间图层都还在你只需要重新写一行路径省得全部重新下载。后续如果想再进一步可以把这个脚本改造成自动输出报告用 matplotlib 画面积时间序列折线图和专题图一起拼成一个 PDF汇报的时候直接拿结果就行。更自动化一点也可以加入文件监控数据目录里新下载了一期年度数据就自动触发处理把“数据分析”变成“数据流水线”。这些都是在现有这套框架上一个模块一个模块加出来的不需要推倒重来。对于处理 MapBiomas 这类时序栅格数据核心思路就是一个字“批”。把所有重复操作都交给 PyQGIS 脚本你只需要把注意力和时间留给那些真正需要人类判断的部分比如变化原因的解读、边界区域的核查。这才是工具应该有的位置。