前一阵子整理旧项目资料翻到一张2000年的土地利用程度网格化专题图当时为了这张图折腾了小一个月。这种图跟常见的按行政区划涂颜色的专题地图完全不是一个思路它把研究区切成一个个大小一致的规则格子再把每个格子里土地利用的程度算成一个数涂到格子面上。做出来之后一眼就能看出开发强度从市区向外怎么衰减、山前地带和平原怎么过渡比传统县域色块图的信息量高出一大截。下面我把这张“2000年石家庄土地利用程度网格化图”的完整思路和操作细节拆开讲。会涉及土地利用数据从哪来、土地利用程度指数怎么算、格网怎么做、图怎么出好看也把我当时踩过的坑一并交代。如果你正在做土地利用、区域规划或者资源环境类的分析制图这套流程可以直接拿去参考。1. 为什么要把专题地图“切成格子”从行政区划图到网格图的思路转变1.1 传统按县区统计的专题图内部差异被“平均”掉了大部分人对专题地图的认知停留在“一个县一个色块”的阶段。比如要表达某区域的土地利用程度常规做法是先按县区统计建设用地占比、耕地占比算出一个综合值然后填到县级行政区划面上。这种图的优点是好读、好解释但问题也很明显行政单元内部的空间差异完全被抹平了。举一个很典型的情况。石家庄西部山区和东部平原在同一个县里可能各占一半西部以林地和未利用地为主东部大片耕地甚至还有零散的建设用地。按县统计时这个县的计算结果落在中间段位出图之后整个县一个颜色看不出来西边和东边完全是两套土地利用结构。这种由行政边界引发的“均质化错觉”对空间分析来说是致命的。我当时做图时也想过用乡镇级别替代县级别来缓解问题但乡镇之间的面积差异仍然很大大乡镇和小乡镇一个值空间可比性还是不够。1.2 网格化地图的核心逻辑让分析单元规则化网格化专题图的做法是把研究区按固定的经纬度间距或地面距离间距划分为规则格网比如1公里乘1公里、5公里乘5公里然后用每个格网内部的地类面积构成来计算一个“土地利用程度指数”最终形成一张由成百上千个小格子组成的专题图。规则格网的好处首先是比较单元完全一致。1公里乘1公里的格网在石家庄全市范围内不管落在市区还是山区面积都是确定的这样算出来的指数才可以直接互相比较。其次是空间分辨能力明显提升。行政区划图只能回答“哪个县程度高”网格图能回答“这个县内部的哪个方向程度高、哪个方向程度低”甚至能看出沿河流和交通线的带状变化。这类思路在人口密度图、经济产出图上早就是常规手段很多社会经济的格网化产品都是这样制作的。土地利用程度做网格化本质上是用同一套空间化逻辑处理地表覆盖特征只不过把“人口数”换成“土地利用程度综合指数”而已。1.3 什么场景适合网格化什么场景不适合网格化不是万能的。网格单元是人为划出来的天然与行政管理脱钩如果做的是耕地占补平衡、建设用地指标分配这类涉及行政辖区口径的话题网格图不适合作为唯一依据因为它不直接对应街道、乡镇或村界。适合的场景集中在空间格局分析和多图层叠加研究比如做土地利用程度与海拔的关系、与道路距离的关系或者做两期指数变化的空间扩散分析网格图都是很好用的载体。网格单元大小可控、属性表干净导出后能直接进入统计软件和Python做下一步分析这是行政区划面数据做不到的。2. 数据准备2000年这期土地利用数据得先落到自己手里2.1 两条获取路径现成数据产品与自己解译做2000年这一时间节点的土地利用分析数据获取有两条路。第一条路是直接用公开的全国土地利用数据产品。这类产品通常基于Landsat影像解译生成比例尺约为1:10万分类体系完整包含耕地、林地、草地、水域、城乡工矿居民用地和未利用地六大类已经做过度量校正和拓扑检查。拿到数据后只需要裁剪出石家庄范围再做一些投影转换就能进入下一步。优点是省时间缺点是数据时相、分类标准和具体来源需要自己确认清楚不能拿来就用。第二条路是自己从遥感影像解译。2000年前后可以获取的Landsat 5 TM和Landsat 7 ETM影像空间分辨率30米用来做1:10万比例尺的土地利用分类没有问题。我在那次项目中用的主要数据源就是这两颗卫星的影像经过挑选后选了夏季生长季、云量低于10%的若干景拼接后裁剪出研究区。两条路对最终网格化制图的影响主要在分类细节。现成产品的地类多边形一般比较平滑自己解译的成果通常更碎、边界更毛糙在后续格网统计时会让邻近格网的指数差异变大。如果追求制图美观建议对分类结果做一个轻度平滑处理。2.2 影像预处理链路几何校正、辐射定标、大气校正自己解译影像的话预处理这关躲不掉。几何校正的目的是让影像坐标与实地位置匹配通常以地形图或更高分辨率的正射影像为基准控制点要均匀分布在整景影像范围内校正误差控制在0.5到1个像元以内。如果影像已经经过系统级几何校正这一步可以省很多力气但也需要验证边界和道路是否匹配。辐射定标和大气校正决定了地表反射率算得准不准。简单说是把传感器记录的DN值先转成辐射亮度再通过大气模型消除大气散射和吸收的影响。山区阴影在夏季影像上很常见如果跳过大气校正直接分类阴影区很容易被错分成水体或建设用地这个误差后期很难消除。影像融合、拼接和裁剪属于常规操作。拼接时要注意不同景影像之间的色调差异多时相之间存在季节差异或光照差异最好先做直方图匹配再拼接否则分类时的训练样本会互相干扰。2.3 土地利用分类方案、样本与精度验证分类体系我建议采用通用的二级分类框架一级六大类耕地、林地、草地、水域、城乡工矿居民用地、未利用地。这个体系下各类地物的光谱可分性最好也方便与公开产品衔接。如果想细化可以在耕地里分水田和旱地但一级六大类已经能满足土地利用程度指数计算需求。分类方法我用的是监督分类加人工目视修正。监督分类的流程是先通过影像和实地经验确定训练样区每类地物至少选取十几块样区覆盖该类型在研究区内的主要分布范围然后跑分类器。最大似然法是常用选择学起来快、结果稳定。分类完成后必须做分类后处理。原始分类结果中有一类很常见的现象叫“椒盐噪声”也就是单个像元被随机分到错误类别像一粒粒的椒盐散布在图上。可以用多数滤波或者众数滤波采用3乘3窗口将孤立像元归并到周围主导类别这样分类图看起来干净很多后续格网统计也更稳定。最后是精度验证。我在研究区内按分层随机方式抽样每类不少于50个检验点以原始影像目视判读结果作为参考生成混淆矩阵计算总体精度和Kappa系数。总体精度85%以上、Kappa不低于0.7是比较通用的合格线。如果不达标需要回到训练样区环节补充样本或修正误分区域这一步不能跳。3. 土地利用程度综合指数一个公式把“开发强度”算明白3.1 指数背后的设计思想土地利用程度综合指数不是简单统计某个地类的面积占比而是试图回答一个更抽象的问题一片土地上的人类活动改造强度有多大。这个设计思路是把土地利用分成四个等级从未利用地到城镇建设用地等级越高代表人类对地表系统的影响越强。未利用地基本没有被改造林草地和水域属于半自然状态农业用地意味着大面积人工干预城镇建设用地则是完全人工化的地表。这种分级方式在实际计算中带来一个很好的性质它把六大类地物压缩成一个可比较的数值避免“这个格网耕地多、那个格网林地多、没法直接比”的困境。你要说清楚一片区域到底是农区还是林区可能需要看两个指标但要说清楚开发强度一个指数就够了。3.2 分级赋值与计算公式国内做土地利用程度研究时最常用的分级赋值体系按地表被人类利用的强度划分成四级每一级对应一个指数值土地利用类型分级指数未利用地裸岩、沙地、盐碱地等1林地、草地、水域2耕地、园地、人工草地3城镇、农村居民点、工矿、交通用地4计算模型写成公式就是L 100 × Σ( Ai × Ci )其中A i 是第i级土地利用程度分级指数C i 是第i级土地利用面积占格网面积的百分比。这样算出来的指数落在100到400之间100表示格网内全部是未利用地400表示格网内全部是建设用地数字越大开发强度越高。举一个例子。假设某个5公里格网内耕地占60%林地占20%水域占10%农村居民点占10%那么这个格网的土地利用程度指数就是L 100 × (3×0.6 2×0.2 2×0.1 4×0.1) 100 × (1.8 0.4 0.2 0.4) 280这个280落在“中等偏高”区间说明该格网以农业利用为主、伴有少量建设用地干扰。如果换成同一个格网里耕地占20%、建设用地占50%、水域占30%指数就会变成100×(3×0.24×0.52×0.3)320明显高于前者。这就是综合指数的压缩表达能力。3.3 计算时要用面积占比不是面积绝对值网格面积相等时用面积绝对值算和用面积占比算结果一样。但实际工作中研究区边界会切掉边缘的格网有些格子是不完整形状面积不再相等。如果忽略这一点直接用各地类面积乘以分级指数再加和会放大边缘区域的影响。正确做法是每个格网先计算总有效面积再算各地类占比。比如一个落在边界上的格网理论面积25平方公里但真正落在研究区内的只有18平方公里其中林地9平方公里、耕地9平方公里那么占比均为50%指数按100×(2×0.53×0.5)250计算。如果不按占比而是直接按面积结果就会和完整格网完全不可比。3.4 为什么不直接用建设用地占比有人可能会问既然建设用地等级最高那么直接用建设用地占比来出图不就行了这确实是最简单的思路但会丢掉太多信息。一个以耕地为主的平原区和一个人为干扰严重的林区建设用地占比都接近零画出来几乎一个颜色而实际上两者的开发强度差别明显。综合指数的价值在于把不同等级的利用类型都纳入计算耕地、林地、建设用地各占多少都能反映到最终数值上。土地利用程度不是一个“有没有建设用地”的二值问题而是连续的地表人工化梯度。这个连续梯度正是网格化专题图最想表达的东西。4. 网格化制图全流程实操从格网生成到专题图输出4.1 格网尺寸怎么选先想清楚要表达什么尺度网格大小是整个制图流程里最需要权衡的参数。用1公里格网图上会出现大量高低值相间的小碎块细节丰富但噪声也大用5公里格网图面平滑但城郊接合部的过渡细节会被掩盖。我做这一版图时采用5公里格网作为主图。原因是石家庄2000年的土地利用格局整体比较宏观西部山区与东部平原的差异是地块尺度以上就能看出来的5公里格网足以表达而且网格数量适中ArcGIS处理和出图渲染都很流畅也方便后续做统计分析。如果你关心的是城乡接合部的扩张方向建议用1到2公里的格网。经验法则是格网边长不应小于原始土地利用数据最小图斑边长的10倍以上否则大量格网内只有零星地类指数波动很剧烈图上会出现明显的“盐碱化”斑点。更稳妥的办法是做两套格网对比比如1公里和5公里各出一张观察格局是否稳定如果核心格局一致说明结果可信。4.2 生成格网前先统一坐标系统这是最容易踩坑的环节我单独拿出来强调。生成格网前所有图层必须统一到同一个投影坐标系。我建议使用适合中国区域的等面积投影或高斯-克里格投影而不是用WGS84经纬度坐标系直接画网格。原因很简单经纬度坐标下1度纬度和1度经度对应的地面距离不同直接生成经纬度格网会得到一个“北窄南宽”或“东窄西宽”的变形格网。这类格网计算面积时误差极大后面的指数计算会跟着错。我的做法是把土地利用矢量数据统一投影到阿尔伯斯等面积投影或西安80高斯投影下再以投影坐标定义格网尺寸和模板范围。在ArcGIS里生成规则格网用的是Create Fishnet工具。设置Fishnet Origin为研究区范围的左下角Y-Axis Coordinate为左下角正上方一点Cell Width和Cell Height设为格网边长Template Extent选择研究区边界这样生成的格网范围与研究区基本吻合。4.3 空间连接算清每个格网里都有什么地类格网生成后下一步要把土地利用数据“灌进”格网里。我用的是Tabulate Intersection工具输入Zone要素为格网图层Zone Field选格网ID字段Class要素为土地利用图层Class Field选地类代码。这个工具会输出一张表每一行是“某个格网ID 某个地类代码 面积”单位是投影坐标下的平方米。这一步输出的是明细表还不能直接用于计算。需要把同一格网内各地类的面积汇总成一行地类代码变成列名再计算各地类占比。习惯用Excel做透视表的人可以在Excel里完成我在项目后期改用Python脚本处理速度更快也更方便重复使用。import pandas as pd # 读取Tabulate Intersection导出的表格 df pd.read_csv(grid_landuse_area.csv) # 按格网ID和地类代码汇总面积 pivot df.pivot_table( indexgrid_id, columnslanduse_code, valuesarea_m2, aggfuncsum, fill_value0 ) # 各类地物对应的分级指数 grade {1: 3, 2: 2, 3: 2, 4: 2, 5: 4, 6: 1} # 计算每个格网的土地利用程度指数 def calc_landuse_index(row): total_area row.sum() if total_area 0: return None weighted_sum 0 for code, area in row.items(): if code in grade: weighted_sum grade[code] * area return 100 * weighted_sum / total_area pivot[landuse_index] pivot.apply(calc_landuse_index, axis1) pivot.to_csv(grid_landuse_index.csv)计算完成后得到一张包含格网ID和指数值的表。在ArcGIS里右键格网图层选择连接表把这张表关联回格网属性表导出为新的要素类网格图就具备了按指数符号化的条件。4.4 分级设色别把指数图做成土地利用分类图符号化的逻辑和土地利用分类图完全不同。土地利用分类图的绿色代表林地、蓝色代表水体颜色有语义对应关系土地利用程度指数图表达的是连续变化强度应该用同一色相从浅到深递进或者用浅色到深棕色的渐变表示从低到高。我推荐的方案是自然间断点分级法分5级。自然间断点通过算法寻找数据分布的天然断点让级内差异最小、级间差异最大比等间距更适合表达指数这种非均匀分布的数据。色带选择浅黄色到深棕色的渐变低值浅、高值深直观且不会和地物分类色混淆。分级后要调整图例文字建议直接显示指数数值区间比如“220-260”并注明格网尺度与数据时相。图名设为“2000年石家庄土地利用程度网格化图”图例标题写“土地利用程度综合指数”。这样一张图拿到手里不需要额外解释就能看懂。4.5 布局与输出比例尺、指北针、网格线、DPI制图布局环节容易被忽略但它决定图件能不能直接放进报告或论文。切换到布局视图后加入比例尺和指北针是基础操作但要注意比例尺必须基于格网的实际投影坐标不能凭感觉拉一个默认比例尺。图面背景建议用浅灰色或白色格网边界线设成细线、浅色避免网格线喧宾夺主。研究区外围的格网应当裁剪掉保留完整的部分或者通过透明度让边缘自然地淡出。导出图片时分辨率至少设置300 DPI。很多人用默认96 DPI导出图面文字在打印后发虚。用ArcGIS导出地图时选择“质量”优先图片格式推荐TIFF或高质量PNG这样后期放入文档不会失真。5. 成图之后怎么看2000年石家庄土地利用程度的空间格局解读5.1 总体格局西低东高中间有极值核心图出来后第一眼看到的就是土地利用程度从西向东升高的总体梯度。西部山区大片格网以林地和草地为主指数普遍落在220到260区间部分以裸岩和未利用地为主的低山地带甚至不到220是全图最低的区域。中部平原的格网以耕地为主导指数稳定在280到320之间整片区域颜色比较均匀。这种均匀性说明当时这一带土地利用结构趋同大面积连片耕地构成了平原地区的基本底色。指数最高的区域集中在中心城区及其周边格网指数普遍超过350局部紧邻核心区的格网接近400。这个区域虽然格网数量不多但在图面上形成一个深色团块非常醒目一眼就能看出城市核心的位置和大致范围。5.2 值得细看的局部异常河流廊道与外围延展除了总体梯度图上有两个局部现象值得停下来分析。第一个是流经市域北部的河流及其沿岸滩地。河流水面和河漫滩在分类中分别属于水域和未利用地分级指数偏低所以沿河一带的格网指数被明显拉低在平原高值背景中形成一条浅色条带。这个条带如果做两期对比还会告诉你河流沿线的开发利用强度有没有随时间变化。第二个现象是中心城区以外的“指状”高值延展。仔细观察这些高值格网的分布方向基本沿对外交通干线往四周辐射路网沿线分布着农村居民点和乡镇企业导致沿线格网指数比两侧纯农业区高出一截。这种指状延展是城镇化早期的典型空间特征也是网格图相对行政区划图最能体现价值的地方县级统计完全看不出这种方向性。5.3 从静态图到变化分析网格图不只是用来“看”的网格化表达的最大优势在于属性表干净、单元规则方便做进一步定量分析。我当时做完2000年这幅图之后又做了2005年和2010年两期把三期指数栅格相减得到了土地利用程度变化图。正差值表示开发强度增强负差值表示开发利用退减图上能清楚看出城市扩展的主方向和各时期的热点区域。如果想继续深入还可以把格网指数与海拔、坡度、到主要道路的距离等因子做栅格叠加和回归分析或者计算全局和局部空间自相关指标检验指数是否存在显著的空间集聚。这些分析在行政单元上做容易出现尺度偏误但在规则格网上操作的合理性要高得多。6. 网格化地图最容易翻车的五个细节投影、边界、误差、配色与出图6.1 投影坐标系与面积统计这一条我在前面已经反复强调但仍然值得单独提醒。直接在经纬度坐标系下创建格网、计算面积结果会随纬度发生形变不同纬度的格网实际面积相差很多所有指数横向比较都会失真。务必在生成格网前把所有数据统一到投影坐标系并且用投影坐标下计算的面积进入指数公式。判断自己有没有踩坑有一个简单办法生成格网后抽查几个格网的Shape_Area字段按理论面积估算一下比如5公里完整格网应是25平方公里左右。如果面积数值明显不符合说明坐标系设置有问题回头检查数据框和要素类的投影定义。6.2 研究区边界上的不完整格网怎么处理规则格网与研究区边界相交时边界外的面积会被算进格网总面积。如果直接按原始格网面积计算占比边缘格网的指数会被严重低估因为大量空白面积拉低了各类地物的占比。处理方式有两种我一般选更保守的“相交面积大于50%才保留”。也就是先用格网与研究区求交得到实际落在研究区内的格网部分计算每个格网的有效面积如果有效面积小于该格网理论面积的50%就把这个格网删掉或标记为空值。另一种方式是保留所有与边界相交的格网但按有效面积计算占比这样边缘格网的数值也可用只是图面上会有些略微偏小的格子。两种方式都能接受关键是必须在成图说明里写清楚规则否则同一张图在不同人手里边缘处理方式不一样结果会有偏差。6.3 分类误差会顺着格网传递到指数上土地利用分类不可能做到100%准确尤其2000年的Landsat影像分辨率有限细窄道路、零星房屋经常漏分。当一个格网内只有三四种地类时某一类误分的面积占比可能达到10%以上对指数的影响是实打实的。格网越小这个问题越突出。一张5公里格网内有25个像元层级的地类混合少量误差会被平均掉1公里格网内可能只有几十个像元个别像元误分就能改变指数等级。如果坚持做小格网建议先对分类结果做Majority滤波再进行格网统计能明显减少椒盐噪声带来的指数波动。检验误差影响的一个办法是在分类图上随机修改一小部分图斑类别重算指数观察指数分级图是否剧烈变化。如果变化不大说明结果稳定如果某些格网跳了两个等级说明你需要更谨慎地对待分类精度。6.4 图例配色不能直接套用地类分类色这是一个看起来很小但很容易犯的错。土地利用分类图的默认配色是林地绿、水体蓝、建设用地红紫这些颜色有明确的语义关联。土地利用程度指数如果直接沿用分类色读者会不由自主地把绿色区域理解为林地把红色区域理解为建设用地从而误读整张图。指数图的正确配色逻辑是单一色相或双色渐变的连续色带。我用浅黄到深棕低值浅、高值深表达的是“开发强度从低到高”的渐变语义而不是某个具体地类。如果要用双色可以选浅蓝到深红但避免绿色系和蓝色系同时出现在非地类语义的场景里。图例分级数建议控制在5到7级太多会增加读图负担太少会丢失空间梯度信息。分类方法和分类间隔要在图例里写明比如“自然间断点分类”这样阅读者能判断分级是否合理。6.5 出图清晰度与要素完整度最后是导出环节。很多人做完图直接截图放到报告里结果图例文字模糊、边界锯齿明显。正确的做法是在布局视图排版加入完整的图名字体层级、比例尺、指北针、图例和数据来源说明然后以300 DPI以上分辨率导出。还有一个小细节容易被忽略格网专题图通常没有底图要素外围的行政边界和附近市县的名称就等于缺失了。如果读者不熟悉这个区域会很难定位图的覆盖范围。我当时在图上叠加了浅色的市级行政区边界线并在图框四角标注经纬度坐标这样读图者能快速判断出图区域的位置。写在最后的一点体会搞完这套流程我最深的感受是网格化专题图真正的价值不在于“把图做得不一样”而在于它逼着你把每一个空间单元都拉到同一把尺子下去衡量。行政边界可能掩盖差异平均数值可能抹掉结构但当你把一片区域切成一排排等大的格子每一格都有自己的土地利用程度指数时空间的连续性、过渡带的位置、高值中心的扩散方向全都清清楚楚。如果你准备在别的城市复现这个方法我的建议是先花半天时间把投影坐标系和数据分类标准理顺再开始画格网。这两件事做扎实了后面的指数计算和出图基本就是流水线操作。格网尺寸方面我也建议多做一套对比图看看结果对尺度敏不敏感这既是对自己数据质量的检验也是论文或报告中很加分的方法学细节。