Эффективный способ расчета площади перекрытия в PostGIS

#postgresql #postgis #psql #overlap

Вопрос:

Я работаю над вычислением перекрытия двух таблиц/слоев в базе данных PostGIS. Один набор представляет собой сетку шестиугольников, для которой я хочу рассчитать долю перекрытия с другим набором полигонов для каждого из шестиугольников. В наборе мультиполигонов также есть несколько перекрывающихся полигонов, поэтому мне нужно их растворить/объединить. До этого я делал это в FME, но для некоторых больших полигонов постоянно не хватало памяти. Вот почему я хочу сделать это в базе данных (и PostGIS должен быть очень способен на это).

Вот что у меня есть до сих пор, и это работает, и память сейчас не на исходе:

 EXPLAIN ANALYZE
WITH rh_union AS (
    SELECT (ST_Dump(ST_Union(rh.geometry))).geom AS geometry
    FROM relevant_habitats rh 
    WHERE rh.assessment_area_id=1
) 
SELECT h.receptor_id, 
SUM(ROUND((ST_Area(ST_Intersection(rhu.geometry, h.geometry)) / ST_Area(h.geometry))::numeric,3)) AS frac_overlap
FROM rh_union rhu, hexagons h
WHERE h.zoom_level=1 AND ST_Intersects(rhu.geometry, h.geometry)
GROUP BY h.receptor_id
 

Поэтому я сначала разбираю мультиполигон и объединяю все, что могу. Затем рассчитайте наложение шестиугольников на многоугольники. Затем вычислите (сумму всех небольших кусочков) площади.

Теперь мой вопрос заключается в следующем:

  • является ли это эффективным способом сделать это? Или был бы лучший способ?

(И дополнительный вопрос: правильно ли сначала использовать ST_Union, а затем ST_Dump?)

—Обновление с ОБЪЯСНЕНИЕМ РЕЗУЛЬТАТОВ АНАЛИЗА для одной области:

 "QUERY PLAN"
"GroupAggregate  (cost=1996736.77..15410052.20 rows=390140 width=36) (actual time=571.303..1137.657 rows=685 loops=1)"
"  Group Key: h.receptor_id"
"  ->  Sort  (cost=1996736.77..1998063.55 rows=530712 width=188) (actual time=571.090..620.379 rows=806 loops=1)"
"        Sort Key: h.receptor_id"
"        Sort Method: external merge  Disk: 42152kB"
"        ->  Nested Loop  (cost=55.53..1848314.51 rows=530712 width=188) (actual time=382.769..424.643 rows=806 loops=1)"
"              ->  Result  (cost=55.25..1321.51 rows=1000 width=32) (actual time=382.550..383.696 rows=65 loops=1)"
"                    ->  ProjectSet  (cost=55.25..61.51 rows=1000 width=32) (actual time=382.544..383.652 rows=65 loops=1)"
"                          ->  Aggregate  (cost=55.25..55.26 rows=1 width=32) (actual time=381.323..381.325 rows=1 loops=1)"
"                                ->  Index Scan using relevant_habitats_pkey on relevant_habitats rh  (cost=0.28..28.75 rows=12 width=130244) (actual time=0.026..0.048 rows=12 loops=1)"
"                                      Index Cond: (assessment_area_id = 94)"
"              ->  Index Scan using idx_hexagons_geometry_gist on hexagons h  (cost=0.29..1846.45 rows=53 width=156) (actual time=0.315..0.624 rows=12 loops=65)"
"                    Index Cond: (geometry amp;amp; (((st_dump((st_union(rh.geometry))))).geom))"
"                    Filter: (((zoom_level)::integer = 1) AND st_intersects((((st_dump((st_union(rh.geometry))))).geom), geometry))"
"                    Rows Removed by Filter: 19"
"Planning Time: 0.390 ms"
"Execution Time: 1372.443 ms"
 

Обновление 2: теперь вывод EXPLAIN ANALYZE на втором select ( SELECT h.receptor_id~ ) и CTE заменен таблицей (temp).:

 "QUERY PLAN"
"GroupAggregate  (cost=2691484.47..20931829.74 rows=390140 width=36) (actual time=29.455..927.945 rows=685 loops=1)"
"  Group Key: h.receptor_id"
"  ->  Sort  (cost=2691484.47..2693288.89 rows=721768 width=188) (actual time=28.382..31.514 rows=806 loops=1)"
"        Sort Key: h.receptor_id"
"        Sort Method: quicksort  Memory: 336kB"
"        ->  Nested Loop  (cost=0.29..2488035.20 rows=721768 width=188) (actual time=0.189..27.852 rows=806 loops=1)"
"              ->  Seq Scan on rh_union rhu  (cost=0.00..23.60 rows=1360 width=32) (actual time=0.016..0.050 rows=65 loops=1)"
"              ->  Index Scan using idx_hexagons_geometry_gist on hexagons h  (cost=0.29..1828.89 rows=53 width=156) (actual time=0.258..0.398 rows=12 loops=65)"
"                    Index Cond: (geometry amp;amp; rhu.geometry)"
"                    Filter: (((zoom_level)::integer = 1) AND st_intersects(rhu.geometry, geometry))"
"                    Rows Removed by Filter: 19"
"Planning Time: 0.481 ms"
"Execution Time: 928.583 ms"
 

Комментарии:

1. Применяются ли пространственные индексы?

2. Оба ( h и rh ) имеют индекс gist по геометрии.

3. @romepa можете ли вы добавить план запроса : explain analyse ?

4. @JimJones, я добавил АНАЛИЗ ОБЪЯСНЕНИЙ. Нужно было придумать, где на самом деле это разместить! Полезно ли это?

5. Не могли бы вы создать таблицу вместо предложения CTE и повторно выполнить функцию?