г точек в многоугольниках

Question

У меня есть миллион точек и большой файл формы - 8 ГБ, который слишком велик для загрузки в память в R в моей системе.Файл формы является однослойным, поэтому заданные x, y попадут не более чем в один многоугольник - если он не находится точно на границе!Каждый многоугольник помечен severity - например, 1, 2, 3.Я использую R на 64-битной машине с Ubuntu с оперативной памятью 12 Гб.с дополнительным столбцом, например x, y, severity?

Answer 1

То, что у вас есть только молоток, не означает, что каждая проблема - это гвоздь.

Загрузите ваши данные в PostGIS, создайте пространственный индекс для ваших полигонов и выполните одно пространственное наложение SQL.Экспорт результатов обратно в R.

Кстати, говорить, что шейп-файл имеет размер 8 Гб, не очень полезная часть информации.Шейп-файлы состоят как минимум из трех файлов: .shp, который является геометрией, .dbf, который является базой данных, и .shx, который соединяет их.Если ваш .dbf имеет размер 8 ГБ, вы можете легко прочитать сами фигуры, заменив его другим .dbf.Даже если .shp равен 8 ГБ, это может быть только три полигона, и в этом случае их будет легко упростить.Сколько полигонов у вас есть, и насколько велика часть .shp шейп-файла?

Answer 2

Я думаю, вам следует предварительно обработать данные и создать структуру, в которой перечислены возможные полигоны для прямоугольных областей в сетке.Таким образом, вы можете уменьшить количество полигонов, которые вы должны будете проверять по точкам, и дополнительная структура поместится в память, так как у нее есть только указатели на полигоны.

Вот изображение для иллюстрации:1004 *

Вы хотите проверить, в каком полигоне находится желтая точка. Обычно вы проверяете по всем полигонам, но с оптимизацией сетки (оранжевые линии не рисуют всю сетку,только одно из его полей) вам нужно только проверить заполненные многоугольники, поскольку все они находятся внутри или частично внутри поля сетки.

Аналогичным способом было бы сохранить не все данные многоугольника в памяти, а толькоограничивающие прямоугольники многоугольников, для которых требуется только 2 вместо 3 пар X / Y для каждого многоугольника (и дополнительный указатель на фактические данные многоугольника), но это не экономит столько места, сколько первое предложение.

Answer 3

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

#   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)

Надеюсь, что рамки могут дать вам то, что вы хотите. Можете ли вы сделать это с имеющимся у вас компьютером / оперативной памятью - это другой вопрос!

Answer 4

У меня нет действительно хорошего ответа, но позвольте мне высказать идею. Можете ли вы перевернуть проблему и вместо того, чтобы спрашивать, к какому полигону подходит каждая точка, вместо этого: «какие точки есть в каждом поли?» Возможно, вы могли бы добавить свой шейп-файл, например, в 2000 округов, а затем постепенно брать каждый округ и проверять каждую точку, чтобы увидеть, находится ли он в этом округе. Если точка находится в данном округе, вы помечаете ее и в следующий раз исключаете из поиска.

В том же духе вы можете разбить шейп-файл на 4 области. Затем вы можете поместить один регион плюс все ваши точки в память. Затем просто повторите 4 раза данные.

Другая идея заключается в использовании ГИС-инструмента для уменьшения разрешения (числа узлов и векторов) шейп-файла. Конечно, это зависит от того, насколько важна точность для вашего варианта использования.

Answer 5

Чтобы ответить на мой собственный вопрос ... и в благодарность всем за помощь - окончательное решение было использовать gdal из python, который был относительно легко адаптирован как для растров, так и для файлов форм.Некоторые из растров занимали около 40 ГБ, а некоторые из файлов форм превышали 8 ГБ - поэтому не было никакого способа, которым они поместились бы в памяти на любой из машин, которые у нас были в то время (Теперь у меня есть доступ к машине с ОЗУ 128 ГБ- но я перешел на новые пастбища!).Код python / gdal смог пометить 27 миллионов точек от 1 минуты до 40 минут в зависимости от размеров многоугольников в шейп-файлах - если было много маленьких полигонов, это было потрясающе быстро - если было несколько больших (250 тысяч точек)) Полигоны в шейп-файлах были потрясающе медленными!Однако, для сравнения, раньше мы запускали его в пространственной базе данных оракула, и для того, чтобы пометить 27 миллионов точек, потребуется около 24 часов, иначе растеризация и пометка займет около часа.Как предложил Spacedman, я попытался использовать postgis на своей машине с ssd, но время оборота было немного медленнее, чем при использовании python / gdal, так как окончательное решение не требовало загрузки шейп-файлов в postgis.Подводя итог, можно сказать, что самым быстрым способом сделать это было использование Python / gdal:

• скопируйте файлы форм и x, y csv в ssd
• измените файл конфигурации для pythonСкрипт, чтобы сказать, где находятся файлы и какой слой пометить для
• , запустить пару слоев параллельно - так как это было ограничено процессором, а не вводом / выводом

Answer 6

Я бы попробовал fastshp . В моих поверхностных тестах он значительно превосходит другие методы для чтения шейп-файлов. И он имеет специальную внутри функцию, которая вполне может удовлетворить ваши потребности.

Код должен быть как-то похож на:

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

где x и y - координаты.

Если это не сработает, я бы выбрал решение PostGIS, упомянутое @ Spacedman.