做绿地格局分析的人大概率都遇到过这种场面同一份土地利用分类图两个人各自算绿地斑块密度结果一个0.8一个3.6谁都觉得自己没错。绿地斑块密度Patch DensityPD看起来就是一个数量除以一个面积但真正从图斑统计走到指数测算里面隔着坐标系、最小制图单元、统计口径、单位换算四道坎。这个指标又是城市绿地系统评估、生态破碎化分析、遥感解译后处理里最常被引用的基础指数之一算错了后面所有结论都跟着歪。这篇文章就把“绿地斑块密度计算”整条链路拆开讲先敲定定义和口径再处理原始数据得到干净的图斑接着走一遍手算过程最后给一个可以直接套用的Python自动化脚本并用一组虚拟数据复盘最容易翻车的三个地方。不管是做课程设计、毕业论文、规划项目还是科研数据处理照着这个流程走踩坑概率会小很多。1. 把“绿地斑块密度”的定义钉死口径不一致时一切对比都是空谈1.1 PD公式很简单但简单不等于不容易错斑块密度最常见的表述是单位景观面积内斑块的数量。公式写成PD N / A其中N是目标类别的斑块总个数A是景观总面积。常用单位是“个/平方千米”少数文献会用“个/公顷”或“个/100公顷”。换算关系很简单1平方千米100公顷所以“个/平方千米”和“个/100公顷”在数值上完全一致。这一点很多人不知道等会儿和Fragstats输出对照时就能用上。为什么这个公式容易错因为它把“斑块”和“面积”都当成了现成参数实际上“斑块N”取决于你怎么切分、怎么去碎斑“面积A”取决于你用什么分母。N和A只要有一个口径漂移PD就不是同一个量。我见过一份报告里写“绿地斑块密度2.5个/ha”另一份写“25.0个/km²”看起来数值差别大实际上这俩是同一个数。单位没写清楚才是真问题。1.2 分母用景观总面积还是绿地总面积差出数倍PD的分母可以是景观总面积也可以是该类别斑块总面积。前者衡量斑块在研究区内的绝对分布密度适合做不同区域间的横向对比后者衡量类别自身被切割成多少份更接近“破碎化程度”。绝大多数景观格局分析默认用前者但很多初学者会把分母写成绿地面积。用一个虚拟数字示例。某研究区景观总面积12.6平方千米绿地斑块178个绿地总面积3.62平方千米用景观总面积做分母PD 178 / 12.6 ≈ 14.1个/km²用绿地总面积做分母PD 178 / 3.62 ≈ 49.2个/km²同一个数据集分母一换结果差4.5倍左右。所以论文或报告的方法部分一定要写明“分母为景观总面积”。否则数据一进入横向比较整个结论都会失真。1.3 最小制图单元MMU决定斑块数量也决定密度量级比分母更隐蔽的是“斑块数量N”的取值。同一块区域用10米分辨率影像提取出的碎斑数量会远大于30米分辨率影像单个像元噪点也会形成独立图斑。这个问题的解决办法是设定一个最小制图单元Minimum Mapping Unit, MMU小于该面积的图斑不参与统计。MMU不是一个固定值取决于数据源和应用场景。生态类研究常用1000平方米量级城市街道尺度可能用几十到几百平方米。这个选择应该被写进方法描述里因为不同MMU会让“N”发生几十个百分点以上的变化PD自然跟着变。算密度前先想清楚“什么算一个斑块”比急着点计算器更重要。2. 从原始分类数据到干净图斑栅格转面、投影和碎斑处理2.1 第一步先建一张“绿地类别对照表”不同项目对“绿地”的定义差别很大。有的把林地、草地、灌丛都算绿地有的把城市公园绿地、防护绿地、附属绿地纳入有的还会把水体单独处理把农田单独拎出去。如果只是机械地把分类图里几个代码选出来后面所有统计都失去了可比性。我建议数据处理前先做一张对照表包含类别名称、原分类代码、是否属于本研究的绿地、备注。例如类别名称分类代码是否计入绿地备注乔木林地21是包含片林、行道树带灌木林地22是连续灌丛草地23是城市草坪、天然草地农田31否视研究目标可单列水体41否通常不计入绿地建设用地51否不透水面这张表建议写进项目文档。因为不同年份的数据如果分类代码略微调整有对照表才能快速同步口径也是后续自动化脚本里最需要维护的配置文件。2.2 栅格分类图转矢量图斑简化面选项要慎用如果手头是栅格分类图第一步通常是栅格转面。在ArcGIS里执行Raster to Polygon时有一个“简化面”相关的选项我的建议是保持关闭或者使用保留原始像元边界的模式。原因是简化面会把锯齿状边界抹平抹平过程中可能改变小斑块的形状甚至面积。转完矢量之后属性表里会多一个VALUE字段存放原始栅格值。这个字段就是类别代码之后所有筛选和分组都以它为主键。如果原始数据已经是矢量则需要先跑一遍“检查几何”和“修复几何”。面状数据里常见的自相交、重叠、缝隙在后续面积统计时会造成总和不准确早清理早消停。2.3 投影坐标系是算面积的前提经纬度坐标不能直接算这是个老生常谈却依然频繁出现的坑。WGS84经纬度坐标下直接算面积的操作结果只是“度²”量纲完全没有物理意义。正确的做法是把数据投影到等面积投影坐标系再计算面积。常用选择小范围用UTM对应分带或高斯-克吕格投影跨多个分带、做全国尺度的分析用Albers等面积投影。判断方法简单粗暴打开图层属性查看坐标系如果单位是“Degree”就一定要做Project如果单位是“Meter”才能继续算面积。一个特别容易忽略的细节是研究区边界和斑块图层必须使用同一套投影。否则用边界面积做分母、用斑块面积做分子单位都不能严格对齐算出来的PD自然不靠谱。2.4 小图斑的清理方式直接写入处理记录定义好MMU之后用属性查询把所有小于阈值的斑块选中并删除操作过程很简单。比如ArcGIS属性表里执行Shape_Area 1000 AND VALUE IN (21, 22, 23)如果绿地类别需要单列就先按类别筛出绿地再做面积过滤如果做全部地类的景观水平密度那么就先保留所有类别只按面积过滤。这里涉及一个研究层面问题删除碎斑和把碎斑合并进相邻大斑哪种处理更合理我个人的建议是如果只是计算密度和面积比例直接删除可复现性最强如果后续要做形态分析、连通性分析可能需要做融合处理。无论选哪种都要把处理方式写进方法文本否则别人无法对你的N进行审计。3. 图斑属性表里的统计博弈把数量、面积、边界长度规整到一张表3.1 面积字段怎么加不同软件的常规操作矢量斑块拿到手后首先给每个图斑添加面积字段。ArcGIS里可以用添加几何属性工具输出单位选“平方米”QGIS里用字段计算器表达式$areaPython里使用geopandas时每个要素的geometry.area会自动按当前投影坐标单位给出面积。值得提示的是如果数据来自不同机构有的面积字段可能叫Area_km2有的叫Area_ha还有的干脆没有。统一命名和单位在这一步一定要做我习惯全部统一成平方米字段名统一为area_m2。字段名一致后脚本可以批量处理多个文件。3.2 分组聚合统计斑块数量和类别总面积统计每个类别的斑块数用属性表的分组统计或数据库查询都能实现。ArcGIS里可以把属性表导出后用Excel透视表也可以直接用“汇总统计数据”工具分组。本质上都是对“类别代码”做分组聚合斑块个数N用COUNT统计类别总面积用SUM统计聚合结果是一张每个类别一行的表。关键来了这张表里的面积之和不一定是研究区总面积。因为研究区里可能有未分类区域、可能有编辑时的缝隙也可能因为数据叠加问题出现重复。计算PD时分母必须单独从研究区边界矢量上取这一条要刻在脑子里。3.3 单位换算一次推导永远不再搞混假设研究区总面积字段study_area_m2、某类斑块个数N、某类总面积area_m2那么斑块密度PD N / (study_area_m2 / 1000000)单位个/km²绿地覆盖比例PLAND area_m2 / study_area_m2 × 100单位%平均斑块面积MPS area_m2 / N单位m²前面说过1平方千米等于100公顷所以“个/km²”和“个/100公顷”数值相等。Fragstats输出的PD单位是每100公顷的斑块数读取结果时直接按“个/km²”理解即可不需要再做任何换算。这个小知识点能省下不少和软件输出结果大眼瞪小眼的时光。3.4 手工计算和Excel透视表的典型流程如果只算一次手动操作完全够用。典型步骤导出属性表到Excel→用透视表统计类别对应的行数和面积总和→把研究区面积作为常量放入单元格→用公式算出PD和PLAND→顺手加一列注释记录MMU、投影、数据版本。这个方法的好处是每一步都可见方便核对。缺点是一旦有几十个区域、几十个年份的数据透视表就要做几十次公式复制时容易漏引用。此时就该进入下一章的自动化方案。4. 自动化一步到位用Python把整条链路封装起来4.1 为什么要自动化口径固定比省时间更重要手工算一次只要十几分钟如果数据量小谈不上瓶颈。但景观格局分析通常不是一次性的可能要对同一城市不同年份、同一片区不同方案反复测算。手工流程每重复一次就有一次漏掉碎斑过滤、错选字段的风险而且多人协作时各种操作习惯都会悄悄混进来。把流程写成脚本等于把“怎么算”这件事实物化、版本化。代码就是方法记录参数就是口径清单。以后数据更新重新跑一遍就能得到新结果还能确保和上一版结果严格同口径。我认为这才是“一步到位”的真实含义不是省略步骤而是把步骤固定下来。4.2 一个可直接套用的Python计算脚本以下脚本使用geopandas输入是已经转换到投影坐标系的矢量面文件比如landuse.shp。绿地类别代码通过GREEN_TYPES列表控制MMU通过MMU_M2控制研究区总面积通过STUDY_AREA_M2控制。import geopandas as gpd # ---- 参数区开始 ---- GREEN_TYPES [21, 22, 23] # 绿地类别代码按自己的对照表配置 MMU_M2 1000.0 # 最小制图单元单位平方米 STUDY_AREA_M2 12_600_000.0 # 研究区总面积单位平方米从研究区边界矢量计算而来 INPUT_FILE landuse.shp # 投影坐标系矢量面文件 EPSG_PROJ EPSG:32650 # 若原始数据无投影或为经纬度先做重投影 # ---- 参数区结束 ---- gdf gpd.read_file(INPUT_FILE) # 确保使用投影坐标系本文以UTM 50N为例 if gdf.crs and gdf.crs.is_geographic: gdf gdf.to_crs(EPSG_PROJ) # 添加面积字段并过滤小斑块 gdf[area_m2] gdf.geometry.area gdf gdf[(gdf[landuse_code].isin(GREEN_TYPES)) (gdf[area_m2] MMU_M2)].copy() # 按类别聚合斑块个数、面积总和 rho gdf.groupby(landuse_code)[area_m2].agg([count, sum]) rho rho.rename(columns{count: patch_count, sum: total_area_m2}) # 计算景观水平的常见密度指标 rho[pd_per_km2] rho[patch_count] / (STUDY_AREA_M2 / 1_000_000.0) rho[pland_percent] rho[total_area_m2] / STUDY_AREA_M2 * 100.0 rho[mean_patch_area_m2] rho[total_area_m2] / rho[patch_count] rho.reset_index(inplaceTrue) print(rho) # 如果需要总绿地口径 total_count rho[patch_count].sum() total_green_area rho[total_area_m2].sum() print(f总体斑块数{total_count}个绿地面积{total_green_area / 1_000_000.0:.2f}平方千米 fPD {total_count / (STUDY_AREA_M2 / 1_000_000.0):.2f}个/平方千米)实际使用时landuse_code字段名要和自己的数据一致绿地代码也要按前文对照表调整。脚本会自动滤掉小于MMU的碎斑并根据研究区总面积计算每个绿地类别的PD、绿地覆盖度、平均斑块面积最后输出总体数值。值得注意脚本把分母统一为STUDY_AREA_M2这就是口径固定最好的体现。4.3 不写代码也可以串成“一步到位”的模型不会Python的分析师可以用ArcGIS模型构建器或QGIS图形化建模器实现相同功能。节点顺序大概是栅格转面如果需要按属性筛选绿地图斑并过滤小面积斑块添加几何属性计算面积汇总统计分组求斑块数和面积添加研究区面积常量字段计算器计算PD、PLAND、MPS模型建好后每次只要替换输入图层和参数就能一键出结果。我仍然推荐同时维护一份Python脚本因为脚本更容易做批量循环和版本管理。两种方式并不冲突都是把重复劳动收敛成一个单元。4.4 用Fragstats交叉验证一次能救回不少低级错误脚本算完建议用Fragstats做一次抽查验证。操作要点输入原始栅格分类图把背景值Background设置为0或-9999以外的实际不参与分析的类别千万不要把非绿地区域全部设为NoData否则NoData不进入景观总面积分母会缩水到只剩下绿地面积PD数值直接膨胀数倍。Fragstats输出的PD如果和脚本结果对不上优先检查三处背景值是否一致、最小图斑是否在栅格预处理阶段被合并、面积单位是否正确。我实操中遇到的无数次“对不上”最后都能归到这三件事上数据处理环节的重合率非常高。5. 用一组虚拟数据把数字跑通复盘最容易翻车的三个地方5.1 两张对照表实战先看结果再解释意义我用一组模拟数据来演示输出形式。假设A城市高新区研究区面积12.6平方千米B城市新城区研究区面积18.4平方千米绿地类别都是林地、灌木、草地MMU设为1000平方米。处理后得到区域绿地斑块数绿地面积(km²)PD(个/km²)绿地覆盖度(%)平均斑块面积(m²)A城市高新区1783.6214.128.720337B城市新城区2457.3613.340.030041只看PD两个区域似乎接近一个14.1、一个13.3。但结合绿地覆盖度就发现B区绿地面积占比明显更高但平均斑块面积也更大。这说明A区的绿地破碎度相对更高B区有更多成片大绿地。如果只报一个PD值很容易把差距讲歪。密度类指标的价值在于配合覆盖度和平均斑块面积一起解读而不是单独看某个数字。5.2 翻车点一碎斑阈值一改数量直接翻倍同样用A城市高新区的模拟数据如果把MMU从1000平方米降为0也就是不过滤碎斑斑块数会从178飙到406PD从14.1变成32.2个/km²。这个变化完全不是因为区域发生了变化只是处理口径变了。这意味着什么任何跨时间、跨区域的PD比较都必须先锁定MMU。否则2020年的分析用了1000平方米2025年的分析用了100平方米哪怕真实现状没有太大改变PD也会呈现“假性倍增”。做这类指标分析时把MMU写进表头或者研究方法里是对自己后续结论的一种保护。5.3 翻车点二没有统一投影面积偏差悄悄混进结果前面提到经纬度不能直接算面积。实际操作里还有一种更隐蔽的口径漂移部分数据是地理坐标系部分数据是投影坐标系汇总时没做统一转换就用软件默认方式叠在一起。结果可能是研究区面积被系统自动换算斑块面积却没同步最终PD偏差几个百分点还很难察觉。批量处理多个文件时建议在脚本里预先判断坐标参考系若为地理坐标系就强制重投影。前面Python代码里已经演示了如何检查。这句gdf.crs.is_geographic实际上是一个很简单的护栏能拦下大量面积统计事故。5.4 翻车点三分母从这个变量悄悄换成了那个变量最后一个坑来自复制代码时的粗心。假设某次计算把“研究区总面积”写成了“绿地总面积”178个斑块除以3.62平方千米得到49.2个/km²和正确值14.1完全不在一个量级。这种错不会让程序报错输出看着也很合理恰恰最难发现。规避方法很简单脚本里把STUDY_AREA_M2做成独立常量不要混在GeoDataFrame的字段里计算完后顺手打印一组中间变量比如总共筛选出多少个斑块、绿地面积大概多少目测一下量级是不是符合直觉。我见过太多“看起来挺正常”的错误结果多打印一行中间变量就能救回来。最后分享一个我自己的操作习惯每次跑完密度计算我都会把输出表里附上计算参数清单包括MMU、投影坐标系、绿地类别代码、数据版本。参数比结果更容易被遗忘而几年后回看这批成果、或者和外部单位对口径时参数齐全能让沟通省下大把时间。再补充一点绿地斑块密度虽然好用但它永远只是描述格局的第一步最好和绿地覆盖度、平均斑块面积、边界密度摆在一起看才能把城市绿地的破碎化和分布特征讲清楚。 
长河君 · 2026/10/22 14:27:06 · 阅读 1.2k 绿地斑块密度计算全攻略：从图斑处理到Python自动化
做绿地格局分析的人大概率都遇到过这种场面同一份土地利用分类图两个人各自算绿地斑块密度结果一个0.8一个3.6谁都觉得自己没错。绿地斑块密度Patch DensityPD看起来就是一个数量除以一个面积但真正从图斑统计走到指数测算里面隔着坐标系、最小制图单元、统计口径、单位换算四道坎。这个指标又是城市绿地系统评估、生态破碎化分析、遥感解译后处理里最常被引用的基础指数之一算错了后面所有结论都跟着歪。这篇文章就把“绿地斑块密度计算”整条链路拆开讲先敲定定义和口径再处理原始数据得到干净的图斑接着走一遍手算过程最后给一个可以直接套用的Python自动化脚本并用一组虚拟数据复盘最容易翻车的三个地方。不管是做课程设计、毕业论文、规划项目还是科研数据处理照着这个流程走踩坑概率会小很多。1. 把“绿地斑块密度”的定义钉死口径不一致时一切对比都是空谈1.1 PD公式很简单但简单不等于不容易错斑块密度最常见的表述是单位景观面积内斑块的数量。公式写成PD N / A其中N是目标类别的斑块总个数A是景观总面积。常用单位是“个/平方千米”少数文献会用“个/公顷”或“个/100公顷”。换算关系很简单1平方千米100公顷所以“个/平方千米”和“个/100公顷”在数值上完全一致。这一点很多人不知道等会儿和Fragstats输出对照时就能用上。为什么这个公式容易错因为它把“斑块”和“面积”都当成了现成参数实际上“斑块N”取决于你怎么切分、怎么去碎斑“面积A”取决于你用什么分母。N和A只要有一个口径漂移PD就不是同一个量。我见过一份报告里写“绿地斑块密度2.5个/ha”另一份写“25.0个/km²”看起来数值差别大实际上这俩是同一个数。单位没写清楚才是真问题。1.2 分母用景观总面积还是绿地总面积差出数倍PD的分母可以是景观总面积也可以是该类别斑块总面积。前者衡量斑块在研究区内的绝对分布密度适合做不同区域间的横向对比后者衡量类别自身被切割成多少份更接近“破碎化程度”。绝大多数景观格局分析默认用前者但很多初学者会把分母写成绿地面积。用一个虚拟数字示例。某研究区景观总面积12.6平方千米绿地斑块178个绿地总面积3.62平方千米用景观总面积做分母PD 178 / 12.6 ≈ 14.1个/km²用绿地总面积做分母PD 178 / 3.62 ≈ 49.2个/km²同一个数据集分母一换结果差4.5倍左右。所以论文或报告的方法部分一定要写明“分母为景观总面积”。否则数据一进入横向比较整个结论都会失真。1.3 最小制图单元MMU决定斑块数量也决定密度量级比分母更隐蔽的是“斑块数量N”的取值。同一块区域用10米分辨率影像提取出的碎斑数量会远大于30米分辨率影像单个像元噪点也会形成独立图斑。这个问题的解决办法是设定一个最小制图单元Minimum Mapping Unit, MMU小于该面积的图斑不参与统计。MMU不是一个固定值取决于数据源和应用场景。生态类研究常用1000平方米量级城市街道尺度可能用几十到几百平方米。这个选择应该被写进方法描述里因为不同MMU会让“N”发生几十个百分点以上的变化PD自然跟着变。算密度前先想清楚“什么算一个斑块”比急着点计算器更重要。2. 从原始分类数据到干净图斑栅格转面、投影和碎斑处理2.1 第一步先建一张“绿地类别对照表”不同项目对“绿地”的定义差别很大。有的把林地、草地、灌丛都算绿地有的把城市公园绿地、防护绿地、附属绿地纳入有的还会把水体单独处理把农田单独拎出去。如果只是机械地把分类图里几个代码选出来后面所有统计都失去了可比性。我建议数据处理前先做一张对照表包含类别名称、原分类代码、是否属于本研究的绿地、备注。例如类别名称分类代码是否计入绿地备注乔木林地21是包含片林、行道树带灌木林地22是连续灌丛草地23是城市草坪、天然草地农田31否视研究目标可单列水体41否通常不计入绿地建设用地51否不透水面这张表建议写进项目文档。因为不同年份的数据如果分类代码略微调整有对照表才能快速同步口径也是后续自动化脚本里最需要维护的配置文件。2.2 栅格分类图转矢量图斑简化面选项要慎用如果手头是栅格分类图第一步通常是栅格转面。在ArcGIS里执行Raster to Polygon时有一个“简化面”相关的选项我的建议是保持关闭或者使用保留原始像元边界的模式。原因是简化面会把锯齿状边界抹平抹平过程中可能改变小斑块的形状甚至面积。转完矢量之后属性表里会多一个VALUE字段存放原始栅格值。这个字段就是类别代码之后所有筛选和分组都以它为主键。如果原始数据已经是矢量则需要先跑一遍“检查几何”和“修复几何”。面状数据里常见的自相交、重叠、缝隙在后续面积统计时会造成总和不准确早清理早消停。2.3 投影坐标系是算面积的前提经纬度坐标不能直接算这是个老生常谈却依然频繁出现的坑。WGS84经纬度坐标下直接算面积的操作结果只是“度²”量纲完全没有物理意义。正确的做法是把数据投影到等面积投影坐标系再计算面积。常用选择小范围用UTM对应分带或高斯-克吕格投影跨多个分带、做全国尺度的分析用Albers等面积投影。判断方法简单粗暴打开图层属性查看坐标系如果单位是“Degree”就一定要做Project如果单位是“Meter”才能继续算面积。一个特别容易忽略的细节是研究区边界和斑块图层必须使用同一套投影。否则用边界面积做分母、用斑块面积做分子单位都不能严格对齐算出来的PD自然不靠谱。2.4 小图斑的清理方式直接写入处理记录定义好MMU之后用属性查询把所有小于阈值的斑块选中并删除操作过程很简单。比如ArcGIS属性表里执行Shape_Area 1000 AND VALUE IN (21, 22, 23)如果绿地类别需要单列就先按类别筛出绿地再做面积过滤如果做全部地类的景观水平密度那么就先保留所有类别只按面积过滤。这里涉及一个研究层面问题删除碎斑和把碎斑合并进相邻大斑哪种处理更合理我个人的建议是如果只是计算密度和面积比例直接删除可复现性最强如果后续要做形态分析、连通性分析可能需要做融合处理。无论选哪种都要把处理方式写进方法文本否则别人无法对你的N进行审计。3. 图斑属性表里的统计博弈把数量、面积、边界长度规整到一张表3.1 面积字段怎么加不同软件的常规操作矢量斑块拿到手后首先给每个图斑添加面积字段。ArcGIS里可以用添加几何属性工具输出单位选“平方米”QGIS里用字段计算器表达式$areaPython里使用geopandas时每个要素的geometry.area会自动按当前投影坐标单位给出面积。值得提示的是如果数据来自不同机构有的面积字段可能叫Area_km2有的叫Area_ha还有的干脆没有。统一命名和单位在这一步一定要做我习惯全部统一成平方米字段名统一为area_m2。字段名一致后脚本可以批量处理多个文件。3.2 分组聚合统计斑块数量和类别总面积统计每个类别的斑块数用属性表的分组统计或数据库查询都能实现。ArcGIS里可以把属性表导出后用Excel透视表也可以直接用“汇总统计数据”工具分组。本质上都是对“类别代码”做分组聚合斑块个数N用COUNT统计类别总面积用SUM统计聚合结果是一张每个类别一行的表。关键来了这张表里的面积之和不一定是研究区总面积。因为研究区里可能有未分类区域、可能有编辑时的缝隙也可能因为数据叠加问题出现重复。计算PD时分母必须单独从研究区边界矢量上取这一条要刻在脑子里。3.3 单位换算一次推导永远不再搞混假设研究区总面积字段study_area_m2、某类斑块个数N、某类总面积area_m2那么斑块密度PD N / (study_area_m2 / 1000000)单位个/km²绿地覆盖比例PLAND area_m2 / study_area_m2 × 100单位%平均斑块面积MPS area_m2 / N单位m²前面说过1平方千米等于100公顷所以“个/km²”和“个/100公顷”数值相等。Fragstats输出的PD单位是每100公顷的斑块数读取结果时直接按“个/km²”理解即可不需要再做任何换算。这个小知识点能省下不少和软件输出结果大眼瞪小眼的时光。3.4 手工计算和Excel透视表的典型流程如果只算一次手动操作完全够用。典型步骤导出属性表到Excel→用透视表统计类别对应的行数和面积总和→把研究区面积作为常量放入单元格→用公式算出PD和PLAND→顺手加一列注释记录MMU、投影、数据版本。这个方法的好处是每一步都可见方便核对。缺点是一旦有几十个区域、几十个年份的数据透视表就要做几十次公式复制时容易漏引用。此时就该进入下一章的自动化方案。4. 自动化一步到位用Python把整条链路封装起来4.1 为什么要自动化口径固定比省时间更重要手工算一次只要十几分钟如果数据量小谈不上瓶颈。但景观格局分析通常不是一次性的可能要对同一城市不同年份、同一片区不同方案反复测算。手工流程每重复一次就有一次漏掉碎斑过滤、错选字段的风险而且多人协作时各种操作习惯都会悄悄混进来。把流程写成脚本等于把“怎么算”这件事实物化、版本化。代码就是方法记录参数就是口径清单。以后数据更新重新跑一遍就能得到新结果还能确保和上一版结果严格同口径。我认为这才是“一步到位”的真实含义不是省略步骤而是把步骤固定下来。4.2 一个可直接套用的Python计算脚本以下脚本使用geopandas输入是已经转换到投影坐标系的矢量面文件比如landuse.shp。绿地类别代码通过GREEN_TYPES列表控制MMU通过MMU_M2控制研究区总面积通过STUDY_AREA_M2控制。import geopandas as gpd # ---- 参数区开始 ---- GREEN_TYPES [21, 22, 23] # 绿地类别代码按自己的对照表配置 MMU_M2 1000.0 # 最小制图单元单位平方米 STUDY_AREA_M2 12_600_000.0 # 研究区总面积单位平方米从研究区边界矢量计算而来 INPUT_FILE landuse.shp # 投影坐标系矢量面文件 EPSG_PROJ EPSG:32650 # 若原始数据无投影或为经纬度先做重投影 # ---- 参数区结束 ---- gdf gpd.read_file(INPUT_FILE) # 确保使用投影坐标系本文以UTM 50N为例 if gdf.crs and gdf.crs.is_geographic: gdf gdf.to_crs(EPSG_PROJ) # 添加面积字段并过滤小斑块 gdf[area_m2] gdf.geometry.area gdf gdf[(gdf[landuse_code].isin(GREEN_TYPES)) (gdf[area_m2] MMU_M2)].copy() # 按类别聚合斑块个数、面积总和 rho gdf.groupby(landuse_code)[area_m2].agg([count, sum]) rho rho.rename(columns{count: patch_count, sum: total_area_m2}) # 计算景观水平的常见密度指标 rho[pd_per_km2] rho[patch_count] / (STUDY_AREA_M2 / 1_000_000.0) rho[pland_percent] rho[total_area_m2] / STUDY_AREA_M2 * 100.0 rho[mean_patch_area_m2] rho[total_area_m2] / rho[patch_count] rho.reset_index(inplaceTrue) print(rho) # 如果需要总绿地口径 total_count rho[patch_count].sum() total_green_area rho[total_area_m2].sum() print(f总体斑块数{total_count}个绿地面积{total_green_area / 1_000_000.0:.2f}平方千米 fPD {total_count / (STUDY_AREA_M2 / 1_000_000.0):.2f}个/平方千米)实际使用时landuse_code字段名要和自己的数据一致绿地代码也要按前文对照表调整。脚本会自动滤掉小于MMU的碎斑并根据研究区总面积计算每个绿地类别的PD、绿地覆盖度、平均斑块面积最后输出总体数值。值得注意脚本把分母统一为STUDY_AREA_M2这就是口径固定最好的体现。4.3 不写代码也可以串成“一步到位”的模型不会Python的分析师可以用ArcGIS模型构建器或QGIS图形化建模器实现相同功能。节点顺序大概是栅格转面如果需要按属性筛选绿地图斑并过滤小面积斑块添加几何属性计算面积汇总统计分组求斑块数和面积添加研究区面积常量字段计算器计算PD、PLAND、MPS模型建好后每次只要替换输入图层和参数就能一键出结果。我仍然推荐同时维护一份Python脚本因为脚本更容易做批量循环和版本管理。两种方式并不冲突都是把重复劳动收敛成一个单元。4.4 用Fragstats交叉验证一次能救回不少低级错误脚本算完建议用Fragstats做一次抽查验证。操作要点输入原始栅格分类图把背景值Background设置为0或-9999以外的实际不参与分析的类别千万不要把非绿地区域全部设为NoData否则NoData不进入景观总面积分母会缩水到只剩下绿地面积PD数值直接膨胀数倍。Fragstats输出的PD如果和脚本结果对不上优先检查三处背景值是否一致、最小图斑是否在栅格预处理阶段被合并、面积单位是否正确。我实操中遇到的无数次“对不上”最后都能归到这三件事上数据处理环节的重合率非常高。5. 用一组虚拟数据把数字跑通复盘最容易翻车的三个地方5.1 两张对照表实战先看结果再解释意义我用一组模拟数据来演示输出形式。假设A城市高新区研究区面积12.6平方千米B城市新城区研究区面积18.4平方千米绿地类别都是林地、灌木、草地MMU设为1000平方米。处理后得到区域绿地斑块数绿地面积(km²)PD(个/km²)绿地覆盖度(%)平均斑块面积(m²)A城市高新区1783.6214.128.720337B城市新城区2457.3613.340.030041只看PD两个区域似乎接近一个14.1、一个13.3。但结合绿地覆盖度就发现B区绿地面积占比明显更高但平均斑块面积也更大。这说明A区的绿地破碎度相对更高B区有更多成片大绿地。如果只报一个PD值很容易把差距讲歪。密度类指标的价值在于配合覆盖度和平均斑块面积一起解读而不是单独看某个数字。5.2 翻车点一碎斑阈值一改数量直接翻倍同样用A城市高新区的模拟数据如果把MMU从1000平方米降为0也就是不过滤碎斑斑块数会从178飙到406PD从14.1变成32.2个/km²。这个变化完全不是因为区域发生了变化只是处理口径变了。这意味着什么任何跨时间、跨区域的PD比较都必须先锁定MMU。否则2020年的分析用了1000平方米2025年的分析用了100平方米哪怕真实现状没有太大改变PD也会呈现“假性倍增”。做这类指标分析时把MMU写进表头或者研究方法里是对自己后续结论的一种保护。5.3 翻车点二没有统一投影面积偏差悄悄混进结果前面提到经纬度不能直接算面积。实际操作里还有一种更隐蔽的口径漂移部分数据是地理坐标系部分数据是投影坐标系汇总时没做统一转换就用软件默认方式叠在一起。结果可能是研究区面积被系统自动换算斑块面积却没同步最终PD偏差几个百分点还很难察觉。批量处理多个文件时建议在脚本里预先判断坐标参考系若为地理坐标系就强制重投影。前面Python代码里已经演示了如何检查。这句gdf.crs.is_geographic实际上是一个很简单的护栏能拦下大量面积统计事故。5.4 翻车点三分母从这个变量悄悄换成了那个变量最后一个坑来自复制代码时的粗心。假设某次计算把“研究区总面积”写成了“绿地总面积”178个斑块除以3.62平方千米得到49.2个/km²和正确值14.1完全不在一个量级。这种错不会让程序报错输出看着也很合理恰恰最难发现。规避方法很简单脚本里把STUDY_AREA_M2做成独立常量不要混在GeoDataFrame的字段里计算完后顺手打印一组中间变量比如总共筛选出多少个斑块、绿地面积大概多少目测一下量级是不是符合直觉。我见过太多“看起来挺正常”的错误结果多打印一行中间变量就能救回来。最后分享一个我自己的操作习惯每次跑完密度计算我都会把输出表里附上计算参数清单包括MMU、投影坐标系、绿地类别代码、数据版本。参数比结果更容易被遗忘而几年后回看这批成果、或者和外部单位对口径时参数齐全能让沟通省下大把时间。再补充一点绿地斑块密度虽然好用但它永远只是描述格局的第一步最好和绿地覆盖度、平均斑块面积、边界密度摆在一起看才能把城市绿地的破碎化和分布特征讲清楚。