实验 16：Z 序曲线（Z-order / Morton code）
为什么做这个实验
数据库里最好用、最成熟的索引是 B+ 树，可它只能对**一个**排序键做范围查找。
"找出某个矩形框内的所有点"是个二维问题，一个 B+ 树区间表达不了。
绕过去的办法是找一个映射
f: (x, y) → z（一个整数）
让**空间上邻近的点，在一维顺序里也尽量邻近**。这样二维的矩形查询就能近似变成"若干个 z 值区间"，
直接复用现成的 B+ 树。这类映射叫**空间填充曲线**（space-filling curve，SFC）；
Z 序（Z-order，也叫 Morton 码）是其中构造最简单的一种——只需要把两个坐标的二进制位**交织**起来。
本实验要弄清楚四件事：
1. 位交织怎么算、怎么反算（它是严格的一一映射，可逆）；
2. 为什么"坐标高位相同 ⇒ 排序后仍邻近"，也就是 Z 序局部性的来源；
3. 用普通一维 B+ 树索引二维数据的**两阶段流程**：先取候选、再精确过滤；
4. 在仓库自带的真实经纬度数据上，把 Z 序与顺序扫描、与 Hilbert 曲线的**候选率**量出来，
并说清楚这个对比的方法学前提。
一句话预告结论：Z 序不是"最好"的曲线，但它是**最便宜**的曲线——位交织一次就算完，
工业界大量系统（Geohash、S2 的部分环节）仍在用它。
运行前准备
• Python 3.9 及以上。
• 必选依赖：`numpy`（整数版 Morton 码与真实数据读取）、`matplotlib`（画图）。
• 可选依赖：`ipycanvas`（本 notebook 最后一个 cell 用它画交互式画布）。
装不上也没关系，改用下面给出的 matplotlib 版本即可。
• 缺包时安装：
```bash
python -m pip install numpy matplotlib
```
（本机已配置环境可用：`python -m pip install numpy matplotlib`）
• 数据文件：`notebooks/latitude_longitude.bin`（仓库自带，197 229 个经纬度点）。
读取数据时要 `os.chdir` 到仓库根目录，或用绝对路径。
• 打开方式：在仓库根目录执行 `jupyter notebook notebooks/Z-Order-Curve.ipynb`。
**关于字符串实现**：下面的教学版用字符串表示 z 码，一眼能看出每一位从哪来。
真实系统用整数（位运算），并且要一次算几百万个点——末尾给出了 numpy 向量化版本。
空间填充曲线家族与位交织规则
常见的空间填充曲线：
• **Z 序**（Z-order / Morton order）：本实验的主角，构造最简单，只需位交织；
• **Gray 码变体**：把 Z 序的码字再套一次格雷码变换，缓解边界跳跃；
• **Hilbert 曲线**：递归的"U"形拼接，通常被认为局部性更好（**但有争议，见末尾**）；
• **Peano 曲线 / Sierpinski 曲线**：其他递归构造。
位序约定（本 notebook 采用）
设递归层级为 L，x、y 各用 L 个二进制位表示，z 码共 2L 位。
记 x 的各位为 a_{L−1} … a₁ a₀、y 的各位为 b_{L−1} … b₁ b₀（下标 0 是最低位），则
z = ( b₂ a₂ b₁ a₁ b₀ a₀ )₂ （以 L = 3 为例）
即：**y 的一位、x 的一位**交替，从高位排到低位；y 的位落在 z 的奇数位，x 的位落在偶数位。
举例（L = 3，x = 2，y = 1）：
x = 2 = 010 y = 1 = 001
y₂ x₂ y₁ x₁ y₀ x₀
z = 0 0 0 1 1 0 = "000110"
网格规模
| 递归层级 L | 网格 | 单元数 | z 码位宽 |
| --- | --- | --- | --- |
| 1 | 2 × 2 | 4 | 2 |
| 2 | 4 × 4 | 16 | 4 |
| 3 | 8 × 8 | 64 | 6 |
| L | 2^L × 2^L | 2^(2L) | 2L |
量化：浮点坐标必须先变成整数
真实经纬度是浮点数，位交织只认整数，所以要先**量化**到整数网格上。本实验用线性映射：
x = round( (lon + 180) / 360 × (2^L − 1) )
y = round( (lat + 90) / 180 × (2^L − 1) )
量化层级越高越精确，但码越长；量化本身也是一次精度损失——这是所有空间索引共同的取舍。
第 1 步：位交织（坐标 → z 码）
`scalar_to_z_code(recursion_level, x, y)` 做三件事：
1. 把 x、y 转成二进制字符串，并用 `zfill(recursion_level)` 左侧补零到 L 位；
2. 断言两者长度都等于 L（防止坐标超出网格范围）；
3. 用 `"".join(map(add, y_binary, x_binary))` 把 y 与 x **逐位拼接**——
`map(add, y, x)` 把两个字符串的对应字符两两相加，于是第 i 次得到 y[i] + x[i]，
拼起来正好是 `y_{L−1} x_{L−1} … y₀ x₀`。这就是位交织的全部实现。
预期看到：`scalar_to_z_code(3, 2, 1)` 输出 `'000110'`，与上一节的手算一致。
注意这是**字符串**实现，慢且占内存，只用于教学。真实系统用整数位运算（末尾给出）。
from operator import add
def scalar_to_z_code(recursion_level: int, x: int, y: int) -> str:
"""
Converts scalar coordinates into a z-code.
Notice that we use a string-representation of z-codes for educational purposes
a real implementation would represent z-codes in a bit-vector or map them to signed integers
:param recursion_level: the recursion_level of the z-code
:param x: x-coordinate
:param y: y-coordinate
:return: z-code as a string
"""
if recursion_level == 0:
# return empty string for recursion_level 0
return ""
# get the binary representation of x and y as a string and prefix them with zeros:
x_binary = bin(x)[2:].zfill(recursion_level)
y_binary = bin(y)[2:].zfill(recursion_level)
assert len(x_binary) == len(y_binary) == recursion_level
# zip the binary representations of x and y and interleave them#
# this is the actual "trick" of z-codes: the bits of the n=2 binary numbers are interleaved
ret = "".join(map(add, y_binary, x_binary))
assert len(ret) == 2 * recursion_level
return ret
scalar_to_z_code(3, 2, 1)第 2 步：逆变换（z 码 → 坐标）
`z_code_to_scalar(z_code)` 把交织好的码拆回去：
• 码长除以 2 就是递归层级 L；
• 从**最高位**开始，每 2 位一组：`z_code[i]` 属于 y，`z_code[i+1]` 属于 x；
• 每组的权重是 2^exponent，exponent 从 L−1 递减到 0。
所以它能还原，是因为位交织是**严格的一一映射**：z 码里从低位起的第 0、2、4… 位就是 x 的各位，
第 1、3、5… 位就是 y 的各位，按原权重累加回去即可，信息一点没丢。
最后一个 cell 是**往返验证**：`z_code_to_scalar(scalar_to_z_code(2, 3, 1))` 应返回原坐标。
如果这里对不上，99% 是位序写反了（把 x 放到了高位）——这是本实验最常见的错误。
def z_code_to_scalar(z_code: str) -> tuple[int, int, int]:
"""
inverse of scalar_to_z_code function
converts z-code to scalar-coordinates:
:param z_code: string-representation of a z-code
:return: recursion_level, scalar x, scalar y
"""
assert len(z_code) % 2 == 0
dimensions: int = 2
recursion_level: int = len(z_code) // dimensions
x: int = 0
y: int = 0
exponent: int = int(recursion_level) - 1
i: int
for i in range(0, len(z_code), 2):
# reconstruct x and y by reversing the interleaving of the z-code
# basically builds sum of 2^i exponents
y += int(z_code[i]) * 2**exponent
x += int(z_code[i + 1]) * 2**exponent
exponent -= 1
return recursion_level, x, y
z_code_to_scalar("000110")z_code_to_scalar(scalar_to_z_code(2, 3, 1))第 3 步：把整张表的 z 码打印出来
辅助函数 `color_code` 只是给码字分段上 ANSI 颜色（终端转义序列），
让相邻两位的来源在屏幕上更容易分辨；它不影响码值本身。
• `print_z_codes(L)`：按 y 从 0 到 2^L−1、x 从 0 到 2^L−1 的顺序打印整个网格的 z 码；
• `compute_z_codes(L)`：返回 `(z_code, x, y)` 三元组列表，供后面画图用。
观察要点（用 L = 3 试）：
• 沿同一行看，变化的是 z 码的**低位**；
• 跨行时变化的是**高位**；
• 整体读下来，码值从小到大的路径正是那个"Z"字形。
这也直观说明了局部性的来源：**码的高位来自坐标的高位**。
两个点的高位相同 ⇒ 它们落在同一个粗粒度的象限里 ⇒ 它们的 z 码也共享高位
⇒ 按数值排序后必然连成一段连续区间，在一维顺序里仍然邻近。
反过来，一旦跨象限，高位一变，一维顺序上就会出现一次**大跳跃**——
Z 曲线在象限边界上的那些长跳线，就是这个效应的可视化，也是 Z 序的主要缺点。
startEsc: str = "\033[1;"
middleEsc: str = "m"
endEsc: str = "\033[1;m"
def color_code(string: str, stepSize: int = 2):
"""color-codes a string in steps of <stepSize> with different colors
:param string: the string to color-code
:param stepSize: the step-size for coloring
"""
color: int = 31
i: int
for i in range(0, len(string), stepSize):
prefix: str = string[i : i + stepSize]
print(startEsc + str(color) + middleEsc + prefix + endEsc, end="")
color += 1
color_code("hello world")def print_z_codes(recursion_level: int):
"""print all z-codes of a given recursion_level
:param recursion_level: the recursion level
"""
y: int
for y in range(2**recursion_level):
x: int
for x in range(2**recursion_level):
z_code: str = scalar_to_z_code(recursion_level, x, y)
color_code(z_code)
print(" ", end="")
print()
print_z_codes(3)def compute_z_codes(recursion_level: int):
"""computes all z-codes for the given recursion_level
@param recursion_level: the recursion level
@return: a list of tuples containing the z-code, scalar x, and scalar y coordinates
"""
results: list[tuple[str, int, int]] = []
y: int
for y in range(2**recursion_level):
x: int
for x in range(2**recursion_level):
z_code: str = scalar_to_z_code(recursion_level, x, y)
results.append((z_code, x, y))
return results第 4 步：把 Z 曲线画出来
`visualize_z_codes(L)` 用 `ipycanvas` 画交互式画布：它先算出整个网格的 z 码，
**按码值排序**，再依次连线，于是笔尖走出的轨迹就是 Z 曲线。
预期看到（L 从 1 到 5 依次画）：每一级递归都在四个象限里重复一个"Z"字，
层级越高，"Z"字越小越密，整体是一条自相似的折线。
**如果 `ipycanvas` 装不上**，用下面这个 matplotlib 版本替代（只依赖 numpy + matplotlib）：
```python
import numpy as np
import matplotlib.pyplot as plt
def morton(x, y, levels=16):
"""向量化的 Morton 码：y 位放奇数位，x 位放偶数位，与教学版位序一致。"""
x = np.asarray(x, dtype=np.int64)
y = np.asarray(y, dtype=np.int64)
z = np.zeros_like(x)
for i in range(levels):
z |= ((y >> i) & 1) << (2 * i + 1)
z |= ((x >> i) & 1) << (2 * i)
return z
def draw_z_curve(L=3):
n = 1 << L
xs, ys = np.meshgrid(np.arange(n), np.arange(n))
xs, ys = xs.ravel(), ys.ravel()
order = np.argsort(morton(xs.astype(np.int64), ys.astype(np.int64), levels=L),
kind="stable")
plt.rcParams["font.sans-serif"] = ["Microsoft YaHei", "SimHei", "DejaVu Sans"]
fig, ax = plt.subplots(figsize=(6, 6))
ax.plot(xs[order], ys[order], "-o", color="#1F4E79", ms=8, lw=2)
ax.set_aspect("equal")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title(f"Z 曲线（递归层级 {L}，{n}×{n} 网格）")
ax.grid(ls=":", alpha=0.4)
plt.tight_layout()
plt.show()
draw_z_curve(3)
```
把 L 改成 1、2、4 各画一次，看递归结构怎么长出来。
from ipycanvas import Canvas, hold_canvas
from IPython.display import display
def visualize_z_codes(recursion_level: int):
"""visualize the data space in z-code order, i.e. draw the z-curve
:param recursion_level: the recursion level
"""
results: list[tuple[str, int, int]] = compute_z_codes(recursion_level)
results.sort()
x_width: int = 300
y_width: int = x_width
canvas = Canvas(size=(x_width, y_width), sync_image_data=True)
canvas.line_width = 1
max_x: int = 2**recursion_level
max_y: int = max_x
mult_x: float = x_width / max_x
mult_y: float = y_width / max_y
start: tuple[str, int, int] = results.pop(0)
start_x: int = start[1]
start_y: int = start[2]
with hold_canvas(canvas):
while len(results) > 0:
next: tuple[str, int, int] = results.pop(0)
next_x: int = next[1]
next_y: int = next[2]
canvas.begin_path()
canvas.move_to((start_x + 0.5) * mult_x, (start_y + 0.5) * mult_y)
canvas.line_to((next_x + 0.5) * mult_x, (next_y + 0.5) * mult_y)
canvas.stroke()
start_x, start_y = next_x, next_y
return canvas
# visualize the z-codes for recursion levels 1 to 5
for _recursion_level in range(1, 6):
display(visualize_z_codes(_recursion_level))在真实数据上量一量：197 229 个经纬度点
仓库自带 `notebooks/latitude_longitude.bin`：
• 197 229 个点；
• **float32 小端**存储，纬度、经度**交替**存放：`lat = raw[0::2]`、`lon = raw[1::2]`；
• 纬度范围 47.28° – 55.07°，经度范围 −2.01° – 50.51°（欧洲中北部一带）。
用 z 码做范围查询：两阶段流程
① 候选阶段（filter）：把查询矩形翻译成一段（或几段）z 区间，用一维索引（B+ 树 / 有序数组 + 二分）取出区间内所有点
② 精化阶段（refine）：对候选点逐个用真实坐标判断是否在矩形内
第 ① 步取到的是矩形的**超集**：一维区间无法精确表达二维矩形，
所以第 ② 步**不可省略**。工业实现会用 `BIGMIN` / `LITMAX` 把大区间切成若干真正落在矩形内的小段，
跳过大段空洞，从而少取候选。
衡量指标：**候选率 = 候选点数 / 真实命中点数**。越接近 1 越好；
顺序扫描的候选率 = 总点数 / 命中点数。
实测结果（同一 2^16 网格，89 个有效查询窗口，取中位数）
| 方案 | 候选数 / 真实命中数 |
| --- | --- |
| 顺序扫描 | 12.64× |
| Z 序索引 | 3.19× |
| Hilbert 序索引 | **1.77×** |
（口径：随机生成小窗口矩形查询；**只统计真实命中数 ≥ 5 的查询**——命中太少的查询不足以体现差异，
本次有效窗口 89 个；Z 序与 Hilbert 序建立在**同一个** 2^16 网格上。数字以本案例库
《00-实测数字口径表》为准，本机重跑会随随机种子略有差异。）
复现要点：量化到 2^16 网格 → 算码 → 按码排序 → 用 `np.searchsorted` 取区间
（右端用 `side="right"`）→ 与真实命中数求比值 → 取中位数。
一个必须讲的方法学教训：比较两条曲线要站在同一把尺子上
本案例库早期版本做过一次对比：Z 码用 2^16 网格、Hilbert 码用 2^10 网格。
结果 Hilbert 看起来领先很多——但那个"优势"是假的。
原因：网格越粗，同一个查询框在码空间里覆盖的区间越短，候选集天然就更少。
也就是说，测出来的差距有相当一部分来自**量化层级不同**，而不是曲线本身的性质。
那个数字已经作废，**不要再引用**。
现在把两者统一到 2^16 网格后，Hilbert 相对 Z 序的领先幅度从约 2.5 倍回落到约 1.8 倍（3.19 → 1.77）。
这才是可复现、可辩护的对比。
结论值得单独记住：**任何对比实验，只要两方的"其他条件"不一致，读出来的差距就不可信。**
比较两条空间填充曲线之前，先确认它们站在同一把尺子上——同网格、同查询集、同过滤阈值、同统计口径。
顺便说一句：别把"Hilbert 一定比 Z 序好"当定论
这个说法**有争议**，取决于用什么度量。Dai & Su（DMTCS 2003）在"平均跨簇距离"与
"平均总跨簇距离"两项统计上得到的结论是 **Z 序家族优于 Hilbert 家族**，与很多教材的说法相反。
本实验在"小窗口范围查询的候选率"这个口径下 Hilbert 更好（1.77× vs 3.19×），
但**换一个度量就可能翻转**。工程上的正确做法是：用自己的真实数据和真实查询分布跑一遍，
按候选率（或干脆按端到端查询耗时）来选，而不是照搬结论。
Z 序在真实系统里的落地形态
| 系统 / 结构 | 做法 | 备注 |
| --- | --- | --- |
| **Geohash** | 基于 Z 序，把经纬度交替二分后编成 base32 字符串 | 前缀相同表示区域相邻，常用作前缀匹配的粗过滤 |
| **Google S2** | 用**六条 Hilbert 曲线**首尾相连覆盖球面（立方体投影的六个面） | 层级 level 0–30，单元 ID 是 64 位整数 |
| **Uber H3** | 六边形层次网格，分辨率 0–15 | 总单元数 c(r) = 2 + 120·7^r；**每层恰有 12 个五边形**（六边形无法完美铺满球面，这是拓扑上不可避免的） |
| **R 树** | 用最小外接矩形（MBR）组织空间对象，不是线性化路线 | 出处：Guttman，《R-trees: A Dynamic Index Structure for Spatial Searching》，SIGMOD 1984 |
| **PostGIS** | GiST / SP-GiST 索引 | 见下方 `ST_DWithin` 的提醒 |
| **GeoParquet** | 列式空间数据文件格式，带 bbox 统计便于裁剪 | 配合 DuckDB / Sedona 做列存空间分析 |
| **Apache Sedona** | 在 Spark 之上提供空间分区、空间索引与空间连接算子 | 海量空间数据的分布式处理路径 |
三者取舍：Geohash / S2 的四边形单元便于位运算与前缀查询；H3 的六边形单元到邻居的距离更均匀
（每个单元有 6 个等距邻居），适合邻域聚合，代价是那 12 个五边形。
**PostGIS 的一个高频坑**：`ST_DWithin` 是**索引感知**（index-aware）的距离阈值谓词，
能用 GiST 索引先做包围盒粗过滤；而 `ST_Distance(...) < d` 只是一个普通表达式，
优化器不知道它等价于"距离阈值"，只能逐行算距离再比较，**不走索引**，大表上退化成全表扫描。
（依据：PostGIS 官方文档 Chapter 5。）
```sql
-- 推荐：能走 GiST 索引
SELECT * FROM restaurants
WHERE ST_DWithin(geom, ST_MakePoint(:lon, :lat)::geography, 500);
-- 不推荐：全表扫描
SELECT * FROM restaurants
WHERE ST_Distance(geom, ST_MakePoint(:lon, :lat)::geography) < 500;
```
思考题与详解
**1. 手算 `scalar_to_z_code(3, 5, 2)`：写出 x、y 的二进制与交织过程，再用代码验证。**
<details><summary>答案要点</summary>
x = 5 = 101 y = 2 = 010
位序 y₂ x₂ y₁ x₁ y₀ x₀
0 1 1 0 0 1
z 码 = "011001"
即 a₂a₁a₀ = 101（x = 5），b₂b₁b₀ = 010（y = 2），按 y 位在前、x 位在后交织得 `0 1 1 0 0 1`。
验证：`print(scalar_to_z_code(3, 5, 2))` 输出 `'011001'`。
</details>
**2. 求 `z_code_to_scalar("011001")` 的结果，并说明逆变换为什么能还原出 x 与 y。**
<details><summary>答案要点</summary>
码长 6 ⇒ L = 3。按 2 位一组拆：`01 | 10 | 01`，权重从 2² 递减到 2⁰：
第 1 组 "01"（权重 4）：y += 0×4 = 0，x += 1×4 = 4
第 2 组 "10"（权重 2）：y += 1×2 = 2，x += 0×2 = 0
第 3 组 "01"（权重 1）：y += 0×1 = 0，x += 1×1 = 1
结果：x = 5，y = 2 → 返回 (3, 5, 2)
能还原的原因是位交织是**一一映射**：从 z 码最低位起，第 0、2、4… 位恰好是 x 的各位，
第 1、3、5… 位是 y 的各位，把它们按原来的权重 2^i 累加回去即可，信息没有丢失。
（若你的结果对不上，先检查位序：本 notebook 约定 **y 在奇数位、x 在偶数位**。）
</details>
**3. 递归层级 L = 16 时，网格边长是多少？单元总数是多少？z 码需要几个二进制位、能否放进一个 64 位整数？**
<details><summary>答案要点</summary>
• 网格边长 = 2^16 = **65 536**
• 单元总数 = 2^16 × 2^16 = 2^32 = **4 294 967 296**（约 43 亿）
• z 码位宽 = 2L = **32 位**，**能**放进 64 位整数（也能放进无符号 32 位整数）
延伸：若 L 取到 32，z 码就是 64 位，此时必须注意有符号整数溢出，要用 `uint64`；
Python 的 `int` 虽然不会溢出，但与 C/NumPy 交互时类型要显式指定。
</details>
**4. 为什么范围查询必须"先取候选、再精确过滤"？举一个具体反例：某点的 z 值落在 `[z(x₀,y₀), z(x₁,y₁)]` 内，但点本身不在查询矩形里。**
<details><summary>答案要点</summary>
取 L = 2（4×4 网格），查询矩形 x ∈ [1, 2]、y ∈ [0, 1]。
z(1, 0)：x = 01, y = 00 → y₁ x₁ y₀ x₀ = 0 0 0 1 = 1
z(2, 1)：x = 10, y = 01 → y₁ x₁ y₀ x₀ = 0 1 1 0 = 6
候选区间 = [1, 6]
区间里的 z = 2 对应二进制 `0010`：`y₁ x₁ y₀ x₀ = 0 0 1 0` ⇒ y = `01` = 1，x = `00` = 0，
即点 **(0, 1)**。它的 x = 0 不在 [1, 2] 内——**落在 z 区间里却不在矩形里**。
原因：**一维区间无法精确表达二维矩形**，z 区间只是矩形的一个超集。所以精化阶段
（逐个用真实坐标判断）不可省略。工业实现用 `BIGMIN` 算法把大区间切成若干"真正落在矩形内"的小段，
跳过大段空洞，从而减少这类无用候选。
</details>
**5. 关于 Hilbert 与 Z 序的优劣，文献结论一致吗？Dai & Su（DMTCS 2003）得到的是什么结论？工程上应该怎么处理？**
<details><summary>答案要点</summary>
**不一致。** Dai & Su（DMTCS 2003）在"平均跨簇距离"与"平均总跨簇距离"两项统计上，
得到的结论是 **Z 序家族优于 Hilbert 家族**；而很多教材与工程经验又说 Hilbert 局部性更好。
分歧来自**度量口径**不同：用哪种距离、按哪种查询分布来统计，结论会变。
本实验在"小窗口范围查询的候选率"这个口径下 Hilbert 更好（1.77× vs 3.19×），
但**换一个度量就可能翻转**。
工程做法：**不要照搬结论，用自己的真实数据 + 真实查询分布实测一遍**，
按候选率（或干脆按端到端查询耗时）来选曲线与层级。
</details>
**6. 有人把 Z 码用 2^16 网格、Hilbert 码用 2^10 网格做对比，得出 Hilbert 大幅领先的结论。这个比较公平吗？请设计一个公平的方案。**
<details><summary>答案要点</summary>
**不公平。** 网格越粗，同一个查询框在码空间里覆盖的区间越短，候选集天然更少。
所以那个"领先"里混进了量化层级的差异，不是曲线本身的性质——
本案例库早期那个数字已经作废，不要引用。统一到同一 2^16 网格后，
Hilbert 的领先幅度从约 2.5 倍回落到约 1.8 倍（3.19 → 1.77）。
公平方案：
1. **同一量化层级**：两者都用 L = 16（或都用 L = 10），量化公式完全一致；
2. **同一查询集**：同一批随机小窗口，同样的"命中数 ≥ 5 才计入"过滤；
3. **同一统计口径**：同样的候选数计算方式（`searchsorted` 的 `side` 一致）、同样的中位数而非均值；
4. 若 Hilbert 的逐点循环太慢，可**抽样**再做，但两侧必须抽**同一批点**。
一句话：**比较两条曲线之前，先确认它们站在同一把尺子上。**
</details>
**7. Z 序在真实系统里有哪些落地形态？请说明 Geohash、Google S2、Uber H3 各自的编码方式，以及 H3 每层有多少个五边形。**
<details><summary>答案要点</summary>
• **Geohash**：基于 **Z 序**，把经纬度交替二分后编成 base32 字符串；
前缀相同表示区域相邻，常用作前缀匹配的粗过滤。
• **Google S2**：把球面（立方体投影后的六个面）用 **六条 Hilbert 曲线首尾相连**覆盖，
层级 level 0–30，单元 ID 是 64 位整数。
• **Uber H3**：**六边形**层次网格，分辨率 0–15，总单元数 c(r) = 2 + 120·7^r；
由于六边形无法完美铺满球面，**每层恰有 12 个五边形**（拓扑上不可避免）。
取舍：Geohash / S2 的四边形单元便于位运算与前缀查询；H3 的六边形单元到邻居的距离更均匀
（每个单元 6 个等距邻居），适合邻域聚合，代价是那 12 个五边形。
另外，不走线性化路线的还有 **R 树**（Guttman，SIGMOD 1984，用最小外接矩形组织空间对象），
以及 PostGIS 的 GiST / SP-GiST 索引、GeoParquet 列式格式、Apache Sedona（Spark 上的空间计算）。
</details>
**8. 在 PostGIS 里想查"距离某点 500 米以内的餐厅"，为什么应该写 `ST_DWithin` 而不是 `ST_Distance < 500`？**
<details><summary>答案要点</summary>
`ST_DWithin(geom1, geom2, d)` 是**索引感知**（index-aware）的谓词：
PostGIS 会先用 GiST 索引的包围盒做粗过滤，再算精确距离，因此**能走索引**。
而 `ST_Distance(geom1, geom2) < d` 只是一个普通表达式：优化器不知道它等价于"距离阈值"，
只能对每一行都算一次距离再比较，**不使用索引**，在大表上退化成全表扫描。
```sql
-- 推荐：能走 GiST 索引
SELECT * FROM restaurants
WHERE ST_DWithin(geom, ST_MakePoint(:lon, :lat)::geography, 500);
-- 不推荐：全表扫描
SELECT * FROM restaurants
WHERE ST_Distance(geom, ST_MakePoint(:lon, :lat)::geography) < 500;
```
（依据：PostGIS 官方文档 Chapter 5。）
这个坑与 Z 序实验是同一个道理：**能不能用上索引，取决于谓词的形式，而不是语义。**
想让优化器用索引，就要写成它能识别的谓词形式。
</details>