Ищем самое быстрое решение для анализа Point in Polygon на 200 миллионов точек [закрыто]

У меня есть CSV, содержащий 200 миллионов наблюдений в следующем формате:

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"

Для каждого набора координат (x1 / y1 и x2 / y2) я хочу назначить американский переписной тракт или блок переписи, в который он попадает (я скачал шейп-файл TIGER тракта переписи здесь: ftp://ftp2.census.gov/ geo / tiger / TIGER2011 / TRACT / tl_2011_08_tract.zip ). Итак, мне нужно выполнить операцию точка-полигон дважды для каждого наблюдения. Важно, чтобы совпадения были очень точными.

Какой самый быстрый способ сделать это, включая время на изучение программного обеспечения? У меня есть доступ к компьютеру с 48 ГБ памяти - в случае, если это может быть уместным ограничением.

Несколько потоков рекомендуют использовать PostGIS или Spatialite (Spatialite выглядит проще в использовании - но так ли он эффективен, как PostGIS?). Если это лучшие варианты, обязательно ли заполнять пространственный индекс (RTree?)? Если да, то как это сделать (например, с помощью шейп-файла Census Tract)? Я был бы чрезвычайно признателен за любые рекомендации, которые включают пример кода (или указатель на пример кода).

Моя первая попытка (до нахождения этого сайта) заключалась в использовании ArcGIS для пространственного соединения (только x1 / y1) подвыборки данных (100 000 точек) в блоке переписи США. Это заняло более 5 часов, прежде чем я убил процесс. Я надеюсь на решение, которое может быть реализовано на всем наборе данных менее чем за 40 часов вычислительного времени.

Извиняюсь за задание вопроса, который был задан ранее - я прочитал ответы, и у меня остается вопрос, как выполнить рекомендации. Я никогда не использовал SQL, Python, C и только однажды использовал ArcGIS - я начинающий.

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

Я сделал еще несколько тестов и две вещи почти в 10 раз увеличили скорость. Сначала я попробовал на более новом компьютере, но все еще довольно обычный ноутбук, может быть, кроме SATA3-дисков ssd.

Тогда запрос ниже занял 18 секунд вместо 62 секунд на старом ноутбуке. Затем я обнаружил, что был совершенно неправ, когда писал, что указатель на таблицу точек не нужен. С этим индексом ST_Intersects вел себя как ожидалось, и все стало очень быстро. Я увеличил количество точек в таблице баллов до 1 миллиона точек и запрос:

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

работает за 72 секунды. Поскольку существует 1249 полигонов, 1249000000 тестов выполняются за 72 секунды. Это составляет около 17000000 тестов в секунду. Или тестирование почти 14000 точек против всех полигонов в секунду.

Из этого теста ваши 400000000 точек для тестирования должны занять около 8 часов без каких-либо проблем с распределением нагрузки на несколько ядер. PostGIS не перестает меня удивлять :-)

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

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

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

Я сделал несколько тестов, чтобы проверить, кажется ли это возможным PostGIS.

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

Из тестов это кажется выполнимым, но вам может понадобиться какое-то креативное решение, чтобы распределить нагрузку более чем на одно ядро ​​процессора.

Я тестировал 4-летний ноутбук с двухъядерным процессором Centrino (примерно 2,2 ГГц), 2 ГБ оперативной памяти. Если у вас 48 ГБ ОЗУ, я думаю, у вас гораздо больше процессора.

Что я сделал, чтобы создать таблицу случайных точек с 100000 точек, как это:

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;

Затем добавив gid, как:

ALTER TABLE t ADD COLUMN GID SERIAL;

Затем работает:

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

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

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

Но если вам нужно проверить обе точки, это может быть сложнее.

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

В моем примере с 50000 точками и двумя процессорами я попробовал:

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

на одной db-сессии одновременно с запуском:

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

на другой сессии дб.

Это заняло около 36 секунд, так что это немного медленнее, чем в первом примере, возможно, в зависимости от одновременной записи на диск. Но поскольку ядра работают одновременно, это заняло не более 36 секунд моего времени.

Чтобы объединить таблицу t1 и t2 в качестве:

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

используя около половины секунды.

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

Стоит отметить, что пример взят из Linux (Ubuntu). Использование Windows будет другой историей. Но у меня все остальные ежедневные приложения работают, поэтому ноутбук довольно сильно загружен. Так что это может очень хорошо имитировать случай Windows, не открывая ничего, кроме pgadmin.

Вероятно, самый простой способ с PostGIS. В Интернете есть несколько руководств по импорту данных точек csv / txt в PostGIS. Link1

Я не уверен насчет выполнения поиска по точкам в многоугольнике в PostGIS; это должно быть быстрее, чем ArcGIS. Пространственный индекс GIST, который использует PostGIS, довольно быстрый. Link2 Link3

Вы также можете протестировать геопространственный индекс MongoDB . Но это требует немного больше времени, чтобы начать. Я считаю, что MongoDB может быть очень быстрым. Я не проверял это с поиском точки-полигона, поэтому не могу быть уверен.