Optymalizacja bardzo dużego punktu w zapytaniu wielokąta

Mam krajowy zestaw danych punktów adresowych (37 milionów) i wielokątny zbiór konturów powodzi (2 miliony) typu MultiPolygonZ, niektóre z wielokątów są bardzo złożone, maks. ST_NPoints wynosi około 200 000. Próbuję zidentyfikować za pomocą PostGIS (2.18), które punkty adresowe znajdują się w wielokącie powodziowym, i zapisać je w nowej tabeli z identyfikatorem adresu i szczegółami ryzyka powodzi. Próbowałem z perspektywy adresu (ST_Within), ale potem zamieniłem to, zaczynając od perspektywy obszaru powodziowego (ST_Contains), uzasadniając to tym, że są duże obszary w ogóle bez ryzyka powodzi. Oba zestawy danych zostały ponownie przerzucone na 4326 i obie tabele mają indeks przestrzenny. Moje zapytanie poniżej działa od 3 dni i nie pokazuje żadnych oznak ukończenia w najbliższym czasie!

select a.id, f.risk_factor_1, f.risk_factor_2, f.risk_factor_3
into gb.addresses_with_flood_risk
from gb.flood_risk_areas f, gb.addresses a
where ST_Contains(f.the_geom, a.the_geom);

Czy istnieje bardziej optymalny sposób na uruchomienie tego? Ponadto, w przypadku długich zapytań tego typu, jaki jest najlepszy sposób monitorowania postępu poza analizowaniem wykorzystania zasobów i pg_stat_activity?

Moje oryginalne zapytanie zakończyło się OK, chociaż przez 3 dni, i zostałem odsunięty na bok z inną pracą, więc nigdy nie musiałem poświęcać czasu na wypróbowanie rozwiązania. Jednak właśnie ponownie to odwiedziłem i zapoznałem się z zaleceniami, jak dotąd jak dotąd. Użyłem następujących:

1. Utworzono 50-kilometrową siatkę nad Wielką Brytanią za pomocą zaproponowanego tutaj rozwiązania ST_FishNet
2. Ustaw SRID wygenerowanej siatki na British National Grid i zbuduj na niej indeks przestrzenny
3. Obciąłem moje dane powodziowe (MultiPolygon) za pomocą ST_Intersection i ST_Intersects (tylko tutaj musiałem użyć ST_Force_2D na geomie, ponieważ shape2pgsql dodał indeks Z
4. Obcięto moje dane punktów przy użyciu tej samej siatki
5. Utworzono indeksy dla wiersza i kolumny oraz indeks przestrzenny dla każdej z tabel

Jestem teraz gotowy do uruchomienia skryptu, będę iterował wiersze i kolumny wypełniające wyniki w nowej tabeli, dopóki nie obejmę całego kraju. Ale właśnie sprawdziłem moje dane powodziowe i niektóre z największych wielokątów wydają się być zagubione w tłumaczeniu! To jest moje zapytanie:

SELECT g.row, g.col, f.gid, f.objectid, f.prob_4band, ST_Intersection(ST_Force_2D(f.geom), g.geom) AS geom 
INTO rofrse.tmp_flood_risk_grid 
FROM rofrse.raw_flood_risk f, rofrse.gb_grid g
WHERE (ST_Intersects(ST_Force_2D(f.geom), g.geom));

Moje oryginalne dane wyglądają tak:

Jednak po wycinaniu wygląda to tak:

To jest przykład „brakującego” wielokąta:

Aby najpierw odpowiedzieć na ostatnie pytanie, zobacz ten posto potrzebie monitorowania postępu zapytań. Problem jest trudny i zostałby złożony w zapytaniu przestrzennym, ponieważ wiedza, że ​​99% adresów zostało już przeskanowanych w celu przechowania w wielokącie powodziowym, który można uzyskać z licznika pętli w implementacji skanowania tabeli, niekoniecznie pomoc, jeśli 1% adresów przecina wielokąt powodziowy z największą liczbą punktów, podczas gdy poprzednie 99% przecinają niewielki obszar. Jest to jeden z powodów, dla których EXPLAIN może czasami być nieprzydatny w przestrzeni, ponieważ daje wskazanie wierszy, które będą skanowane, ale z oczywistych powodów nie bierze pod uwagę złożoności wielokątów (a zatem dużej części czasu wykonywania) wszelkich zapytań typu przecięcie / skrzyżowanie.

Drugi problem polega na tym, że jeśli spojrzysz na coś takiego

EXPLAIN 
SELECT COUNT(a.id) 
FROM sometable a, someothertable b
WHERE ST_Intersects (a.geom, b.geom)

zobaczysz coś takiego, po pominięciu wielu szczegółów:

_st_intersects(a.geom, b.geom)
   ->  Bitmap Index Scan on ix_spatial_index_name  (cost...rows...width...))
   Index Cond: (a.geom && geom)

Ostateczny warunek, &&, oznacza sprawdzenie pola granicznego przed wykonaniem dokładniejszego przecięcia rzeczywistych geometrii. Jest to oczywiście rozsądne i leży u podstaw działania R-Trees. Jednak w przeszłości pracowałem również nad danymi powodziowymi w Wielkiej Brytanii, więc znam strukturę danych, jeśli (wielokątne) wielokąty są bardzo rozległe - problem ten jest szczególnie dotkliwy, jeśli rzeka płynie, powiedzmy, 45 stopnie - otrzymujesz ogromne ramki ograniczające, które mogą wymusić sprawdzenie dużej liczby potencjalnych skrzyżowań na bardzo złożonych wielokątach.

Jedynym rozwiązaniem, jakie udało mi się wymyślić dla „mojego zapytania działało przez 3 dni i nie wiem, czy mamy 1% czy 99%” problemu, jest zastosowanie pewnego rodzaju podziału i podboju dla manekinów Podejście, które mam na myśli, podziel obszar na mniejsze części i uruchom je osobno, albo w pętli w pliku plpgsql, albo jawnie w konsoli. Ma to tę zaletę, że przecina złożone wielokąty na części, co oznacza, że ​​kolejny punkt w kontrolach wielokątów działa na mniejszych wielokątach, a obwiednia wielokątów jest znacznie mniejsza.

Udało mi się uruchomić zapytania w ciągu jednego dnia, dzieląc Wielką Brytanię na 50 km na 50 km bloków, po zabiciu zapytania, które działało przez ponad tydzień w całej Wielkiej Brytanii. Nawiasem mówiąc, mam nadzieję, że powyższe zapytanie to CREATE TABLE lub UPDATE, a nie tylko SELECT. Kiedy aktualizujesz jedną tabelę, adresy, oparte na byciu w wieloboku powodziowym, będziesz musiał przeskanować całą aktualizowaną tabelę, adresy tak czy inaczej, tak więc posiadanie indeksu przestrzennego na niej wcale nie pomaga.

EDYCJA: Na podstawie tego, że obraz jest wart tysiąca słów, oto obraz niektórych danych o powodzi w Wielkiej Brytanii. Jest jeden bardzo duży wielobok, którego obwiednia pokrywa cały ten obszar, więc łatwo jest zobaczyć, jak na przykład, przecinając najpierw wielokąt powodziowy z czerwoną siatką, kwadrat w południowo-zachodnim rogu zostałby nagle tylko przetestowany na niewielki podzbiór wielokąta.