Najszybszy sposób na przestrzenne połączenie punktu CSV za pomocą wielokąta Shapefile

Mam plik CSV o wartości 1 miliarda punktów i plik kształtu z około 5000 wielokątów. Jaki byłby najszybszy sposób na przestrzenne połączenie punktów i wielokątów? Dla każdego punktu muszę uzyskać zawierający identyfikator wielokąta. (Wielokąty się nie pokrywają).

Zazwyczaj ładuję oba zestawy danych do PostGIS. Czy istnieje szybszy sposób na wykonanie pracy?

Szukam rozwiązania typu open source.

Jeśli „najszybszy” obejmuje kwotę swoim czasie, który spędził, rozwiązanie będzie zależeć od tego, jakie oprogramowanie są wygodne i można używać sprawnie. Poniższe uwagi konsekwentnie koncentrują się na pomysłach na osiągnięcie najszybszych możliwych czasów obliczeniowych .

Jeśli używasz programu w puszkach, prawie na pewno najlepsze, co możesz zrobić, to wstępnie przetworzyć wielokąty, aby skonfigurować strukturę danych punkt-w-wielokącie, taką jak drzewo KD lub quadtree, którego wydajność zwykle wynosi O (log (V ) * (N + V)) gdzie V jest całkowitą liczbą wierzchołków w wielokątach, a N jest liczbą punktów, ponieważ struktura danych zajmie co najmniej wysiłek O (log (V) * V), a następnie utworzy należy sondować dla każdego punktu według kosztu jednostkowego O (log (V)).

Możesz zrobić znacznie lepiej, najpierw łącząc siatki wielokątów, wykorzystując założenie braku nakładania się. Każda komórka siatki znajduje się albo całkowicie we wnętrzu wielokąta (w tym we wnętrzu „uniwersalnego wielokąta”), w którym to przypadku oznacz komórkę identyfikatorem wielokąta, albo zawiera jedną lub więcej krawędzi wielokąta. Koszt tej rasteryzacji, równej liczbie komórek siatki, do których się odwołujesz podczas rasteryzacji wszystkich krawędzi, wynosi O (V / c), gdzie c jest rozmiarem komórki, ale stała domyślna w notacji big-O jest niewielka.

(Jedną z zalet tego podejścia jest to, że możesz wykorzystać standardowe procedury graficzne. Na przykład, jeśli masz system, który (a) narysuje wielokąty na ekranie wirtualnym, używając (b) odrębnego koloru dla każdego wielokąta i (c) pozwala możesz odczytać kolor każdego piksela, którym chcesz się zająć, już go wykonałeś.)

Po ustawieniu tej siatki należy wstępnie przeskanować punkty, obliczając komórkę zawierającą każdy punkt (operacja O (1) wymagająca tylko kilku zegarów). O ile punkty nie są skupione wokół granic wielokąta, zwykle pozostawia tylko około punktów O (c) z niejednoznacznymi wynikami. Całkowity koszt budowy sieci i wstępnej kontroli wynosi zatem O (V / c + 1 / c ^ 2) + O (N). Musisz użyć innej metody (takiej jak te zalecane do tej pory), aby przetworzyć pozostałe punkty (to znaczy bliskie granicom wielokątów), kosztem O (log (V) * N * c) .

Gdy c maleje, coraz mniej punktów będzie w tej samej komórce siatki z krawędzią, a zatem coraz mniej będzie wymagało późniejszego przetwarzania O (log (V)). Przeciwdziałając temu, istnieje potrzeba przechowywania komórek siatki O (1 / c ^ 2) i spędzania czasu O (V / c + 1 / c ^ 2) na rasteryzacji wielokątów. Dlatego będzie optymalny rozmiar siatki c. Przy jego użyciu całkowity koszt obliczeniowy wynosi O (log (V) * N), ale stała domyślna jest zwykle znacznie mniejsza niż przy zastosowaniu procedur standardowych, ze względu na szybkość O (N) wstępnego badania przesiewowego.

20 lat temu przetestowałem to podejście (wykorzystując równomiernie rozmieszczone punkty w całej Anglii i na morzu i wykorzystując stosunkowo surową siatkę około 400 000 komórek oferowaną przez bufory wideo tamtych czasów) i uzyskałem przyspieszenie o dwa rzędy wielkości w porównaniu z najlepszym opublikowanym algorytmem, jaki mogłem odnaleźć. Nawet gdy wielokąty są małe i proste (jak trójkąty), masz praktycznie pewność przyspieszenia rzędu wielkości.

Z mojego doświadczenia wynika, że ​​obliczenia były tak szybkie, że cała operacja była ograniczona szybkościami we / wy danych, a nie procesorem. Przewidując, że We / Wy może być wąskim gardłem, można osiągnąć najszybsze wyniki, przechowując punkty w możliwie skompresowanym formacie, aby zminimalizować czas odczytu danych. Zastanów się także, jak należy przechowywać wyniki, aby ograniczyć zapisy na dysku.

Z mojej strony prawdopodobnie załadowałbym dane CSV do pliku shp , a następnie napisałem skrypt Pythona używając shapefile i foremnie, aby uzyskać zawierający identyfikator wielokąta i zaktualizować wartość pola.

Nie wiem, czy narzędzia geotools i JTS są szybsze niż plik kształtu / foremny ... Nie mam czasu na testowanie!

edycja : Nawiasem mówiąc, konwersja csv do formatu pliku kształtu prawdopodobnie nie jest wymagana, ponieważ wartości można łatwo sformatować w celu przetestowania z obiektami przestrzennymi z pliku kształtu wielokąta.

Skończyło się na konwersji wielokątów na raster i próbkowaniu ich w punktach. Ponieważ moje wielokąty nie zachodziły na siebie, a wysoka dokładność nie była konieczna (wielokąty reprezentowały klasy zagospodarowania przestrzennego, a ich granice i tak były uważane za raczej niepewne), było to najbardziej efektywne czasowo rozwiązanie, jakie mogłem wymyślić.

Chciałbym szybko napisać mały program java opartą na czytniku Shapefile z geotools i operacja zawiera od WST . Nie wiem, jak szybko to może być ...

Użyj Spatialite .

Pobierz GUI. Możesz otworzyć zarówno plik Shapefile, jak i CSV jako tabele wirtualne. Oznacza to, że tak naprawdę nie importujesz ich do bazy danych, ale pojawiają się one jako tabele i możesz szybko do nich dołączyć i zapytać je w dowolny sposób.

Możesz to zrobić dość szybko, używając OGR w C / C ++ / Python (Python powinien być najwolniejszy z 3). Zapętlaj wszystkie wielokąty i ustaw filtr na punktach, zapętlaj przefiltrowane punkty, a będziesz wiedział, że każdy punkt, przez który przechodzisz, będzie należeć do bieżącego wielokąta. Oto przykładowy kod w Pythonie przy użyciu OGR, który będzie przechodził przez wielokąty i odpowiednio filtrował punkty. Kod C / C ++ będzie wyglądał dość podobnie do tego, i wyobrażam sobie, że uzyskasz znaczny wzrost prędkości w porównaniu do Pythona. Musisz dodać kilka wierszy kodu, aby zaktualizować plik CSV w miarę postępów:

from osgeo import ogr
from osgeo.gdalconst import *

inPolyDS = ogr.Open("winnipeg.shp", GA_ReadOnly)
inPolyLayer = inPolyDS.GetLayer(0)
inPointDS = ogr.Open("busstops.vrt", GA_ReadOnly)   
inPointLayer = inPointDS.GetLayerByName("busstops")

inPolyFeat = inPolyLayer.GetNextFeature()
while inPolyFeat is not None:
  inPtFeat = inPointLayer.GetNextFeature()
  while inPtFeat is not None:
    ptGeom = inPtFeat.GetGeometryRef()
    # Do work here...

    inPtFeat = inPointLayer.GetNextFeature()

  inPolyFeat = inPolyLayer.GetNextFeature()

Plik VRT (busstops.vrt):

<OGRVRTDataSource>
  <OGRVRTLayer name="busstops">
    <SrcDataSource>busstops.csv</SrcDataSource>
    <GeometryType>wkbPoint</GeometryType>
    <LayerSRS>WGS84</LayerSRS>
    <GeometryField encoding="PointFromColumns" x="X" y="Y" reportSrcColumn="FALSE" />
  </OGRVRTLayer>
</OGRVRTDataSource>

Plik CSV (busstops.csv):

FID,X,Y,stop_name
1,-97.1394781371062,49.8712241633646,Southbound Osborne at Mulvey

Plik CSVT (busstops.csvt, OGR potrzebuje go do identyfikacji typów kolumn, w przeciwnym razie nie wykona filtra przestrzennego):

Integer,Real,Real,String

może spróbować csv2shp csv2shp

ciekawy, w jakiej branży znajduje się CSV o wartości miliarda punktów?