网格化专题地图制图全流程：土地利用程度指数计算与可视化

网格化专题地图第一次看到成图的时候我确实被这种表达方式震住了。它不是以往那种“一个行政单位一块色块”的常规区划专题图而是把整个研究区像棋盘一样切成成百上千个等面积的格子每个格子用颜色告诉读者这一平方公里的土地被人类活动开发和利用到了什么程度。去年我在整理2000年华北平原某城市的土地利用数据时就用这套方法做了一张土地利用程度网格化图不少同行看了图都来问制作思路。这篇就完整复盘一下什么样的数据能喂进来、网格切多大最合适、土地利用程度指数怎么算、颜色和分级怎么设计以及我实际踩过的几个坑。无论你是GIS初学者还是想把传统专题图做出更多空间细节的从业者这套流程都可以直接照搬。1. 网格化专题地图的制图逻辑它凭什么“不一样”1.1 常规分级统计图的三块短板先说说我们最熟悉的“按行政区填色”也就是统计专题图里常说的分级统计图法。这种做法太常见了找一份县级或乡镇级行政边界把统计数据挂到属性表上然后用分级配色填色完事。没错它看起来整齐清爽尤其是数据齐、边界顺的时候放在报告里确实很体面。但如果你真的把图放大到局部去看问题马上就暴露了。第一个问题是行政单元内部的平均化。一个县的土地利用程度可能是县城周边极高、乡镇驻地中等、大片农业区偏低结果一平均整个县成了同一个颜色。你根本看不出来高强度开发的区域到底集中在哪更看不出城市辐射衰减的梯度。第二个问题是面积悬殊导致视觉误导。同一个市域范围内面积大的乡镇在视觉上天然占据主导哪怕它的指数不高大色块也会让人觉得它“很重要”。第三个问题是跨期不可比。行政边界隔几年就可能调整一次2000年的乡界和2020年往往不一样两期图放在一起做变化分析时边界一变动统计口径就乱了。这三点放在一起是我放弃常规做法的直接理由。我当时要呈现的并不是“哪个行政区强、哪个弱”而是土地利用程度在地理空间上是如何连续变化的、高强度开发的斑块到底分布在哪里。这种表达需求行政区划这种“人造边界”反而成了干扰。1.2 网格化地图用统一的“空间刻度”重新丈量土地网格化专题地图的思路很简单——不按行政边界统计而是把研究区用规则格网切碎每个网格成为独立的统计单元再对每个网格内部的地类面积比例进行综合计算最后按指数分级填色。这套逻辑最吸引我的地方在于它的“等积性”。每个格子的面积基本相同不存在大色块碾压小色块的视觉偏差网格不随行政区划变动而变化多年遥感数据做出来可以直接叠图对比更关键的是它能表达出行政边界完全展现不了的空间细节——城乡交错带上的梯度、沿着道路和河流的延伸、城市周边的建设用地扩散方向网格图上都非常直观。有人可能会问那直接用栅格数据做连续表面不是更细吗理论上是但多数业务场景不需要逐像元级别的连续表达。30米分辨率栅格图细节极多打印出来在小比例尺下反而成了一片噪点而且统计口径太细和土地利用类型斑块的实际语义对不上。网格化的尺度介于“像元”和“行政单元”之间既能看出空间分异又不至于被随机噪声干扰。1.3 土地利用程度指数把“开发强度”变成一个数要做土地利用程度网格化图核心指标就是土地利用程度综合指数。这个指数的思路简单说是把土地利用类型按照人类活动干扰强度的不同分成若干级别再按面积加权综合成一个数值。我采用的是比较经典的四级分类一级指数1未利用地级包括沙地、盐碱地、裸地、沼泽等自然状态几乎未受人工改造的地类二级指数2林草水级包括林地、草地、水域等受轻度利用但生态功能主导的地类三级指数3农业级包括耕地、园地等以农业生产为目标、人类持续干预的地类四级指数4城镇聚落级包括城镇建设用地、农村居民点、工矿用地等高强度开发的地类。综合指数的计算公式写出来是L 100 × ΣAᵢ × Cᵢ其中Aᵢ是分级指数1到4Cᵢ是第i级地类面积占网格总面积的比例取值在0到1之间。这个公式并不复杂。它实际上等于把每个网格里不同开发强度的地类按比例加权最后乘100是为了让数值更好读。如果一个网格全是未利用地L就是100全是城镇建设用地L就是400如果一半是耕地一半是林地L就是100×(3×0.52×0.5)250落在中段。数值越高代表受人类活动影响越强土地被开发利用的程度越高。这个指数的妙处在于它没有量纲上的突兀感区间稳定在100到400之间做成图例非常直观。而不像直接统计“建设用地面积占比”那样只反映单一地类它能体现一个网格里多种地类组合后的综合开发水平。2. 数据准备与核心参数选择2.1 土地利用数据的来源、格式与预处理做这类图的前提是先有一份可信的土地利用数据。我做的是2000年那一期当时主要用的是Landsat TM遥感影像的人机交互解译成果——逐像元识别地类、勾画图斑、拓扑建库最后交出的是带地类属性字段的矢量面图层。如果你手头是栅格分类数据比如30米分辨率的土地利用栅格也完全可以做网格化只不过统计方式要换成“分区统计”而不是“面积叠加”。但无论是矢量还是栅格有几点必须提前确认地类属性字段要规范。至少得有一个字段能区分出未利用地、林草水、农业、城镇聚落这四类如果原始分类更细也没关系事后在属性表里重分类归并即可。图形要素要干净。矢量面之间不能有大量重叠或缝隙否则后期做面积统计时会出现双算或漏算。建议在建图前用拓扑检查工具过一遍。时间口径要统一。2000年的数据就全用2000年前后的影像解译别混入2005年的补充图斑。网格化图对空间细节敏感混入异期数据会让局部格子值出现明显跳动。我当时拿到的数据是Shapefile格式的矢量面图层地类代码比较细比如水田、旱地、有林地、灌木林地、河流、城镇、农村居民点等大概十几个二级类。正式开工前我先在属性表里加了一个新字段“LU_Class”用条件赋值把二级类归并成1到4四个级别这样后面算指数时可以直接用。2.2 网格尺寸怎么选1公里和5公里的区别网格大小是整个制图流程里最需要斟酌的参数。切得太大空间细节被抹平网格图看起来就像打了马赛克的行政区图切得太小格子数量爆炸式增长制图渲染变慢不说边缘区域还会出现大量的碎斑块、噪声值读图体验反而变差。这里有一个尺度权衡的基本逻辑网格尺寸应当和“现象本身的尺度”匹配同时考虑数据源的分辨率。土地利用类型图斑一般最小的也有几十亩如果网格只有200米见方很可能一个格子只覆盖了半块农田或者一角村庄数据噪声极大如果网格大到10公里见方城市边缘的梯度就被掩盖了整个研究区也才几十个格子做出来跟统计图差别不大。我做2000年那个城市案例时先用5公里网格跑了一遍预览效果比较平滑但看不出城乡接合部的细节然后切成1公里网格细节出来了总格子数在8000左右单机处理仍然流畅。最终成图用的就是1公里。这里我可以给大家一个经验参数平原地区的城市型研究区1公里网格通常是最稳妥的起点如果是省域或流域级大范围5公里或10公里比较合适如果是中心城区的精细研究500米网格也可以尝试但要做好去噪的心理准备。另外有个小技巧网格尺寸最好取整数公里这样后期报告里写“以一平方公里为基本统计单元”会很清爽而且不同图幅间对比时口径也统一。2.3 坐标系与面积统计第一个连环坑这一点我放在参数部分说因为它真的是最容易出问题的环节而且出了问题往往藏得很深。做网格化图的整个统计过程无论是网格面积、地类面积还是面积占比全部依赖投影坐标系下的平面面积。如果你用的是WGS84经纬度坐标直接拿$area或者计算几何去求面积得到的不是一个可靠的地面面积值尤其在北方高纬度地区误差会大到离谱。正确做法是提前把数据投影到一个适合研究区的等积投影或高斯投影坐标系上。2000年那个案例用的是西安80坐标系的高斯-克吕格投影中央经线按研究区位置选定这样每个格子面积和地类面积的比例关系才准确。如果你用QGIS在“处理”之前先做一次“重投影图层”ArcGIS则用“投影”工具。还有一点必须提醒创建网格的范围、原数据、输出结果的坐标系必须保持一致。我见过不止一次有人建渔网时默认用了WGS84而土地利用数据是投影坐标系两套数据叠在一起空间关系完全对不上后面空间连接出来的统计结果当然全是错的。所以开始建网格之前先打开图层属性确认所有数据都在同一个坐标系里这一步花不了两分钟但能帮你避免后面浪费一晚上。2.4 属性字段设计与指数分级准备网格化图的属性表是整个分析的中枢建议在建网格之后就把字段结构规划好。我当时建了这么几个字段字段名类型说明MeshID整型网格编号方便后期排序和检索MeshArea双精度网格总面积单位平方米UnusedArea / ForestArea / AgriArea / BuiltArea双精度四个级别的土地利用面积单位平方米Ratio1 到 Ratio4双精度各级别面积占网格面积的比值LUDI双精度土地利用程度综合指数这组字段看着多其实用途很清晰。前几个是为了算面积比例最后一个是成图的分级依据。另有一件事我建议顺手做掉给每个网格算完指数后顺便建一个“等级”字段把指数值分成低、较低、中等、较高、高五档这样后续符号化渲染可以直接按等级字段来配颜色不用每次调整分级区间。3. 实操全流程从面状要素到网格化指数图3.1 第一步投影统一与数据体检正式干活之前先把数据整理到一个干净的状态。我在QGIS里的操作顺序是这样的在图层列表中加载土地利用矢量面和准备作为底图的行政边界。打开图层属性在“信息”页签确认当前坐标系。如果原始数据是WGS84就先使用工具箱里的“重投影图层”工具目标坐标系选适合当地的高斯投影或阿尔伯斯等积投影输出为新的Shapefile。这个重投影后的图层就是后续所有操作的唯一数据源。然后做一次简单体检打开属性表确认LU_Class字段已经被正确赋值1到4都有记录没有空值或异常值用“工具箱→矢量几何→检查有效”扫一遍看看有没有自相交、空几何之类的非法要素。如果有极少量的碎小图斑不用删因为网格统计时它们对面积比例的影响微乎其微但如果有大面积要素几何错误就一定要修复否则空间连接时会出现莫名其妙的漏算。3.2 第二步创建网格渔网在QGIS中创建网格用的是“矢量→研究工具→创建网格”这个工具英文界面叫Create Grid。参数设置里有几个关键选项网格类型选“矩形多边形”而不是“矩形线”或“点”因为我们需要的是多边形网格来做面积统计网格范围选择“继承自图层”并指定土地利用数据图层作为范围来源这样网格刚好覆盖整个研究区不会多出一圈无数据的空白格水平间距和垂直间距填10000米也就是1公里。注意这里的单位是地图单位必须在投影坐标系下才能直接填米输出的是一个新的多边形图层要素数量就是研究区面积除以1平方公里。建完网格之后我习惯先打开属性表看一眼总数。如果研究区大约8000平方公里网格数量应该接近8000。如果数量差得多比如只有几百个那多半是间距单位填错了——把10000米填成了10000度之类的低级错误。3.3 第三步空间连接与面积汇总网格有了接下来就是把土地利用信息“灌”进每个格子里。这一步我用的是按位置连接属性的方法也就是空间连接。在QGIS处理工具箱里找到“按位置连接属性”工具设置如下目标图层选网格图层连接图层选重投影后的土地利用矢量面几何谓词选“相交”只要土地利用面与网格有任何公共部分就纳入连接连接类型选“一对一”目标要素每个网格对应多条连接要素时把它们的属性合并起来。执行空间连接后每个网格会带着所有与之相交的地类要素的属性一起出现在输出表里。这里有个非常重要的操作因为一对多连接会生成多条记录直接在这个结果上算面积是不行的需要先做一步“按字段汇总”——把每个网格内、按LU_Class分别求面积总和。QGIS里可以用“汇总统计”工具分组字段选MeshID统计字段选面积分类依据是LU_Class也可以用“按属性值分组统计”直接得到每个网格、每个地类级别各自的面积合计。假如不想用空间连接后再汇总的路径也可以直接把网格图层和土地利用面做“相交”处理得到每个网格与每个地类面相交的碎块再按MeshID和LU_Class分组汇总面积。两条路都能走通我个人更推荐先行空间连接再汇总因为中间步骤少、逻辑更直观。3.4 第四步字段计算写出土地利用程度指数面积汇总完成后属性表里应该有了每个网格的四类面积。在此基础上继续加字段按比例和指数公式做计算。这一步在QGIS字段计算器里就能完成。先算面积比例。新增四个字段Ratio1到Ratio4分别用对应级别的面积除以网格总面积。比如Ratio4的表达式就是“BuiltArea / MeshArea”。这里要注意MeshArea必须用平面投影面积QGIS里可以新建字段用$area来计算经纬度坐标系下$area返回的不是真实地面面积这一步就废了。然后算综合指数LUDI表达式为100 * (1 * Ratio1 2 * Ratio2 3 * Ratio3 4 * Ratio4)表达式的含义就是把每个网格内四类土地的开发级别按面积比例加权再乘100定标。算完之后可以先做一次统计检查LUDI的最小值应该在100附近最大值不超过400。如果出现大面积数值为100或400的极端情况说明地类归并可能有问题或者网格尺度太小很多网格只覆盖了单一地类。3.5 ArcGIS下的等价操作如果你常用的是ArcGIS Desktop或ArcGIS Pro操作路径也对应得上。建网格用“数据管理工具→采样→创建渔网”设置模板范围时选择土地利用数据的图层范围像元宽度和高度填10000空间连接用“分析工具→叠加→空间连接”目标要素是网格连接要素是土地利用面连接操作选“一对一”或“一对多”后再汇总面积比例和指数计算用“字段计算器”表达式和QGIS里的写法基本一致。两套工具本质完全相同区别只是菜单位置和名称表述。对于不常用ArcGIS的读者优先推荐QGIS毕竟免费开源用起来没有授权压力而且处理几千个网格这种量级性能完全够用。4. 符号化与出图让网格图真正“读得懂”4.1 分级方法自然断点法还是等间距指数值算出来后到渲染成图那一步还有一个很影响成图效果的决策——怎么分级。直接拿连续数值给每个网格一个独一无二的颜色理论上信息量最大但打印出来或者缩放到小比例尺观看时人的视觉很难连续区分几十种过渡色图面反而显得杂乱。实用的做法是把LUDI分成5个等级。分级间隔可以用等间距也就是每60分一档也可以用自然断点法让类内差异最小、类间差异最大。我个人的经验是如果研究区土地利用结构比较简单、指数分布相对集中用等间距分级更稳定图例也好解释如果数据分布有明显的不均匀比如大量网格集中在中低值区域少数极高值网格把色带拉伸得很厉害就改用自然断点法。自然断点法的分级界线看起来不那么“整齐”但图面更符合数据本身的分布规律。5%左右的网格不管落在哪一档用自然断点法时都会自动被拉开。如果某个数值区间恰好是空档分级后依然会保留一级没有网格分布的颜色——没关系这是正常现象说明研究区里没有那一类强度的土地利用状态报告里还能反过来说明土地利用的空间分化特征。4.2 配色与图例设计网格图跟常规专题图一样颜色必须符合阅读直觉。土地利用程度这个指标直觉上是从“低”到“高”的单一维度我倾向于用从浅到深、从冷到暖的渐变色带来表达低值用浅黄或浅绿中值用橙色高值用深红。这样看图的人不需要读数值第一眼就能感觉到哪里开发强、哪里弱。具体配色上我参考了ColorBrewer色带用YlOrRd这一组浅黄→黄橙→橙→橙红→深红。它在绝大多数显示器、打印机上表现稳定而且对红绿色弱的人相对友好。无数据网格统一设置为浅灰色和正式数据区域形成对比避免被误读成低值区。图例需要注意的细节是必须标注指标名称和数值范围比如“土地利用程度综合指数”下面注明“指数越高表示受人类活动影响越强”。图例里最好也附上一句简短的分类说明例如“城镇聚落级地类占比高的网格取高值”这样非专业人士拿到图也能马上明白表达意图。4.3 地图整饰与信息完整性网格化专题图在打印和报告场景下地图整饰比普通专题图要求更高。因为网格图的小方块多、颜色变化细如果整饰不到位整个版面会显得拥挤和凌乱。我的版面基本元素是主图、图例、比例尺、指北针、数据来源说明、坐标系说明。主图尽量占版面主体图例放在右下或左下空白区域比例尺放在图的下方中央数据来源说明放在图的上角或下角。网格图因为没有行政边界线读者对空间位置的感知会弱一些所以我通常会叠加一个细线画的研究区边界轮廓再叠加几条县级或乡镇边界参考线用浅灰色显示既不影响网格主体又提供了方位参考。数据来源说明这一段千万别省。我当时写的是“土地利用数据源自2000年Landsat TM影像人机交互解译成果统计单元为1km×1km网格坐标系采用高斯-克吕格投影指数计算方法为土地利用程度综合指数模型”。这段信息看着是小事但决定了这张图能不能作为学术成果或规划依据被使用。4.4 进一步玩法转成栅格做变化检测网格化图的优势之一是可以和栅格数据结合。如果你在QGIS里把网格图层转成栅格按LUDI字段赋值分辨率设为1000米就得到了一张连续的指数表面栅格图。再用同样方法做2000年和2020年的两张网格指数栅格做一次栅格减法得到变化差值图哪里的土地利用程度提高了、提高了多少一目了然。这种做法的好处是突破了行政边界的天然限制。城市扩张、农村居民点扩展、生态退耕这些过程它们在空间上本来就不是按行政边界发生的用网格做统计单元才能把它们真实地表达成“空间过程”。这也是我为什么坚持做网格化而不愿意退回常规行政区专题图的最重要原因。5. 常见问题与排查技巧实录5.1 面积怎么和行政面积对不上这是网格化做出来后最容易碰到的疑问。研究区总面积明明是8000平方公里把所有网格面积加起来却少了或多了几十平方公里。原因一般有两个一是网格覆盖范围是矩形而研究区是自然边界矩形四个角上的网格一定有一部分越出研究区这些网格算总面积时仍然按完整矩形面积计算而地类面积只统计到研究区边界为止二是原始土地利用数据里本身可能有小面积空洞或未分类区域导致地类面积之和不足。这个现象是正常的不必纠结。网格化图的统计单元是“网格”不是“行政区”每个网格的统计口径就是网格范围超出研究区边界的部分不计入地类面积但依然占有网格面积——在做面积比例时分母必须用完整网格面积而不是研究区截断后的面积。只要你在报告里明确写清“统计单元为完整格网边缘网格存在部分无图区域”数据就是自洽的。不要为了让网格总面积等于行政区面积去修改网格面积那反而破坏了统计单元的一致性。5.2 网格和图层边缘对不齐出现细条碎格当你用“继承图层范围”创建网格时有时候生成结果和研究区的边界会有一点错位或者边缘出现一些极细长的碎块。这往往是因为网格起点默认从数据范围最小坐标开始而数据范围本身可能不是整公里倍数导致最后一行或一列只有半格宽度。要处理这个问题有两个思路。一是接受现状因为这些边缘网格面积占整个研究区比例很小对全局图面影响不大二是手动指定网格起点把创建网格工具里的“网格起点X坐标/起点Y坐标”设置为研究区合理范围内的整公里数让网格从规整的坐标开始铺这样边缘的碎格会少一些。更重要的是如果后期要做两个时期的对比建议两期的网格起点必须一致否则网格在空间上无法一一对应变化检测根本没法做。5.3 土地利用程度指数全是极端值图面像黑白噪声如果做完LUDI后发现大部分网格不是100就是400中间值很少做出来的图像撒了一把盐那基本可以判断网格尺寸过小了。格网一旦小于地类图斑的平均尺度很多格子就只能覆盖单一种类的地物指数自然只会取极端值。解决办法很简单把网格尺寸放大比如从1公里调到2.5公里或5公里重新跑一遍。还有一种情况是地类归并层级太粗把差异很大的地类压进了同一个级别比如把城镇、农村居民点全部归为城镇聚落级导致城镇周边网格几乎全取400。这时可以考虑把农村居民点和城镇区分开甚至临时引入二级权重调整指数的相对高低让中间梯度显现出来。5.4 空间连接后属性重复一个网格出现多条记录你可能会在空间连接的时候发现输出的表里同一个MeshID出现了好几行每个网格对应的是与之相交的多个地类面。这一步是正常的甚至是我们需要的中间状态——每个网格和多少块地类相交就会有多少条记录面积聚合要以MeshID分组完成。但如果你忘了分组聚合直接把连接结果看成“每个网格一行”后面算出来的指数就会错得离谱。我的经验是在空间连接步骤之后马上用“统计”工具按MeshID分组把各级地类面积合计进独立字段形成一个干净的一格一行的新表。后续所有计算都基于这个聚合表进行别再回头触碰原始连接结果。5.5 出图时网格边界太明显像马赛克网格化专题图的成图如果多边形边界线太粗或颜色太深视觉上会像马赛克干扰对颜色本身的阅读。我在最终出图时把网格边界线的颜色设为与当前色带相近、但更深的灰度线宽极细仅为0.1到0.2毫米更多时候直接把边界线设为无只保留颜色填充。如果你需要在地图上有“格子感”可以用一个很不起眼的浅灰色细线来画网格边界保证能看出统计单元的轮廓又不至于盖过填色信息。另外在导出图片时建议把输出分辨率设置到300dpi以上矢量格式PDF作为备选这样网格线不会因为压缩产生锯齿打印出来也更干净。6. 一点个人体会这套流程做完之后我最大的感触是网格化专题图并不是什么高深的技术更多是制图思维上的转换——从“行政单元思维”转向“空间连续思维”。同样是土地利用程度这个指标用行政区表达和用网格表达讲出来的故事完全不同。网格图里能看到城市周边开发的梯度衰减能看到河流沿线的绿色廊道也能看到村镇居民点像星星一样散落在农业平原上这些信息在常规区划图里几乎都会被平均掉。最后分享一个小技巧做网格化图时别急着追求“完美参数”。先用粗网格跑通全流程确认数据和逻辑都没问题再逐渐加密网格比较效果。因为网格大小直接影响总格子数量和后续所有计算量先用5公里网格把流程走通1公里网格去出成图能省下不少试错时间。按这套方法做出来的土地利用程度网格化图放到报告里或者学术论文里都会比千篇一律的行政区划图多出明显的空间说服力。