空间数据清洗与预处理

做空间数据分析这些年，我踩过最多的坑不在建模环节，反而在数据准备阶段。说白了一句话：垃圾进，垃圾出。你算法再漂亮，拿到一份缺胳膊断腿、拓扑关系全乱的数据，结果也是废的。

这一章咱们就专门聊聊数据清洗的四个硬骨头：缺失值处理、异常值检测、数据一致性检查、拓扑错误修复。我个人的习惯是，拿到任何空间数据集，先按这四个维度过一遍，再谈后续分析。这样能省掉至少一半的返工时间。

缺失值处理

空间数据里的缺失值，和普通表格数据不一样。它可能是某个栅格像元没有值（NoData），也可能是矢量数据的属性字段是空的。我记得有一次做城市绿地分析，到手的数据有一半的绿地斑块缺少“面积”属性，当时差点崩溃。

处理缺失值，常用的方法就这么几类：

• 删除法：直接扔掉有缺失的记录或像元。如果你缺失比例小于5%，且缺失是随机的，这招最省事。
• 均值/中位数插补：用周围邻居的均值填进去。对栅格数据尤其好用，像 ArcGIS 的 Nibble 工具就是干这个的。
• 克立金插值（Kriging）：适合有空间自相关特征的数据。我曾为土壤采样点的缺失值做过一次普通克立金插值，效果比简单均值好太多。

下面是一个 QGIS 中使用栅格计算器做简单邻域插补的脚本片段：

# QGIS Python 控制台：栅格 NoData 填充
import processing
from qgis.core import QgsRasterLayer

# 假设你已经加载了一个名为 "dem.tif" 的 DEM 图层
input_raster = "dem.tif"
output_raster = "dem_filled.tif"

# 使用 Fill NoData 工具，基于周围有效像元
processing.run("gdal:fillnodata", {
    'INPUT': input_raster,
    'BAND': 1,
    'DISTANCE': 10,
    'OUTPUT': output_raster
})

异常值检测

空间数据里的异常值往往不按常理出牌。你想想看，一个高程点突然比周围高出100米，或者一个湖泊的周长面积比离谱到爆——这些都要揪出来。

我常用的检测手段有以下几种：

方法	适用场景	说明
箱线图（IQR）	属性值异常	超过1.5倍四分位距的视为异常
空间局部Z值	空间分布异常	计算每个点与其邻居的标准化差值
几何规则检查	矢量图形异常	如狭长多边形、自相交环
领域上下文检查	专题数据	土地利用类型与周围不协调

记得有一次做土地利用变化分析，发现一块“森林”斑块的大小在5年内暴涨了20倍，结合遥感影像一看，原来是当年分类错误——把新造林也算进去了。这就是典型的异常值藏在属性里的案例。

用 Python 检测空间数据的局部异常，可以这么写：

import geopandas as gpd
import numpy as np
from shapely.geometry import Point

# 读取空间数据
gdf = gpd.read_file("building_footprints.shp")

# 计算每个面要素的面积
gdf['area'] = gdf.geometry.area

# 用 IQR 方法检测面积异常值
Q1 = gdf['area'].quantile(0.25)
Q3 = gdf['area'].quantile(0.75)
IQR = Q3 - Q1

lower_bound = Q1 - 1.5 * IQR
upper_bound = Q3 + 1.5 * IQR

outliers = gdf[(gdf['area'] < lower_bound) | (gdf['area'] > upper_bound)]
print(f"检测到 {len(outliers)} 个面积异常要素")

数据一致性检查

一致性检查说白了就是看数据内部的“逻辑”对不对。我给你列几个最常见的检查项：

• 属性域一致性：字段取值是否在合理范围内，比如“土地利用类型”不能填“3.14”这类浮点数。
• 几何合法性：多边形的环方向是否正确？有没有自相交？
• 时空一致性：同一位置不同时段的属性变化是否合理，比如1990年的城区面积不可能比2020年还大。

这里有一张我整理的检查清单表格，大家可以照着做：

检查项	检查方式	典型问题
字段类型	对比 Schema 定义	数字字段存了字符串
唯一值检查	统计 distinct count	ID 字段重复
空间参考	对比 EPSG 编码	多个图层 CRS 不一致
范围检查	Min/Max 统计	高程出现负数
逻辑关系	字段间比对	面积小于0或大于行政边界总面积

# 用 GeoPandas 检查属性域
valid_landuse = [1, 2, 3, 4, 5]  # 假设有效编码
gdf = gpd.read_file("landuse.shp")

inconsistent = gdf[~gdf['landuse_code'].isin(valid_landuse)]
print(f"发现 {len(inconsistent)} 条无效土地利用编码")

# 输出异常记录到独立文件
inconsistent.to_file("landuse_inconsistent.shp")

拓扑错误修复

拓扑错误是矢量数据的老大难。面重叠、缝隙、悬线、伪节点——这些错误不修掉，后续的叠加分析、空间连接全是错的。

我经历过的惨痛教训：有一次做城市路网分析，拓扑检查做了三遍才把悬线全清掉。为什么？因为部分路段的端点差了0.001个度，肉眼根本看不出来，但网络分析就是跑不通。

常见的拓扑错误类型：

1. 悬线（Dangles）：线要素端点没有与其他线连接。多出现在交通网络数据。
2. 伪节点（Pseudo Nodes）：一个节点只有两条线连接，这通常是数字化错误。
3. 缝隙（Gaps）：相邻面之间有空隙，多出现在行政区或地块数据。
4. 重叠（Overlaps）：面要素交叉覆盖，常见于多源数据融合时。

下面是一个用 PostGIS 修复拓扑错误的通用流程：

-- 1. 创建拓扑结构
SELECT CreateTopology('my_topology', 4326);

-- 2. 将要素加入拓扑层
SELECT TopoGeo_AddPolygon('my_topology', 'land_parcels', 'geom', 'gid');

-- 3. 检查拓扑错误
SELECT * FROM my_topology.edge WHERE ST_IsValid(geom) = false;

-- 4. 自动修复悬线和缝隙
-- 注意：ST_Snap 可以将不匹配的节点吸附到容差范围内
UPDATE land_parcels
SET geom = ST_Snap(geom, ref_geom, 0.001)
FROM (SELECT ST_Collect(geom) AS ref_geom FROM land_parcels) AS ref;

做个总结：数据清洗没有捷径，但有流程可依。我的习惯是：miss-check-consistency-fix，四步走一遍，再用空间可视化扫一眼。如果你赶时间，至少要把拓扑错误修掉——这是很多后续分析出错的根本原因。