PostGIS 有多种地理哈希函数。如果你的字符串足够长,搜索会变得更快(每个盒子+它的 8 个邻居的碰撞更少),但插入新点时 geohash 的生成速度会变慢。
问题还在于您想要多准确。随着纬度的增加,纬度/经度“距离”会恶化,因为经度从赤道的约 110 公里缩小到两极的 0,而纬度始终约为 110 公里。在中纬度 45 度处,经度接近 79 公里,距离误差为 2 倍 (sqr(110/79))。为您提供纬度/经度对之间的真实距离的球面距离的计算成本非常昂贵(需要进行大量三角学计算),然后您的地理哈希将无法工作(除非您将所有点转换为平面坐标)。
可能有效的解决方案如下:
-
CREATE INDEX hash8 ON tablename(substring(hash_column FROM 1 FOR 8))
。这为您提供了两倍于分辨率的框的索引,这有助于查找点并减少搜索相邻散列框的需要。
- On
INSERT
对于一个点,使用 PostGIS 将其长度为 9(大约 10m 分辨率)的 geohash 计算到 hash_column 中。你可以使用BEFORE INSERT TRIGGER
here.
在一个函数中:
- 给定一个点,通过查找所有 geohash 值缩短为 8 个字符(等于给定点 8 字符 geohash)的点(因此是上面的索引)来找到最近的点。
- 使用球坐标计算到每个遇到的点的距离,保留最近的点。但由于您只是寻找最近的点(至少最初是这样),因此不要使用球坐标搜索距离,而是使用下面的优化,这应该会使搜索速度更快。
- 计算给定点是否比最近的计算点更接近由 8 字符 geohash 确定的框的边缘。如果是这样,请在其 8 个邻居的所有点上使用 7 字符 geohash 重复该过程。这可以通过计算到各个框边和角的距离并仅评估相关的邻居散列框来高度优化;我把这个留给你去修改。
无论如何,这都不会特别快。如果您确实要处理数十亿个点,您可能需要考虑聚类,它对地理哈希有一个相当“自然”的解决方案(将您的表分解为substring(hash_column FROM 1 FOR 2)
例如,给你四个象限)。只要确保您考虑到跨境搜索即可。
两项优化可以相当快地制作:
First,“标准化”您的球面坐标(意思是:随着纬度的增加补偿经度长度的减少),以便您可以使用“伪笛卡尔”方法搜索最近的点。这只适用于点距离很近的情况,但由于您正在处理很多点,这应该不是问题。更具体地说,这应该适用于长度为 6 或更长的 geohash 框中的所有点。
假设 WGS84 椭球体(用于所有 GPS 设备),地球长轴 (a) 为 6,378,137 米,椭圆率 (e2) 为 0.00669438。一秒经度的长度为
longSec := Pi * a * cos(lat) / sqrt(1 - e2 * sqr(sin(lat))) / 180 / 3600
or
longSec := 30.92208078 * cos(lat) / sqrt(1 - 0.00669438 * sqr(sin(lat)))
对于一秒的纬度:
latSec := 30.870265 - 155.506 * cos(2 * lat) + 0.0003264 + cos(4 * lat)
使本地坐标系成为“正方形”的校正因子是将经度值乘以longSec/latSec
.
Secondly,由于您正在寻找最近的点,因此不要搜索距离,因为平方根的计算成本很高。相反,如果愿意的话,可以搜索平方根(平方距离)内的术语,因为它具有选择最近点的相同属性。
在伪代码中:
CREATE FUNCTION nearest_point(pt geometry, ptHash8 char(8)) RETURNS integer AS $$
DECLARE
corrFactor double precision;
ptLat double precision;
ptLong double precision;
currPt record;
minDist double precision;
diffLat double precision;
diffLong double precision;
minId integer;
BEGIN
minDist := 100000000.; -- a large value, 10km (squared)
ptLat := ST_Y(pt);
ptLong := ST_X(pt);
corrFactor := 30.92208078 * cos(radians(ptLat)) / (sqrt(1 - 0.00669438 * power(sin(radians(ptLat)), 2)) *
(30.870265 - 155.506 * cos(2 * radians(ptLat)) + 0.0003264 + cos(4 * radians(ptLat))));
FOR currPt IN SELECT * FROM all_points WHERE hash8 = ptHash8
LOOP
diffLat := ST_Y(currPt.pt) - ptLat;
diffLong := (ST_X(currPt.pt) - ptLong) * corrFactor; -- "square" things out
IF (diffLat * diffLat) < (minDist * diffLong * diffLong) THEN -- no divisions here to speed thing up a little further
minDist := (diffLat * diffLat) / (diffLong * diffLong); -- this does not happen so often
minId := currPt.id;
END IF;
END LOOP;
IF minDist < 100000000. THEN
RETURN minId;
ELSE
RETURN NULL;
END IF;
END; $$ LANGUAGE PLPGSQL STRICT;
不用说,这在 C 语言函数中会快得多。另外,不要忘记进行边界检查以查看是否需要搜索相邻的 geohash 框。
顺便说一句,“空间纯粹主义者”不会在 8 字符 geohash 上建立索引并从那里进行搜索;相反,它们将从 9 个字符的哈希开始,并从那里向外工作。然而,初始哈希盒中的“未命中”(因为没有其他点或者您靠近哈希盒一侧)的成本很高,因为您必须开始计算到相邻哈希盒的距离并提取更多数据。在实践中,您应该使用大约是典型最近点大小两倍的散列框;该距离是多少取决于您的点集。