r punktów w wielokątach

Question

Mam milion punktów i plik dużego kształtu - 8 GB - który jest zbyt duży, aby można go było wczytać do pamięci w R w moim systemie. Plik kształtu jest jednowarstwowy, więc dany x, y uderzy najwyżej w jeden wielokąt - o ile nie jest dokładnie na granicy! Każdy wielokąt jest oznaczony jako severity - np. 1, 2, 3. Używam R na 64-bitowym urządzeniu ubuntu z pamięcią RAM 12 GB.

Jaki jest najprostszy sposób, aby móc „tag” ramka danych do wielokąta severity abym uzyskać data.frame z dodatkową kolumnę, tj x, y, severity?

Answer 1

Tylko dlatego, że wszystko, co masz, to młotek, nie oznacza, że ​​każdy problem jest gwóźdź.

Załaduj dane do PostGIS, utwórz indeks przestrzenny dla wielokątów i wykonaj pojedynczą nakładkę przestrzenną SQL. Eksportuj wyniki do R.

Przy okazji, mówiąc, że shapefile ma rozmiar 8 Gb, nie jest to bardzo przydatna informacja. Pliki kształtów są tworzone z co najmniej trzech plików .shp, który jest geometrią, .dbf, który jest bazą danych, oraz .shx, który łączy te dwa. Jeśli plik .dbf ma wielkość 8 Gb, możesz łatwo odczytać kształty, zastępując je innym plikiem .dbf. Nawet jeśli .shp ma wielkość 8 Gb, może to być tylko trzy wielokąty, w którym to przypadku może być łatwo je uprościć. Ile masz wielokątów i jak duża jest .shp część kształtu?

Answer 2

Myślę, że należy wstępnie przetworzyć dane i utworzyć strukturę, która zawiera listę możliwych wielokątów dla prostokątnych obszarów w siatce. W ten sposób możesz zredukować wielokąty, które musisz sprawdzić w stosunku do punktów, a dodatkowa struktura będzie pasować do pamięci, ponieważ ma tylko wskaźniki do wielokątów.

Oto obraz, aby zilustrować.

Chcesz sprawdzić, które wielokąta żółty punkt w Byłbyś zazwyczaj sprawdzić na wszystkich wielokątów, ale z optymalizacją sieci (linie pomarańczowy, nie narysował całej siatki, tylko jednego z jej pól) wystarczy sprawdzić wypełnione wielokąty, ponieważ wszystkie one są wewnątrz lub częściowo wewnątrz pola siatki.

Podobnym sposobem nie byłoby przechowywanie wszystkich danych wielokąta w pamięci, ale tylko wielokąty ograniczające pola, które wymagałyby tylko 2 zamiast 3 par X/Y dla każdego wielokąta (i dodatkowy wskaźnik do rzeczywistych danych wielokąta), ale nie oszczędza to tyle miejsca, co pierwsza sugestia.

Answer 3

Nie mam naprawdę dobrej odpowiedzi, ale pozwól mi rzucić pomysł. Czy potrafisz odwrócić problem i zamiast pytać, do którego punktu pasuje każdy punkt, zamiast tego "jakie punkty znajdują się w każdym poli?" Może zdołasz zniszczyć swój shapefile np. Do 2000 hrabstw, a następnie stopniowo zdobywać każde hrabstwo i sprawdzać każdy punkt, aby sprawdzić, czy jest on w tym hrabstwie. Jeśli punkt znajduje się w danym hrabstwie, oznaczasz go tagiem i następnym razem wyłączysz go z wyszukiwania.

Wzdłuż tych samych linii można podzielić plik kształtu na 4 regiony. Następnie możesz umieścić pojedynczy region i wszystkie swoje punkty w pamięci. Następnie po prostu powtórz czterokrotne przetwarzanie danych.

Innym pomysłem byłoby użycie narzędzia GIS do obniżenia rozdzielczości (liczby węzłów i wektorów) kształtu pliku shape. To oczywiście zależy od tego, jak ważna jest dokładność w twoim przypadku użycia.

Answer 4

Byłem zainteresowany, aby to zobaczyć i zastanawiałem się, czy zrobiłeś jakiekolwiek postępy na tym froncie. Ponieważ zadałeś pytanie, wyobrażam sobie twój sprzęt komputerowy i oprogramowanie, z którego możesz zrobić to stosunkowo prosta operacja poprawiła się nieco do punktu, w którym rozwiązanie (jeśli nadal potrzebne!) Może być dość proste, chociaż może to zająć dużo czasu przetwarzać milion punktów. Możesz wypróbować coś takiego:

# Load relevant libraries 
library(sp) 
library(maptools) 
library(spatstat) 

# Read shapefile. Hopefully you have a .prj file with your .shp file 
# otherwise you need to set the proj4string argument. Don't inlcude 
# the .shp extension in the filename. I also assume that this will 
# create a SpatialPolygonsDataFrame with the "Severity" attribute 
# attached (from your .dbf file). 
myshapefile <- readShapePoly("myshapefile_without_extension",  proj4string=CRS("+proj=latlong +datum=WGS84")) 


# Read your location data in. Here I assume your data has two columns X and Y giving  locations 
# Remeber that your points need to be in the same projection as your shapefile. If they aren't 
# you should look into using spTransform() on your shapefile first. 
mylocs.df <- read.table(mypoints.csv, sep=",", h=TRUE) 

# Coerce X and Y coordinates to a spatial object and set projection to be the same as 
# your shapefile (WARNING: only do this if you know your points and shapefile are in 
# the same format). 
mylocs.sp <- SpatialPoints(cbind(mylocs.df$X,mylocs.df$Y),  proj4string=CRS(proj4string(myshapefile)) 

# Use over() to return a dataframe equal to nrows of your mylocs.df 
# with each row corresponding to a point with the attributes from the 
# poylgon in which it fell. 
severity.df <- over(mylocs.sp, myshapefile) 

Mam nadzieję, że ten framework da ci to, czego chcesz. To, czy możesz to zrobić za pomocą dostępnego teraz komputera/pamięci RAM, to już inna sprawa!

Answer 5

Dałbym fastshp pakiet spróbować. W moich pobieżnych testach znacząco bije other methods dla shapefiles reading. I ma wyspecjalizowaną funkcję inside, którą moja dobrze pasuje do twoich potrzeb.

Kodeks powinien być w jakiś sposób podobny do:

shp <- read.shp(YOUR_SHP_FILE, format="polygon")) 
inside(shp, x, y) 

gdzie x i y są współrzędnymi.

Jeśli to nie zadziała, wybrałbym rozwiązanie PostGIS wspomniane przez @Spacedman.

Answer 6

Aby odpowiedzieć na moje własne pytanie ... i dziękuję wszystkim za pomoc - Ostatecznym rozwiązaniem było użycie gdala z Pythona, który był stosunkowo łatwo dostosowany zarówno do rastrów, jak i plików kształtów. Niektóre rastry pobiegły do ​​około 40 gb, a niektóre pliki kształtu przekroczyły 8 gb - więc nie było możliwości, by zmieściły się w pamięci na którejkolwiek z maszyn, które mieliśmy w tym czasie (Teraz mam dostęp do maszyny z 128-bitowym ramkiem - ale przeniosłem się na nowe pastwiska!). Kod Pythona/Gdala był w stanie oznaczyć 27 milionów punktów od 1 minuty do 40 minut w zależności od rozmiarów wielokątów w plikach shape - jeśli było dużo małych wielokątów, było to oszałamiająco szybko - gdyby istniało masywne (250k punktów) wielokąty w shapefiles był oszałamiająco wolny! Jednak, aby to porównać, używaliśmy go poprzednio w przestrzennej bazie danych oracle i zajęłoby to około 24 godzin + by oznaczyć 27 milionów punktów, albo rasteryzacja i tagowanie zajęłoby około godziny. Jak zasugerował Spacedman, spróbowałem użyć postgis na moim komputerze z ssd, ale czas na odwrócenie był trochę wolniejszy niż użycie pythona/gdala, ponieważ ostateczne rozwiązanie nie wymagało załadowania plików shape do postgis. Więc podsumować, najszybszym sposobem, aby to zrobić używał Python/gdal:

• kopiować pliki kształt i X, Y CSV SSD
• zmodyfikować plik konfiguracyjny dla skryptu Pythona powiedzieć gdzie znajdują się pliki byli i którego warstwa oznaczyć przed
• perspektywie kilku warstw równolegle - jak to było cpu ogranicza zamiast i/o ograniczone