Poszukuję najszybszego rozwiązania do analizy Point in Polygon na 200 milionów punktów [zamknięte]

Mam plik CSV zawierający 200 milionów obserwacji w następującym formacie:

id,x1,y1,x2,y2,day,color
1,"-105.4652334","39.2586939","-105.4321296","39.2236632","Monday","Black"
2,"-105.3224523","39.1323299","-105.4439944","39.3352235","Tuesday","Green"
3,"-104.4233452","39.0234355","-105.4643990","39.1223435","Wednesday","Blue"

Dla każdego zestawu współrzędnych (x1 / y1 i x2 / y2) chcę przypisać US Census Tract lub Census Tract, w który wchodzi (pobrałem plik kształtu TIGER spisu ludności tutaj: ftp://ftp2.census.gov/ geo / tiger / TIGER2011 / TRACT / tl_2011_08_tract.zip ). Muszę więc dwukrotnie wykonać operację punkt-wielokąt dla każdej obserwacji. Ważne jest, aby mecze były bardzo dokładne.

Jaki jest najszybszy sposób to zrobić, w tym czas na naukę oprogramowania? Mam dostęp do komputera z 48 GB pamięci - na wypadek, gdyby było to istotne ograniczenie.

Kilka wątków zaleca używanie PostGIS lub Spatialite (Spatialite wygląda na łatwiejszy w użyciu - ale czy jest tak wydajny jak PostGIS?). Jeśli są to najlepsze opcje, czy konieczne jest wypełnienie indeksu przestrzennego (RTree?)? Jeśli tak, jak to zrobić (np. Używając pliku kształtu Traktu Spisu)? Byłbym bardzo wdzięczny za wszelkie zalecenia, które zawierają przykładowy kod (lub wskaźnik do przykładowego kodu).

Moja pierwsza próba (przed znalezieniem tej strony) polegała na użyciu ArcGIS do wykonania połączenia przestrzennego (tylko x1 / y1) podpróbki danych (100 000 punktów) w US Census Block. Zajęło mi to 5 godzin, zanim zabiłem ten proces. Mam nadzieję na rozwiązanie, które można wdrożyć dla całego zestawu danych w czasie krótszym niż 40 godzin obliczeniowych.

Przepraszamy za zadane pytanie - przeczytałem odpowiedzi i zastanawiam się, jak wdrożyć zalecenia. Nigdy nie korzystałem z SQL, Python, C i korzystałem z ArcGIS tylko raz - jestem kompletnym początkującym.

ST_DWithin był szybszy w moim teście niż ST_Intersects. Jest to zaskakujące, zwłaszcza że przygotowany algorytm geometrii powinien uruchamiać takie przypadki. Myślę, że jest szansa, że ​​będzie to o wiele szybsze, niż pokazałem tutaj.

Zrobiłem kilka testów i dwie rzeczy prawie 10-krotnie podwoiły prędkość. Najpierw próbowałem na nowszym komputerze, ale wciąż dość zwyczajnym laptopie, może poza dyskami SSD SATA3.

Poniższe zapytanie zajęło 18 sekund zamiast 62 sekund na starym laptopie. Następnie odkryłem, że całkowicie się myliłem, kiedy pisałem, że indeks w tabeli punktów nie jest konieczny. Z tym indeksem ST_Intersects zachowywał się zgodnie z oczekiwaniami i wszystko stało się bardzo szybkie. Zwiększyłem liczbę punktów w tabeli punktów do 1 miliona punktów i zapytanie:

CREATE TABLE points_ct AS
SELECT imported_ct.gid as ct_gid, t.gid as point_id 
FROM imported_ct , t WHERE ST_Intersects(imported_ct.geom , t.geom);

działa w ciągu 72 sekund. Ponieważ istnieje 1249 wielokątów, 1249000000 testów wykonuje się w 72 sekundy. To daje około 17000000 testów na sekundę. Lub testowanie prawie 14000 punktów w stosunku do wszystkich wielokątów na sekundę.

Z tego testu twoje 400000000 punktów do przetestowania powinno zająć około 8 godzin bez żadnego problemu z rozłożeniem obciążenia na kilka rdzeni. PostGIS nigdy nie przestaje mnie imponować :-)

Po pierwsze, aby zwizualizować wynik, możesz dodać geometrię punktów do tabeli wynikowej, na przykład otworzyć ją w QGIS i nadać jej unikalne wartości w polu import_ct.

Po drugie, tak, możesz również uzyskać punkty wypadające poza dowolnym wielokątem, używając połączenia prawego (lub lewego) w następujący sposób:

CREATE TABLE points_ct AS
SELECT imported_ct.gid as ct_gid, t.gid as point_id 
FROM imported_ct right join t ON ST_Intersects(imported_ct.the_geom , t.geom);

Zrobiłem kilka testów, aby sprawdzić, czy wydaje się to możliwe PostGIS.

Po pierwsze coś, czego nie rozumiem. Masz dwa punkty na rząd. Czy oba punkty są zawsze w tym samym wielokącie? Następnie wystarczy wykonać obliczenia na jednym z punktów. Jeśli mogą znajdować się w dwóch różnych wielokątach, potrzebny będzie sposób połączenia jednego rzędu punktów z dwoma wielokątami.

Z testów wydaje się to wykonalne, ale możesz potrzebować kreatywnego rozwiązania, aby rozłożyć obciążenie na więcej niż jeden rdzeń procesora.

Testowałem na 4-letnim laptopie z dwurdzeniowym procesorem centrino (chyba około 2,2 GHz), 2 GB pamięci RAM. Jeśli masz 48 BG RAM, myślę, że masz też o wiele więcej mocy procesora.

Stworzyłem losową tabelę punktów z 100 000 punktami:

CREATE TABLE t AS
WITH r AS
(SELECT ST_Extent(the_geom)::geometry ext FROM imported_ct)
SELECT ST_Point(x,y) AS geom FROM 
(SELECT GENERATE_SERIES(1,100000)) s,
(SELECT ST_Xmin(ext)+(random()*(ST_Xmax(ext)-ST_Xmin(ext))) x, ST_Ymin(ext)+(random()*(ST_Ymax(ext)-ST_Ymin(ext))) y FROM r
) f;

Następnie dodaj gid jak:

ALTER TABLE t ADD COLUMN GID SERIAL;

Następnie uruchom:

CREATE TABLE points_ct AS
SELECT imported_ct.gid as ct_gid, t.gid as point_id FROM imported_ct , t WHERE ST_Dwithin(imported_ct.the_geom , t.geom,0);

zajmuje około 62 sekund (porównaj z wynikiem ArcGIS z taką samą ilością punktów). Rezultatem jest tabela łącząca punkty w mojej tabeli t z gid w tabeli z obszarem spisu.

Przy tej prędkości osiągniesz 200 punktów w ciągu około 34 godzin. Tak więc, jeśli wystarczy sprawdzić jeden punkt, mój stary laptop może to zrobić z jednym rdzeniem.

Ale jeśli musisz sprawdzić oba punkty, może być trudniej.

Następnie można ręcznie rozłożyć obciążenie na więcej niż jeden rdzeń, uruchamiając wiele sesji z bazą danych i uruchamiając różne zapytania.

W moim przykładzie z 50000 punktami i dwoma rdzeniami procesora próbowałem:

CREATE TABLE t1 as
SELECT imported_ct.gid as ct_gid, t.gid as point_id FROM imported_ct , t WHERE t.gid >50000 and  ST_Dwithin(imported_ct.the_geom , t.geom,0);

na jednej sesji db w tym samym czasie, co uruchomienie:

CREATE TABLE t2 as
SELECT imported_ct.gid as ct_gid, t.gid as point_id FROM imported_ct , t WHERE t.gid <=50000 and  ST_Dwithin(imported_ct.the_geom , t.geom,0);

w innej sesji db.

Zajęło to około 36 sekund, więc jest nieco wolniejsze niż pierwszy przykład, prawdopodobnie w zależności od zapisu na płycie w tym samym czasie. Ale ponieważ rdzenie rdzeniowe działają w tym samym czasie, nie zajęło mi to więcej niż 36 sekund.

Aby połączyć tabelę t1 i t2, spróbuj:

CREATE TABLE t3 AS 
SELECT * FROM t1
UNION ALL
SELECT * FROM t2;

używając około pół sekundy.

Tak więc przy świeższym sprzęcie i rozkładzie obciążenia na wiele rdzeni powinno to być absolutnie możliwe, nawet jeśli rzeczywisty świat będzie wolniejszy niż przypadek testowy.

Warto zauważyć, że przykład pochodzi z systemu Linux (Ubuntu). Korzystanie z systemu Windows będzie inną historią. Ale mam wszystkie inne codzienne aplikacje, więc laptop jest dość mocno obciążony. To może symulować obudowę systemu Windows całkiem dobrze, nie otwierając niczego poza pgadmin.

Prawdopodobnie najłatwiej jest z PostGIS. W Internecie znajduje się kilka samouczków dotyczących importowania danych punktów csv / txt do PostGIS. Link1

Nie jestem pewien wydajności wyszukiwania punkt-w-wielokącie w PostGIS; powinno być szybsze niż ArcGIS. Indeks przestrzenny GIST wykorzystywany przez PostGIS jest dość szybki. Link2 Link3

Możesz również przetestować indeks geoprzestrzenny MongoDB . Ale to wymaga trochę więcej czasu, aby zacząć. Wierzę, że MongoDB może być naprawdę szybki. Nie testowałem tego przy wyszukiwaniu wielokątów, więc nie jestem pewien.