第六十八 位置

位置（Location）是万物在地球上的坐标。从 GPS 卫星定位到手机地图导航，从外卖配送调度到共享单车围栏，从物流追踪到气象预报，几乎所有现代应用都离不开“在哪里“这个问题。而“在哪里“的答案，离不开坐标系——尤其是 WGS-84 与 GCJ-02 这一对“相爱相杀“的坐标体系。

一、地理坐标基础

1.1 经纬度

地球表面的位置用**经度（Longitude）和纬度（Latitude）**表示：

   北纬 90°（北极）
        │
  西经 0°┼──────────────► 东经 180°
   -180° │  本初子午线
        │
   南纬 -90°（南极）

1.2 经纬度表示法

#![allow(unused)]
fn main() {
// DMS 转十进制度
fn dms_to_decimal(deg: f64, min: f64, sec: f64) -> f64 {
    deg + min / 60.0 + sec / 3600.0
}
// 116°24'26.6" → 116 + 24/60 + 26.6/3600 ≈ 116.4074
}

1.3 椭球体与基准面

地球不是正球体，而是近似椭球体。描述地球形状需要定义椭球参数：

不同的椭球体 + 基准面构成不同的坐标系，同一个经纬度在不同坐标系下对应地球上不同的点。

二、WGS-84 坐标系

2.1 什么是 WGS-84

WGS-84（World Geodetic System 1984，世界大地测量系统 1984）是 GPS 全球定位系统使用的坐标系，由美国国防部建立，是全球最通用的地理坐标基准。

2.2 GPS 如何定位

至少 4 颗卫星:
卫星位置（已知） + 信号传播时间 → 距离
   │
多颗卫星距离求交 → 三维坐标 (x, y, z)
   │
转换为经纬度 (WGS-84 系)

2.3 WGS-84 的应用

• 全球定位与导航（车载、手机、船舶、航空）
• 科学测绘、大地测量
• 国际标准数据交换
• 开源地图（OpenStreetMap）

三、GCJ-02 坐标系

3.1 为什么需要 GCJ-02

GCJ-02（国测局坐标，俗称“火星坐标系“）是中国国家测绘局发布的地理坐标系，对 WGS-84 坐标施加了非线性加密偏移（通常几十到几百米）。中国境内公开地图产品必须使用 GCJ-02，这是法律法规要求，也是地图厂商（高德、腾讯等）默认使用的坐标系。

WGS-84（真实位置） ──加密偏移──> GCJ-02（火星坐标）
      │                              │
   GPS 直接读出                   地图上显示的位置
   与地图不符                       与真实位置有偏移

在中国境内使用 GPS 模块直接读取的 WGS-84 坐标叠加到高德/腾讯地图上，会发现位置偏移几十米到几百米——这就是坐标系不一致造成的。

3.2 加密偏移原理

GCJ-02 偏移并非简单的固定平移，而是与经纬度相关的非线性函数：

GCJ-02 纬度 = WGS-84 纬度 + 偏移量(lat, lng)
GCJ-02 经度 = WGS-84 经度 + 偏移量(lat, lng)

偏移量随位置变化: 北京、上海、新疆的偏移方向和大小都不同

3.3 BD-09（百度坐标）

百度地图在 GCJ-02 基础上又做了一次加密偏移，形成 BD-09：

WGS-84 ──国测局加密──> GCJ-02 ──百度加密──> BD-09
  (GPS)            (高德/腾讯)         (百度地图)

3.4 坐标转换的必要性

四、坐标转换算法

4.1 WGS-84 → GCJ-02（示意核心）

#![allow(unused)]
fn main() {
/// WGS-84 转 GCJ-02（火星坐标）核心算法（简化示意）
const A: f64 = 6378245.0;        // 长半轴
const EE: f64 = 0.00669342162296594323; // 偏心率平方

fn out_of_china(lng: f64, lat: f64) -> bool {
    !(72.004 <= lng && lng <= 137.8347 && 0.8293 <= lat && lat <= 55.8271)
}

fn transform_lat(x: f64, y: f64) -> f64 {
    let mut ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y
        + 0.1 * x * y + 0.2 * (x.abs().sqrt());
    ret += (20.0 * (6.0 * x * std::f64::consts::PI).sin()
        + 20.0 * (2.0 * x * std::f64::consts::PI).sin()) * 2.0 / 3.0;
    ret += (20.0 * (y * std::f64::consts::PI).sin()
        + 40.0 * (y / 3.0 * std::f64::consts::PI).sin()) * 2.0 / 3.0;
    ret += (160.0 * (y / 12.0 * std::f64::consts::PI).sin()
        + 320.0 * (y * std::f64::consts::PI / 30.0).sin()) * 2.0 / 3.0;
    ret
}

fn transform_lng(x: f64, y: f64) -> f64 {
    let mut ret = 300.0 + x + 2.0 * y + 0.1 * x * x
        + 0.1 * x * y + 0.1 * (x.abs().sqrt());
    ret += (20.0 * (6.0 * x * std::f64::consts::PI).sin()
        + 20.0 * (2.0 * x * std::f64::consts::PI).sin()) * 2.0 / 3.0;
    ret += (20.0 * (x * std::f64::consts::PI).sin()
        + 40.0 * (x / 3.0 * std::f64::consts::PI).sin()) * 2.0 / 3.0;
    ret += (150.0 * (x / 12.0 * std::f64::consts::PI).sin()
        + 300.0 * (x / 30.0 * std::f64::consts::PI).sin()) * 2.0 / 3.0;
    ret
}

/// 返回 (GCJ-02 经度, GCJ-02 纬度)
fn wgs84_to_gcj02(wgs_lng: f64, wgs_lat: f64) -> (f64, f64) {
    if out_of_china(wgs_lng, wgs_lat) {
        return (wgs_lng, wgs_lat); // 境外不偏移
    }
    let d_lat = transform_lat(wgs_lng - 105.0, wgs_lat - 35.0);
    let d_lng = transform_lng(wgs_lng - 105.0, wgs_lat - 35.0);
    let rad_lat = wgs_lat / 180.0 * std::f64::consts::PI;
    let magic = (rad_lat).sin();
    let magic = 1.0 - EE * magic * magic;
    let sqrt_magic = magic.sqrt();
    let d_lat = (d_lat * 180.0) / ((A * (1.0 - EE)) / (magic * sqrt_magic) * std::f64::consts::PI);
    let d_lng = (d_lng * 180.0) / (A / sqrt_magic * rad_lat.cos() * std::f64::consts::PI);
    (wgs_lng + d_lng, wgs_lat + d_lat)
}
}

4.2 GCJ-02 → BD-09

#![allow(unused)]
fn main() {
/// GCJ-02 转 BD-09（百度坐标）
fn gcj02_to_bd09(gcj_lng: f64, gcj_lat: f64) -> (f64, f64) {
    let x_pi = std::f64::consts::PI * 3000.0 / 180.0;
    let z = (gcj_lng * gcj_lng + gcj_lat * gcj_lat).sqrt()
        + 0.00002 * (gcj_lat * x_pi).sin();
    let theta = (gcj_lat).atan2(gcj_lng) + 0.000003 * (gcj_lng * x_pi).cos();
    let bd_lng = z * theta.cos() + 0.0065;
    let bd_lat = z * theta.sin() + 0.006;
    (bd_lng, bd_lat)
}
}

4.3 反向转换

GCJ-02 → WGS-84 没有解析逆函数，通常用迭代逼近：

猜测 wgs = gcj
循环:
    gcj' = wgs84_to_gcj02(wgs)
    wgs += (gcj - gcj')   // 用误差修正猜测
直到误差足够小

注意：坐标转换应遵循各地法律法规；中国境内地图产品必须使用 GCJ-02。

五、距离计算

5.1 Haversine 公式

球面上两点间的大圆距离：

$$d = 2R \cdot \arcsin\left(\sqrt{\sin^2\left(\frac{\Delta\phi}{2}\right) + \cos\phi_1 \cdot \cos\phi_2 \cdot \sin^2\left(\frac{\Delta\lambda}{2}\right)}\right)$$

const EARTH_RADIUS_M: f64 = 6371000.0;

/// 两点间的球面距离（米）
fn haversine(lat1: f64, lng1: f64, lat2: f64, lng2: f64) -> f64 {
    let to_rad = |d: f64| d * std::f64::consts::PI / 180.0;
    let d_lat = to_rad(lat2 - lat1);
    let d_lng = to_rad(lng2 - lng1);
    let a = (d_lat / 2.0).sin().powi(2)
        + to_rad(lat1).cos() * to_rad(lat2).cos() * (d_lng / 2.0).sin().powi(2);
    2.0 * EARTH_RADIUS_M * a.sqrt().asin()
}

fn main() {
    // 北京天安门 → 上海外滩
    let dist = haversine(39.9042, 116.4074, 31.2397, 121.4998);
    println!("距离约 {:.0} 公里", dist / 1000.0);
}

5.2 各距离算法对比

5.3 经验法则

纬度一度距离基本恒定，经度一度距离随纬度余弦缩小，所以“经纬度各差 0.01“不代表实际距离相同。

六、空间数据格式

6.1 GeoJSON

GeoJSON 是地理数据的标准 JSON 格式，被绝大多数地图服务与数据库支持：

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": { "name": "天安门" },
      "geometry": {
        "type": "Point",
        "coordinates": [116.4074, 39.9042]
      }
    },
    {
      "type": "Feature",
      "properties": { "name": "故宫范围" },
      "geometry": {
        "type": "Polygon",
        "coordinates": [[[116.397, 39.915], [116.407, 39.915],
                        [116.407, 39.909], [116.397, 39.909],
                        [116.397, 39.915]]]
      }
    }
  ]
}

6.2 常见几何类型

6.3 WKT 与 WKB

七、数据库存储

7.1 空间数据库对比

7.2 PostGIS 示例

-- 建表（带空间字段）
CREATE TABLE shops (
    id SERIAL PRIMARY KEY,
    name TEXT,
    location GEOMETRY(Point, 4326)  -- 4326 = WGS-84
);

-- 插入数据
INSERT INTO shops (name, location)
VALUES ('天安门', ST_SetSRID(ST_MakePoint(116.4074, 39.9042), 4326));

-- 查询附近 3 公里内的门店（按距离排序）
SELECT name, ST_DistanceSphere(location,
       ST_SetSRID(ST_MakePoint(116.4100, 39.9000), 4326)) AS dist
FROM shops
WHERE ST_DWithin(
    location,
    ST_SetSRID(ST_MakePoint(116.4100, 39.9000), 4326),
    3000)
ORDER BY dist;

7.3 空间索引原理：GeoHash

GeoHash 把经纬度编码为字符串，前缀相同的点地理上相邻，从而把二维搜索降为一维前缀匹配：

北京: wx4g0ec19x3d
上海: wtw3sj5x4u5  （前缀不同 → 相距远）

WGS-84: 116.4074, 39.9042
→ GeoHash: wx4g0ec19x3d（精度约米级）

#![allow(unused)]
fn main() {
// GeoHash 编码示意（二进制交错经纬度）
fn geohash_encode(lng: f64, lat: f64, precision: usize) -> String {
    const BASE32: &[u8] = b"0123456789bcdefghjkmnpqrstuvwxyz";
    let mut bits = String::new();
    let mut lng_range = (-180.0, 180.0);
    let mut lat_range = (-90.0, 90.0);
    for i in 0..(precision * 5) {
        if i % 2 == 0 {
            // 经度位
            let mid = (lng_range.0 + lng_range.1) / 2.0;
            if lng > mid { bits.push('1'); lng_range.0 = mid; }
            else { bits.push('0'); lng_range.1 = mid; }
        } else {
            // 纬度位
            let mid = (lat_range.0 + lat_range.1) / 2.0;
            if lat > mid { bits.push('1'); lat_range.0 = mid; }
            else { bits.push('0'); lat_range.1 = mid; }
        }
    }
    bits.as_bytes()
        .chunks(5)
        .map(|chunk| {
            let mut v = 0u8;
            for &b in chunk { v = (v << 1) | (b - b'0'); }
            BASE32[v as usize] as char
        })
        .collect()
}
}

7.4 Redis GEO

# 添加位置
GEOADD shops 116.4074 39.9042 天安门
GEOADD shops 121.4998 31.2397 外滩

# 查询附近 5 公里
GEOSEARCH shops FROMLONLAT 116.4100 39.9000 BYRADIUS 5 km ASC

八、地图软件与坐标系

8.1 国内主流地图

8.2 国际地图

8.3 地图坐标系使用要点

九、Rust 地理生态

/*
[dependencies]
geo = "0.32.0"
*/
use geo::{Point, Haversine, Distance};

fn main() {
    let a = Point::new(-72.1235, 42.3521);
    let b = Point::new(-72.1260, 42.45);
    
    // 使用 Haversine 度量空间计算距离
    let d = Haversine.distance(a, b);
    
    println!("Distance: {} meters", d);
}

十、总结与练习

本章小结

练习建议

1. 坐标转换工具：用 Rust 实现 WGS-84 ↔ GCJ-02 ↔ BD-09 完整转换，用北京、上海、乌鲁木齐三个城市验证偏移量差异。
2. 距离计算：实现 Haversine 与 Vincenty 算法，比较两个相距 1000 公里城市的结果差异。
3. GeoHash 实现：实现 GeoHash 编解码，验证“前缀越长、相邻点前缀越相似“的特性。
4. PostGIS 实战：安装 PostgreSQL + PostGIS，导入一批点数据，实现“附近门店“查询并对比有无空间索引的性能。
5. GeoJSON 可视化：用 geojson crate 生成一个包含点、线、面的 GeoJSON 文件，在 OpenStreetMap 或任意地图工具中可视化。
6. 坐标系调研：用手机 GPS 与高德地图对比同一地点的坐标差异，估算 GCJ-02 在你所在城市的偏移量。