☰
栅格数据空间叠置分析全流程:原理、步骤与常见问题排查
2026/10/4 1:28:40 网站建设 项目流程

做栅格数据的空间叠置分析,看着是GIS里的基本功,但真到了具体题目或者实际业务里,坑是真不少。早几年我带项目组做某区域地下水位评估,拿到一份全国尺度的地下水埋深栅格,要和区县行政边界做叠置统计,结果一跑就出问题——要么统计出来的区县平均埋深明显不符合实际,要么边界处全是空值,折腾了整整两天才发现是环境设置里的"捕捉栅格"没指定,范围对齐出了偏差。

这篇文章我就从头到尾捋一遍栅格数据空间叠置分析的完整流程,把原理、步骤、参数设置、实战案例、问题排查全部讲透。不管你是刚学GIS的学生,还是做遥感、环境评价、城市规划,需要经常跟栅格打交道的人,这篇内容都适合你。我会结合一个地下水位栅格和行政区面数据做叠置提取的完整案例,把这些年踩过的坑一次性说清楚。

1. 叠置分析到底在解决什么问题

1.1 栅格叠置与矢量叠置的本质差异

先说个很多人容易混淆的点。矢量叠置和栅格叠置,虽然都叫"叠置分析",但背后的逻辑完全不同。

矢量叠置处理的是点、线、面这些几何对象,核心操作是求交、裁剪、合并、擦除,比如拿一块规划用地面去切土地利用图斑,本质上是几何运算,算的是交点、重合范围、拓扑关系。而栅格叠置处理的是覆盖整个研究区的像元阵列,核心操作是对每一个像元做算术运算、逻辑判断、条件筛选,比如"降雨量栅格"和"土地利用栅格"叠置,实际上是把两个栅格在同一像元位置上的值取出来做计算,得到一个新的值。

这个差异决定了工具选型完全不同。ArcGIS里矢量叠置用Overlay工具箱,栅格叠置则用Spatial Analyst里的栅格计算器、重分类、提取分析等工具。很多人一上来就在矢量工具里瞎找,方向就错了。

1.2 哪些业务场景离不开栅格叠置

栅格叠置分析的应用面非常广,几乎每个跟空间评价沾边的行业都会用到。

最典型的就是土地适宜性评价。比如某城市要做建设用地适宜性评价,需要将坡度栅格、地质灾害易发性栅格、生态保护区栅格、断层缓冲栅格叠加起来,按照权重计算综合得分——这个过程就是最经典的栅格叠置,而且通常涉及多图层加权叠加。

第二个常见场景是水文与环境分析。比如把一个区域的降水量栅格、土壤类型栅格、地形湿度指数栅格叠加,计算径流潜力或者水土流失风险等级。这几年热度很高的全国城市形态栅格数据集构建,本质上也是对夜间灯光、人口密度、建筑覆盖等多源栅格做叠置与分类。

第三个场景就是题目里提到的"地下水位栅格数据叠置行政区shp"。这种需求很典型,比如要算每个区县的平均地下水埋深、各行政区内的水位低于某个阈值的面积占比,就需要把地下水位栅格和行政区面数据叠置起来,提取并做分区统计。这类操作在环境监测、资源调查、灾害评估中几乎天天要用到。

明白了栅格叠置的应用场景,你会意识到一个问题:它不是一个单一工具能搞定的"点一下"操作,而是一条完整的分析流水线,从数据准备到参数设置,每一步错了结果就废了。

2. 动手前必须想清楚的三件事

2.1 坐标系不统一,后面全是白干

这是我在实际项目中踩过最狠的一个坑。有一次同事拿来自行下载的栅格数据,投影信息丢了,结果和当地的shapefile叠加时,图层显示在各不相关的位置,以为是软件坏了,最后检查才发现是坐标系不一致。

栅格叠置分析的前提是所有参与运算的数据都位于同一个坐标系下,否则像元与像元之间根本没有对齐关系,计算出来的结果毫无意义。具体来说,有三层要检查:

第一,检查投影信息是否完整。在ArcGIS里右键图层属性,切到"源"选项卡查看空间参考,如果显示Unknown,就要想办法补齐投影信息。如果数据本身是从正规渠道下载的,到元数据里找原始投影。如果实在找不回来,只能通过地理配准控制点来校正,但这是下下策,误差很大。

第二,统一使用投影坐标系还是地理坐标系。我的建议是:如果分析区域面积较小(比如一个城市、一个流域),统一使用投影坐标系,因为平面计算的面积、距离都是准确的。如果分析区域大(比如全国尺度的地下水位栅格统计),要么用适合全国范围的多圆锥投影或兰勃特等角圆锥投影,要么干脆统一到WGS84地理坐标系,但要清楚地理坐标系下栅格的单位是度,计算面积时需要额外转换。

第三,注意动态投影的陷阱。把两个坐标系不同的图层叠加显示时,ArcGIS会自动做动态投影,地图上看起来位置是吻合的,但栅格计算器执行运算时,如果数据源本身坐标系没有真正统一,动态投影反而会掩盖问题。所以我判断的标准很简单:不要在显示层面看对齐,要看数据源层面的坐标系定义是否一致。

2.2 分辨率、范围与像元对齐

坐标系统一了,第二个门槛是像元对齐。栅格数据的本质是一个规则的像元矩阵,两幅栅格要叠置,必须保证同样位置上对应的是同一片地面区域。

这里有三个参数必须强制一致:像元大小(Cell Size)、像元位置(Snap Raster)、分析范围(Analysis Extent)。

像元大小很容易理解,如果一幅栅格是30米分辨率,另一幅是90米分辨率,直接叠置时,90米的像元会对齐到30米像元的格网上,需要做重采样。这个我建议不要依赖工具的默认行为,而是显式地去设置。

像元位置最容易忽略。两幅30米分辨率的栅格,如果像元的起点不一致,叠置时边缘部分会出现半个像元宽的错位,或者数据被整体平移。解决方法是设置环境里的"捕捉栅格"(Snap Raster),让分析输出栅格的像元格网和参考栅格完全重合。

分析范围同理。两个数据范围不一致时,叠置结果的范围内会出现NoData区域,直接影响后续统计结果。我在ArcGIS中每次做叠置前,都会在环境设置里指定"处理范围",确保输出的栅格恰好覆盖目标范围。

2.3 重采样方法怎么选

前面提到分辨率不一致时要重采样,但重采样方法的选择是有讲究的。

常用的有四种:最邻近法(Nearest Neighbor)、双线性内插(Bilinear)、三次卷积(Cubic Convolution)和众数法(Majority)。

我的经验规则是:分类数据(土地利用类型、土壤类型、岩性、生态功能区)一律用最邻近法或众数法,因为双线性和三次卷积会在类别边界处生成介于两个类别之间的"插值值",比如土地利用类型0和1之间出现0.5,这没有物理意义。连续数据(高程、气温、降雨量、地下水位埋深)可以用双线性或三次卷积,它们会让表面更平滑。但要注意,三次卷积运算量大,而且有时候会在数值突变区域产生过冲,比如某地区地下水位在断层附近骤降,三次卷积插值后可能出现低于实际最低值的异常点。

所以我的默认方案是:所有用于叠置分析的分类栅格用最邻近,连续栅格优先用双线性。除非有明确理由,不要轻易改。

3. 完整操作流程:从数据准备到结果出图

3.1 工具选型:ArcGIS、QGIS还是Python

执行栅格叠置分析的工具有很多,我平时用的最多的是ArcGIS的Spatial Analyst扩展模块,尤其是栅格计算器(Raster Calculator)和地图代数(Map Algebra)表达式,因为直观、灵活,适合交互式探索分析。

QGIS也能做,内置的栅格计算器(Raster Calculator)功能几乎等价,而且开源免费,适合预算有限的项目。但QGIS在处理超大栅格时的性能比ArcGIS略逊一筹,尤其是做地图代数表达式涉及多图层运算时,内存管理上有点吃亏。

Python脚本(GDAL、rasterio、numpy)是我的杀手锏。当数据量很大、图层很多、或者需要批量处理几十个年份的栅格时,写脚本比在GUI里点来点去可靠得多。rasterio配合numpy做叠置分析非常顺手,本质就是读数组、做运算、写数组。GDAL则更偏底层,适合处理位深、压缩、坐标系转换等细节。

我的建议是:做一次性的探索分析,用ArcGIS或QGIS的栅格计算器;如果这个叠置分析要被反复执行,或者要嵌入到自动化流程中,果断用Python脚本。

3.2 七个步骤走通一次叠置分析

下面是我在实际项目里总结出来的标准流程,按这个顺序走,基本不会出大问题。

第一步,数据盘点。把所有参与叠置的栅格和面数据加载到工程里,逐个查看属性、坐标系、分辨率、范围、NoData情况。这个环节我通常会写一个简单的记录表,避免后面设置参数时才发现某图层有问题。

第二步,统一坐标系。把面数据(shp)和栅格数据都转到一个共同的坐标系下。我习惯先确定目标坐标系,然后统一通过Project工具把矢量投影过去,栅格用Project Raster工具重投影。

第三步,统一范围与像元。在环境设置里指定处理范围为研究区范围,指定捕捉栅格为基准栅格,指定输出像元大小。这一步是防止错位的核心。

第四步,重采样。如果各栅格分辨率不一致,用Resample工具重采样到统一的像元大小,分类数据选最邻近,连续数据选双线性。

第五步,执行叠置运算。根据分析需求,写地图代数表达式或者调用对应工具。这一步是核心,具体的计算方法要看场景。

第六步,验证结果。输出栅格后,叠加底图检查范围、对齐情况、值域是否合理。我一般会分三段验证:看最大值最小值是否在合理范围,看NoData占比是否异常,看空间分布是否与已知实际情况吻合。

第七步,输出与归档。把结果导出为带坐标系的GeoTIFF格式,并写好元数据说明,包括坐标系、分辨率、生成时间、参与图层。这一条很多新手忽略,但项目做多了你就知道,数据没有元数据,等于没有身份证明。

3.3 地图代数与条件函数的正确用法

栅格计算器看起来简单,但写表达式时容易犯低级错误。

ArcGIS里地图代数的基本语法是表达式写在"Map Algebra expression"框里,例如:

"precip_raster" * 0.4 + "slope_raster" * 0.3 + "soil_raster" * 0.3

这个表达式表示三幅栅格按照权重叠加。要注意的是,双引号里是图层名称,不能带路径,如果图层名有空格,需要用引号括起来。

更复杂一点的条件判断,用Con函数。比如统计地下水位埋深大于5米的区域:

Con("water_table" > 5, 1, 0)

意思是,当water_table像元值大于5时输出1,否则输出0。这样你就得到一幅二值栅格,接着可以统计面积为多少。

多层条件嵌套也是常见需求,比如:

Con("water_table" < 2, "high_risk", Con("water_table" < 5, "medium_risk", "low_risk"))

这里给每层Con都设置真/假两个分支,嵌套起来就能做多级分类。不过嵌套层数一多,表达式就难读,也容易出错。我的习惯是拆成多个中间栅格,分步计算,每一步验证一下中间结果,别追求一行写完。

还有一种常见操作是逻辑运算,比如选出同时满足两个条件的区域:

("landuse" == 1) AND ("water_table" > 3)

这种表达式的结果是1(真)和0(假),后续可以直接用重分类或者乘面积来计算符合条件的区域占比。

4. 实操案例:地下水位栅格与行政区划叠置提取

4.1 案例背景与数据说明

假设我们要解决一个典型的实际业务问题:获得某省地下水位埋深栅格数据(单位:米),需要统计每个区县的平均地下水位、最低水位以及水位埋深大于5米(说明地下水资源量较大)区域的面积占比。

读入数据前先梳理一下手头有什么:

  • 地下水位栅格数据:某省的连续表面,分辨率100米,坐标系为CGCS2000_3_Degree_GK_CM_117E(这里用这个例子),范围覆盖全省。
  • 行政区划shp:包含区县边界要素,属性表里有"省/市/县名称"字段。坐标系也是CGCS2000,但是地理坐标系GCS_CGCS2000。
  • 最终需要输出一份Excel统计表,以及一张分区统计专题图。

这个场景几乎是"地下水位栅格数据shp"热词背后的典型需求,很多人做这类工作一上来就想用"栅格转面然后和shp做相交",其实效率很低。正确的做法是叠置提取加分区统计。

4.2 详细操作步骤

我以ArcGIS为例,一步步说清楚。

第一步,坐标系先统一。前面的数据说明里已经看出问题:栅格是投影坐标系,shp是地理坐标系,所以必须先统一。我用Project工具,把shp从GCS_CGCS2000投影到CGCS2000_3_Degree_GK_CM_117E,与栅格一致。这里注意,如果你不确定用哪个投影带,就用栅格自己的坐标系,让矢量去迁就栅格,因为矢量投影转换不会损失数据精度,栅格重投影反而可能丢信息。

第二步,环境设置。打开ArcToolbox,右键"环境"设置:

  • 处理范围(Processing Extent):把"范围"下拉选为"与图层landuse相同",或者填自定义坐标范围。这里我选择与地下水位栅格相同。
  • 像元大小(Cell Size):选择"与图层water_table相同",保持100米。
  • 捕捉栅格(Snap Raster):选water_table图层。这一步是保持像元对齐的关键。
  • 坐标系:选择"与图层water_table相同"。

这三个环境参数一起设置,能保证后续所有栅格操作的输出格网完全一致。

第三步,执行掩膜提取。这一步相当于用行政区面数据去裁剪地下水位栅格,把研究区之外的像元裁掉。工具路径:"Spatial Analyst Tools" -> "Extraction" -> "Extract by Mask"。输入栅格选water_table,掩膜数据选投影后的区县shp,输出范围会自动限定在区县边界内。这里要注意,掩膜数据必须是面要素,而且只保留你真正需要的那些区县,如果shp里包含无关的行政区,建议先用选择工具筛选出来再输出一个新图层。

第四步,做分区统计。这是整个流程的核心。"Spatial Analyst Tools" -> "Zonal" -> "Zonal Statistics as Table"。参数这样设置:

  • 输入栅格或要素区域数据(Input raster or feature zone data):选择区县shp。
  • 区域字段(Zone field):选择区县名称字段。
  • 输入值栅格(Input value raster):选择提取出来的地下水位栅格。
  • 统计类型:选择ALL,这样会输出最小值、最大值、平均值、标准差等一套统计量。
  • 勾选"忽略NoData"(Ignore NoData in calculations)。

输出是一张属性表,包含每个区县的面积和各统计值。这张表可以直接导出为dBF或通过表格转Excel工具输出为xlsx。

第五步,计算大于5米区域的面积占比。这里不能直接用Zonal Statistics,因为它算不出"满足条件的像元面积"。需要先在栅格计算器里生成二值栅格:

Con("water_table_extracted" > 5, 1, 0)

然后对二值栅格做分区统计,统计类型用SUM,这样每个区县得到的SUM值就是满足条件的像元数量,乘以单个像元面积(100米×100米=10000平方米)就能换算成面积。把这个面积除以区县总面积就得到占比。

如果你想减少手动换算,也可以先知道一个像元的面积,在结果表中新建字段,用字段计算器写公式转换。

第六步,可视化与出图。把提取后的地下水位栅格用分级配色显示,基于区县边界做注记,最后在布局视图中加点图例、指北针、比例尺,导出图片。这里我建议用分位数分级或者自然断点分级,可以让区县之间的差异更清晰。

4.3 结果怎么看

拿到Zonal Statistics的输出表,第一件事不是急着汇报,而是要做合理性检验。

我一般按三步走:先看每个区县的像元数量和面积是否和实际行政区面积大致相符,面积差距超过10%就要回去查掩膜提取和坐标系的问题。再看平均水位值,跟当地水文部门的监测井实测数据做个抽样对比,如果个位数区县偏差过大,很可能是局部区域有数据空洞被NoData填充了,需要进一步处理。最后看标准差,如果某个区县的标准差异常大,说明该区县内部水位变化剧烈,这种情况下平均值的代表性要打问号,最好结合最小/最大水位一起分析。

另外还有一个小技巧:如果在输出表中发现某个区县完全没有记录,十有八九是区县边界和栅格范围没有真正相交,或者坐标转换后边界漂移出了栅格覆盖范围,检查一下这个区县的位置就明白了。

5. 常见问题与排查技巧实录

5.1 结果全为空或NoData

这是最经典的问题。现象是叠置后的栅格整个区域都是NoData,或者Zonal Statistics输出全为空。

排查顺序是这样的:

第一,查范围设置。我遇到最多的情况是环境里的处理范围和输入的掩膜范围压根没重叠,比如处理范围设置成了别的图层,导致输出区域其实在数据覆盖范围之外。把环境设置打开,确认范围没错。

第二,查NoData的定义。栅格数据里-9999、-99这些可能是"填充值"(Fill Value),而不是有效数值。如果这些填充值没有被正确标识为NoData,计算时会把-9999当作真实数值参与运算,结果自然惨不忍睹。

第三,查裁剪方向。面数据和栅格的坐标范围虽然重叠,但如果一个用的是经纬度坐标系,另一个用投影坐标系,且两者在数据层面的坐标系定义都写着"未知"或"WGS84",那么即使地图上看起来叠在一起,实际空间位置也可能偏了几百公里。出现结果全空,一定要回去检查坐标系,这是最根本的原因。

5.2 边缘锯齿和错位

现象是叠置结果栅格和实际边界有半个像元左右的错位,或者裁剪出的栅格边缘呈现锯齿状,但边界位置总是差一点。

这个问题的根源几乎都在于"捕捉栅格"没有设置。当你用两个不同来源的栅格做叠置时,如果没有指定Snap Raster,工具会用自己默认的像元格网起始位置,可能会和基准栅格错开半个像元。

解决办法:在环境设置里把"捕捉栅格"选为作为基准的那幅栅格。如果已经做过错了,重新设置后重新计算即可。这里提醒一下,Snap Raster和输出像元大小是配套的,改了捕捉栅格后要检查输出像元大小是否仍然正确。

另外一个容易被忽视的场景:用矢量面裁剪栅格时,面的边界本身就是折线,裁剪出来的栅格边界必然会出现锯齿。这是分辨率决定的,不是错误。如果出图需要平滑边界,后期可以稍微做一下平滑处理,但数据精度不能因此受影响。

5.3 内存溢出与运行缓慢

栅格叠置涉及大规模数值计算,尤其是多图层加权叠加或者高分辨率大范围数据时,内存不足和运行缓慢非常常见。

有三种处理思路:

第一,缩小处理范围,分块计算。不要一次性拿全国范围的100米分辨率栅格做运算,先用面数据按省或市裁成多个块,分别计算后再拼接。

第二,降低不必要的数据精度。如果做叠置分析的目的只是做趋势判断,可以考虑把浮点型栅格转为整型,或者用比较低的位深(比如8-bit unsigned),这样内存占用会大幅下降。但要注意,转整型前要明确是否会损失所需精度。

第三,改用Python脚本,用rasterio的windowed read/write功能分块读取和写入。我自己处理过好几GB的全国地下水埋深栅格,靠这个方式轻松跑完,而且出问题时还能轻松调试。

5.4 栅格裁剪不掉面数据

有朋友问我"栅格裁剪掉面数据"到底怎么理解,其实这里指的是用面数据去裁剪栅格,结果发现裁剪完的栅格范围并没有被限制在面以内。

这个问题的常见原因有三个。一是你用的裁剪工具是"Clip"而不是"Extract by Mask",Clip工具按范围框裁切,只认外接矩形,所以矩形之外、面之内的部分会变成NoData,或者面之外但矩形之内的部分反而被保留了。正确做法是Extract by Mask,它按面要素精确裁剪。二是掩膜面数据里有多余的、范围超大的要素,比如你想裁剪某个县,但shp里还包含整个省的轮廓,裁剪结果就被省的范围覆盖了。三是掩膜面本身有多个部分,或者带有洞,裁剪结果会遵循这些拓扑细节,如果觉得结果不对,先检查掩膜要素的几何结构。

6. 写在最后的实操心得

做栅格叠置分析这么多年,我最大的体会是:这个活儿技术门槛不高,但对细节的抠挖程度极高。欧洲的同行说过一句话我觉得很有道理:"GIS分析的结果,10%取决于算法,90%取决于你对数据的理解。"坐标系、分辨率、NoData、环境参数——这些看起来琐碎的东西,才是真正决定分析成败的变量。

在我个人操作中,最实用的一个小习惯是:每次做栅格叠置前,先在草稿纸上把参与图层的信息写成一张清单,包括坐标系、分辨率、范围、NoData设置。画两分钟,能省下来半天排查问题的时间。

另外再分享一个扩展方向。如果你处理的是常年连续的地下水位栅格数据,比如十年的月均值栅格,完全可以把叠置分析批量化。用Python配合rasterio写一个循环,把每年的水位栅格都按行政区做分区统计,最终输出一张"区县×年份"的水位统计表,这样构建出来的时间序列数据集,价值远远超过单次分析。很多热词中提到的"全国城市形态栅格数据集""地下水位栅格数据shp",本质上都是这种多源多时相栅格叠置组合的产物。掌握了基础逻辑,你完全可以按照自己的需求去构建这类数据集。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询