如何递归地查找两个表之间的相交地理
Posted
技术标签:
【中文标题】如何递归地查找两个表之间的相交地理【英文标题】：How to find intersecting geographies between two tables recursively 【发布时间】：2017-06-13 11:09:49 【问题描述】：我正在运行 Postgres 9.6.1 和 PostGIS 2.3.0 r15146 并且有两个表。geographies 可能有 150,000,000 行，paths 可能有 10,000,000 行：
CREATE TABLE paths (id uuid NOT NULL, path path NOT NULL, PRIMARY KEY (id))
CREATE TABLE geographies (id uuid NOT NULL, geography geography NOT NULL, PRIMARY KEY (id))
给定表geographies 的ids 数组/集，查找所有相交路径和几何图形的“最佳”方法是什么？
换句话说，如果一个初始的geography 有一个对应的相交path，我们还需要找到这个path 相交的所有其他 geographies。从那里，我们需要找到这些新发现的geographies 相交的所有其他paths，依此类推，直到我们找到所有可能的相交。
初始地理 ID（我们的输入）可能在 0 到 700 之间。平均约为 40。 最小交叉点为 0，最大约为 1000 个。平均可能在 20 个左右，通常少于 100 个连接。
我已经创建了一个执行此操作的函数，但我对 PostGIS 中的 GIS 和一般的 Postgres 是新手。我已经发布了my solution as an answer to this question。
我觉得应该有比我想出的更雄辩和更快的方法来做到这一点。
【问题讨论】：
您是否考虑将pathIds 和geographyIds 存储在表中而不是数组中？
当然。我现在有一些工作可以返回单个table(id uuid, type character varying)，其中type 是path 或geography。
考虑在Geographic Information Systems 询问有知识渊博的 PostGIS 人员。
感谢@TobySpeight。我已经请版主提出问题。
我添加了自己的解决方案作为答案。
【参考方案1】：
Your function 可以彻底简化。
设置
我建议您将列paths.path 转换为数据类型geography（或至少geometry）。 path is a native Postgres type 与 PostGIS 功能和空间索引不兼容。您必须使用 path::geometry 或 path::geometry::geography (resulting in a LINESTRING internally) 才能使其与 ST_Intersects() 等 PostGIS 函数一起使用。
我的答案基于这些改编的表格：
CREATE TABLE paths (
id uuid PRIMARY KEY
, path geography NOT NULL
);
CREATE TABLE geographies (
id uuid PRIMARY KEY
, geography geography NOT NULL
, fk_id text NOT NULL
);所有内容都适用于数据类型geometry 的两列。 geography 通常更准确但也更昂贵。使用哪个？ Read the PostGIS FAQ here.
解决方案 1：优化您的功能
CREATE OR REPLACE FUNCTION public.function_name(_fk_ids text[])
RETURNS TABLE(id uuid, type text)
LANGUAGE plpgsql AS
$func$
DECLARE
_row_ct int;
_loop_ct int := 0;
BEGIN
CREATE TEMP TABLE _geo ON COMMIT DROP AS -- dropped at end of transaction
SELECT DISTINCT ON (g.id) g.id, g.geography, _loop_ct AS loop_ct -- dupes possible?
FROM geographies g
WHERE g.fk_id = ANY(_fk_ids);
GET DIAGNOSTICS _row_ct = ROW_COUNT;
IF _row_ct = 0 THEN -- no rows found, return empty result immediately
RETURN; -- exit function
END IF;
CREATE TEMP TABLE _path ON COMMIT DROP AS
SELECT DISTINCT ON (p.id) p.id, p.path, _loop_ct AS loop_ct
FROM _geo g
JOIN paths p ON ST_Intersects(g.geography, p.path); -- no dupes yet
GET DIAGNOSTICS _row_ct = ROW_COUNT;
IF _row_ct = 0 THEN -- no rows found, return _geo immediately
RETURN QUERY SELECT g.id, text 'geo' FROM _geo g;
RETURN;
END IF;
ALTER TABLE _geo ADD CONSTRAINT g_uni UNIQUE (id); -- required for UPSERT
ALTER TABLE _path ADD CONSTRAINT p_uni UNIQUE (id);
LOOP
_loop_ct := _loop_ct + 1;
INSERT INTO _geo(id, geography, loop_ct)
SELECT DISTINCT ON (g.id) g.id, g.geography, _loop_ct
FROM _paths p
JOIN geographies g ON ST_Intersects(g.geography, p.path)
WHERE p.loop_ct = _loop_ct - 1 -- only use last round!
ON CONFLICT ON CONSTRAINT g_uni DO NOTHING; -- eliminate new dupes
EXIT WHEN NOT FOUND;
INSERT INTO _path(id, path, loop_ct)
SELECT DISTINCT ON (p.id) p.id, p.path, _loop_ct
FROM _geo g
JOIN paths p ON ST_Intersects(g.geography, p.path)
WHERE g.loop_ct = _loop_ct - 1
ON CONFLICT ON CONSTRAINT p_uni DO NOTHING;
EXIT WHEN NOT FOUND;
END LOOP;
RETURN QUERY
SELECT g.id, text 'geo' FROM _geo g
UNION ALL
SELECT p.id, text 'path' FROM _path p;
END
$func$;
呼叫：
SELECT * FROM public.function_name('foo,bar');
比你现有的快得多。
要点
您基于整个集合的查询，而不是仅针对集合的最新添加。在不需要的情况下，每个循环都会变得越来越慢。我添加了一个循环计数器 (loop_ct) 来避免多余的工作。
确保在geographies.geography 和paths.path 上有空间GiST 索引：
CREATE INDEX geo_geo_gix ON geographies USING GIST (geography);
CREATE INDEX paths_path_gix ON paths USING GIST (path);
由于 Postgres 9.5 index-only scans 将成为 GiST 索引的一个选项。您可以添加 id 作为第二个索引列。好处取决于许多因素，您必须进行测试。 但是，uuid 类型没有合适的运算符 GiST 类。安装扩展程序 btree_gist 后，它将与 bigint 一起使用：
Postgres multi-column index (integer, boolean, and array)
Multicolumn index on 3 fields with heterogenous data types
在g.fk_id 上也有一个合适的索引。同样，(fk_id, id, geography) 上的多列索引如果您可以从中获得仅索引扫描，则可能会有所回报。默认 btree 索引，fk_id 必须是第一个索引列。特别是如果您经常运行查询并且很少更新表并且表行比索引宽得多。
您可以在声明时初始化变量。重写后只需要一次。
ON COMMIT DROP 在事务结束时自动删除临时表。所以我明确删除了删除表。但是，如果您在 same 事务中调用该函数两次，则会出现异常。在函数中，我将检查临时表是否存在并在这种情况下使用TRUNCATE。相关：
How to check if a table exists in a given schema
使用GET DIAGNOSTICS 获取行数，而不是针对该计数运行另一个查询。
Count rows affected by DELETE
您需要GET DIAGNOSTICS。 CREATE TABLE 不设置FOUND（如手册中所述）。
在填表之后添加索引或 PK / UNIQUE 约束会更快。而不是在我们真正需要它之前。
ON CONFLICT ... DO ... 是自 Postgres 9.5 以来更简单、更便宜的 UPSERT 方式。
How to UPSERT (MERGE, INSERT ... ON DUPLICATE UPDATE) in PostgreSQL?
对于命令的简单形式，您只需列出索引列或表达式（如ON CONFLICT (id) DO ...）并让 Postgres 执行唯一索引推断以确定仲裁约束或索引。后来我通过直接提供约束进行了优化。但是为此，我们需要一个实际的 约束 - 唯一索引是不够的。相应地修复。 Details in the manual here.
它可能有助于手动ANALYZE 临时表以帮助 Postgres 找到最佳查询计划。 （但我认为您的情况不需要它。）
Are regular VACUUM ANALYZE still recommended under 9.1?
_geo_ct - _geographyLength > 0 是 _geo_ct > _geographyLength 的一种尴尬且更昂贵的说法。但现在完全没有了。
不要引用语言名称。只需LANGUAGE plpgsql。
您的函数参数是varchar[]，用于fk_id 数组，但您后来评论说：
它是一个代表地理区域的
bigint字段（它实际上是一个预先计算的s2cellid 在级别15）。
我不知道 s2cell id 在第 15 级，但理想情况下，您传递一个匹配数据类型的数组，或者如果这不是默认为 text[] 的选项。
还因为您发表了评论：
总是正好有 13 个
fk_ids 传入。
这似乎是 VARIADIC 函数参数的完美用例。所以你的函数定义是：
CREATE OR REPLACE FUNCTION public.function_name(_fk_ids <b>VARIADIC</b> text[]) ...详情：
Pass multiple values in single parameter解决方案 2：使用递归 CTE 的普通 SQL
很难将rCTE 包裹在两个交替循环中，但可以通过一些 SQL 技巧来实现：
WITH RECURSIVE cte AS (
SELECT g.id, g.geography::text, NULL::text AS path, text 'geo' AS type
FROM geographies g
WHERE g.fk_id = ANY($kf_ids) -- your input array here
UNION
SELECT p.id, g.geography::text, p.path::text
, CASE WHEN p.path IS NULL THEN 'geo' ELSE 'path' END AS type
FROM cte c
LEFT JOIN paths p ON c.type = 'geo'
AND ST_Intersects(c.geography::geography, p.path)
LEFT JOIN geographies g ON c.type = 'path'
AND ST_Intersects(g.geography, c.path::geography)
WHERE (p.path IS NOT NULL OR g.geography IS NOT NULL)
)
SELECT id, type FROM cte;
就是这样。 您需要与上述相同的索引。您可以将其包装到 SQL 函数中以供重复使用。
主要补充点
强制转换为text 是必要的，因为geography 类型不是“可散列的”（geometry 也是如此）。 (See this open PostGIS issue for details.) 通过转换为 text 来解决它。仅凭(id, type)，行是唯一的，为此我们可以忽略geography 列。回滚到geography 以加入。不应该花费太多额外费用。
我们需要两个LEFT JOIN，所以不要排除行，因为在每次迭代中，两个表中只有一个可能会贡献更多行。
最后的条件确保我们还没有完成：
WHERE (p.path IS NOT NULL OR g.geography IS NOT NULL)
这是有效的，因为重复的结果被排除在临时 中间表。 The manual:
对于
UNION（但不是UNION ALL），丢弃重复的行和 复制任何先前的结果行。包括所有剩余的行 递归查询的结果，并将它们放在一个临时的 中间表。
那么哪个更快？
rCTE 可能比用于小型结果集的函数更快。函数中的临时表和索引意味着更多的开销。不过，对于大型结果集，函数可能会更快。只有使用您的实际设置进行测试才能给您一个明确的答案。*
见the OP's feedback in the comment。【讨论】：
@DavidMurdoch：我很感兴趣：您是否碰巧测试了这两种解决方案并比较了小/大（多次迭代）结果集的性能？ （感谢慷慨的赏金，顺便说一句！） 我做到了。我选择 rCTE，因为它更快、更容易管理。我最终返回的列不仅仅是 id 和几何图形（geometry 看起来比geography 快得多）。奇怪的是，我发现对于从CASE 中选择的某些列，它会将执行时间增加 10 倍！我现在不记得具体细节了，但我觉得很奇怪。
性能测试可能很棘手，有很多可能的影响因素。可能CASE 使 Postgres 切换到不同的查询计划。可能是服务器配置问题 (cost settings, table statistics?) 因子 10 听起来很奇怪。查询计划通常有助于理解。如果您需要答案，可能是另一个问题的材料。
我再次尝试复制但无法弄清楚。我在办公室展示了其他一些开发人员，他们也很困惑。如果我弄清楚了，我会尝试制作一个简化的测试用例并发布一个关于它的问题。【参考方案2】：
我认为最好在此处发布我自己的解决方案，即使它不是最佳的。
这是我想出的（使用史蒂夫钱伯斯的建议）：
CREATE OR REPLACE FUNCTION public.function_name(
_fk_ids character varying[])
RETURNS TABLE(id uuid, type character varying)
LANGUAGE 'plpgsql'
COST 100.0
VOLATILE
ROWS 1000.0
AS $function$
DECLARE
_pathLength bigint;
_geographyLength bigint;
_currentPathLength bigint;
_currentGeographyLength bigint;
BEGIN
DROP TABLE IF EXISTS _pathIds;
DROP TABLE IF EXISTS _geographyIds;
CREATE TEMPORARY TABLE _pathIds (id UUID PRIMARY KEY);
CREATE TEMPORARY TABLE _geographyIds (id UUID PRIMARY KEY);
-- get all geographies in the specified _fk_ids
INSERT INTO _geographyIds
SELECT g.id
FROM geographies g
WHERE g.fk_id= ANY(_fk_ids);
_pathLength := 0;
_geographyLength := 0;
_currentPathLength := 0;
_currentGeographyLength := (SELECT COUNT(_geographyIds.id) FROM _geographyIds);
-- _pathIds := ARRAY[]::uuid[];
WHILE (_currentPathLength - _pathLength > 0) OR (_currentGeographyLength - _geographyLength > 0) LOOP
_pathLength := (SELECT COUNT(_pathIds.id) FROM _pathIds);
_geographyLength := (SELECT COUNT(_geographyIds.id) FROM _geographyIds);
-- gets all paths that have paths that intersect the geographies that aren't in the current list of path ids
INSERT INTO _pathIds
SELECT DISTINCT p.id
FROM paths p
JOIN geographies g ON ST_Intersects(g.geography, p.path)
WHERE
g.id IN (SELECT _geographyIds.id FROM _geographyIds) AND
p.id NOT IN (SELECT _pathIds.id from _pathIds);
-- gets all geographies that have paths that intersect the paths that aren't in the current list of geography ids
INSERT INTO _geographyIds
SELECT DISTINCT g.id
FROM geographies g
JOIN paths p ON ST_Intersects(g.geography, p.path)
WHERE
p.id IN (SELECT _pathIds.id FROM _pathIds) AND
g.id NOT IN (SELECT _geographyIds.id FROM _geographyIds);
_currentPathLength := (SELECT COUNT(_pathIds.id) FROM _pathIds);
_currentGeographyLength := (SELECT COUNT(_geographyIds.id) FROM _geographyIds);
END LOOP;
RETURN QUERY
SELECT _geographyIds.id, 'geography' AS type FROM _geographyIds
UNION ALL
SELECT _pathIds.id, 'path' AS type FROM _pathIds;
END;
$function$;
【讨论】：
我对您的回答中的ST_Intersects(g.geography, p.path) 感到困惑。在 PostGIS 2.3.1 中没有ST_Intersects() 的变体可以接受数据类型path，并且path 的唯一注册转换目标是point、polygon 和geometry。 paths.path真的是数据类型path吗？
你完全正确。 p.path 实际上是我数据库中的 geometry 类型！对于那个很抱歉！将类型更改为geography 没有问题，因此您的答案看起来会是赢家。对于此函数，是否可以安全地假设使用 geography 类型会比 geometry 更高效，因为我会避免强制转换？
演员阵容很便宜。在geography 类型上操作通常更精确但也更昂贵。我建议对两个 geometry 或两个 geography 值进行所有测试。理想情况下，您有匹配的类型。如果基础表具有不同的类型，则在查询中显式转换为一种类型 (geography::geometry)，并为相应的列提供匹配的表达式索引，例如 CREATE INDEX ON geographies USING GIST ((geography::geometry)); 注意额外的括号。选哪个？阅读：postgis.net/docs/PostGIS_FAQ.html#idm1018【参考方案3】：
Sample plot and data from this script
它可以是与聚合函数的纯关系。此实现使用一张path 表和一张point 表。两者都是几何。这一点比通用地理更容易创建测试数据并进行测试，但它应该易于适应。
create table path (
path_text text primary key,
path geometry(linestring) not null
);
create table point (
point_text text primary key,
point geometry(point) not null
);
一种保持聚合函数状态的类型：
create type mpath_mpoint as (
mpath geometry(multilinestring),
mpoint geometry(multipoint)
);
状态构建函数：
create or replace function path_point_intersect (
_i mpath_mpoint[], _e mpath_mpoint
) returns mpath_mpoint[] as $$
with e as (select (e).mpath, (e).mpoint from (values (_e)) e (e)),
i as (select mpath, mpoint from unnest(_i) i (mpath, mpoint))
select array_agg((mpath, mpoint)::mpath_mpoint)
from (
select
st_multi(st_union(i.mpoint, e.mpoint)) as mpoint,
(
select st_collect(gd)
from (
select gd from st_dump(i.mpath) a (a, gd)
union all
select gd from st_dump(e.mpath) b (a, gd)
) s
) as mpath
from i inner join e on st_intersects(i.mpoint, e.mpoint)
union all
select i.mpoint, i.mpath
from i inner join e on not st_intersects(i.mpoint, e.mpoint)
union all
select e.mpoint, e.mpath
from e
where not exists (
select 1 from i
where st_intersects(i.mpoint, e.mpoint)
)
) s;
$$ language sql;
聚合：
create aggregate path_point_agg (mpath_mpoint) (
sfunc = path_point_intersect,
stype = mpath_mpoint[]
);
此查询将返回一组包含匹配路径/点的multilinestring, multipoint 字符串：
select st_astext(mpath), st_astext(mpoint)
from unnest((
select path_point_agg((st_multi(path), st_multi(mpoint))::mpath_mpoint)
from (
select path, st_union(point) as mpoint
from
path
inner join
point on st_intersects(path, point)
group by path
) s
)) m (mpath, mpoint)
;
st_astext | st_astext
-----------------------------------------------------------+-----------------------------
MULTILINESTRING((-10 0,10 0,8 3),(0 -10,0 10),(2 1,4 -1)) | MULTIPOINT(0 0,0 5,3 0,5 0)
MULTILINESTRING((-9 -8,4 -8),(-8 -9,-8 6)) | MULTIPOINT(-8 -8,2 -8)
MULTILINESTRING((-7 -4,-3 4,-5 6)) | MULTIPOINT(-6 -2)
【讨论】：
以上是关于如何递归地查找两个表之间的相交地理的主要内容，如果未能解决你的问题，请参考以下文章