第三章 空间数据获取与预处理：坐标系统、数据清洗、投影变换
做了这么多年空间数据分析，我最大的体会是：数据准备阶段决定了整个项目的成败。模型再牛，算法再先进，原始数据乱七八糟，结果能好到哪去？今天咱们就聊聊空间数据预处理那几个绕不开的坎——坐标系统、数据清洗、投影变换。说白了，就是把你手里的原始数据，变成能直接喂给模型的好数据。
3.1 坐标系统——空间数据的"地址簿"
坐标系统是什么？就是给地球上的点一个精确的门牌号。问题是，这套门牌号有好多种写法。WGS84、CGCS2000、西安80、北京54……新手经常搞混。
我个人习惯，拿到数据的第一件事——先看坐标系元数据。别急着分析，先搞清楚你的点到底在哪儿。
地理坐标系用经纬度，单位是度；投影坐标系用米，单位是米。两者之间差了投影变换。你想想看，把椭球面的经纬度直接当平面坐标用，距离量出来全是错的。
| 类型 | 常见代表 | 单位 | 适用场景 |
|---|---|---|---|
| 地理坐标系 | WGS84、CGCS2000 | 度（°） | 全球定位、底图显示 |
| 投影坐标系 | UTM 50N、Gauss-Kruger | 米（m） | 面积量算、缓冲区分析 |
我在项目中遇到过一档子事：合作伙伴给的数据，说是WGS84，我一算距离，偏差了四百多米。后来发现他们把CGCS2000当WGS84用了。两个坐标系在大部分地区差异不到一米，但在高精度场景下，那就是灾难。
print(gdf.crs) 看一眼，确认坐标系全不全、对不对。
下面这张图，就是本章的核心知识体系。你可以把它当作地图，指导后面的学习路径。
3.2 数据清洗——别让脏数据毁了你的分析
空间数据尤其容易脏。GPS信号漂移、传感器故障、人工录入错误……什么情况都有。我处理过的数据集中，几乎每一批都有问题。
常见的脏数据有三种：
- 缺失坐标——经纬度字段为空，或者填了 0,0（这个在（0,0）点在非洲西海岸，明显不对）
- 重复点——同一位置被反复记录，特别是轨迹数据里特别常见
- 离群点——坐标值突然跳到几百公里外，多半是GPS受到干扰
来看看我是怎么处理的。下面这段代码，基本是我每个项目都会用到的清洗模板。
import pandas as pd
import numpy as np
# 加载原始数据
df = pd.read_csv("raw_poi_data.csv")
# 1. 移除缺失坐标
df = df.dropna(subset=["lat", "lon"])
df = df[(df["lat"] != 0) & (df["lon"] != 0)]
# 2. 重复点检测（保留第一条）
df = df.drop_duplicates(subset=["lat", "lon"])
# 3. 离群点过滤（基于3-sigma原则）
lat_mean, lat_std = df["lat"].mean(), df["lat"].std()
lon_mean, lon_std = df["lon"].mean(), df["lon"].std()
mask_lat = np.abs(df["lat"] - lat_mean) < 3 * lat_std
mask_lon = np.abs(df["lon"] - lon_mean) < 3 * lon_std
df_clean = df[mask_lat & mask_lon]
print(f"原始记录数：{len(df)}")
print(f"清洗后记录数：{len(df_clean)}")我曾经在一条共享单车轨迹数据上吃过亏——没做离群点过滤，结果模型训练出来，预测的骑行距离平均偏大 40%。后来发现是几个点漂移到了隔壁城市，把整体均值拉偏了。那之后，我每次必做离群点检测。
你可能会问："过滤掉了会不会丢失信息？" 会。所以我的做法是：把过滤掉的点单独存一个文件，留作核查。别直接扔掉，万一那是新出现的模式呢？
3.3 投影变换——把地球"压平"的艺术
地球是个椭球体。你不可能把橘子皮铺平而不撕裂。投影变换就是解决这个问题的——用一种数学方法，把椭球面转换到平面上。
不同的投影方法，适合不同的区域和分析目的：
| 投影方法 | 特点 | 最佳用途 |
|---|---|---|
| UTM 投影 | 分带投影，每带变形小 | 中等尺度区域分析 |
| Albers 等面积 | 面积保持准确 | 面积统计、密度计算 |
| Mercator | 角度保真，面积变形大 | 航海、底图可视化 |
| Lambert 等角 | 角度保真，适合中纬度 | 气象、航空 |
我个人建议：
- 做距离、缓冲区分析 → 用 UTM 或 Gauss-Kruger（单位是米）
- 做面积统计 → 用 Albers 等面积投影
- 只是画个底图看看 → 用 Web Mercator（EPSG:3857）
下面是一个用 geopandas 做投影变换的例子。注意看，代码里指定了源坐标系和目标坐标系。
import geopandas as gpd
# 读取数据（假设是 WGS84 经纬度）
gdf = gpd.read_file("samples.geojson")
print("原始 CRS:", gdf.crs)
# 重投影到 UTM 50N（适用于东经114°-120°区域）
gdf_utm = gdf.to_crs("EPSG:32650")
print("投影后 CRS:", gdf_utm.crs)
# 现在可以算距离了
gdf_utm["distance_km"] = gdf_utm.geometry.length / 1000
print(gdf_utm[["id", "distance_km"]].head())utm_zone = int((lon + 180) / 6) % 60 + 1 快速算出来。北半球用 326+带号，南半球用 327+带号。
3.4 综合实战——从原始 CSV 到干净的空间数据集
把前面三个知识点串起来。假设你拿到一份商家 POI 的 CSV 文件，里面只有经纬度，目标是做成一个投影后的空间数据集，用于后续的聚类分析。
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
# 1. 读取原始数据
df = pd.read_csv("shops.csv")
# 2. 数据清洗
df = df.dropna(subset=["lat", "lon"])
df = df[(df["lat"] != 0) & (df["lon"] != 0)]
df = df.drop_duplicates(subset=["lat", "lon"])
# 3. 转为地理空间数据（WGS84）
geometry = [Point(xy) for xy in zip(df["lon"], df["lat"])]
gdf = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326")
# 4. 重投影到 UTM 50N
gdf_utm = gdf.to_crs("EPSG:32650")
# 5. 验证投影结果
print(gdf_utm.crs)
print(gdf_utm.head())
# 6. 保存为 GeoPackage
gdf_utm.to_file("shops_clean_utm.gpkg", layer="pois", driver="GPKG")shapely.geometry.Point 接收的是 (x, y)，也就是 (lon, lat)。搞反了，点会跑到海里去。
嗯，这里要注意：上面的代码用的是 Point(lon, lat)，先经度后纬度。这是地理空间领域的标准顺序。我见过不少项目因为这个顺序问题，分析做了半天，结果发现点全在非洲。
回顾一下这一章的核心：坐标系统是数据的"语言"，数据清洗是"去杂质"，投影变换是"换工具"。三者缺一不可。你想想看，如果坐标系搞错了，后面的分析就像用中文翻译器去读法语——结果全是乱码。