Пространственная кластеризация с PostGIS?

97

Я ищу алгоритм пространственной кластеризации для использования его в базе данных с поддержкой PostGIS для точечных объектов. Я собираюсь написать функцию plpgsql, которая принимает расстояние между точками в пределах одного кластера в качестве входных данных. На выходе функция возвращает массив кластеров. Наиболее очевидным решением является создание буферных зон на указанном расстоянии вокруг объекта и поиск объектов в этом буфере. Если такие функции существуют, продолжайте создавать буфер вокруг них и т. Д. Если такие функции не существуют, это означает, что создание кластера завершено. Может быть есть какие-нибудь умные решения?

drnextgis
источник
4
Существует огромное разнообразие методов кластеризации из-за различной природы данных и разных целей кластеризации. Чтобы получить представление о том, что там есть, и прочитать, что другие делают для кластеризации матриц расстояний, поищите на сайте CV @ SE . На самом деле, «выбор метода кластеризации» является почти точной копией вашего и дает хорошие ответы.
whuber
8
+1 к вопросу, потому что найти реальный пример PostGIS SQL вместо ссылок на алгоритмы невозможно для чего-либо, кроме базовой кластеризации сетки, особенно для более экзотических кластеров, таких как MCL
wildpeaks

Ответы:

112

Существует как минимум два хороших метода кластеризации для PostGIS: k-средства (через kmeans-postgresqlрасширение) или кластеризация геометрии в пределах порогового расстояния (PostGIS 2.2)


1) k -значит сkmeans-postgresql

Установка: Вам нужно иметь PostgreSQL 8.4 или выше на хост-системе POSIX (я не знаю, с чего начать для MS Windows). Если вы установили это из пакетов, убедитесь, что у вас также есть пакеты разработки (например, postgresql-develдля CentOS). Скачать и извлечь:

wget http://api.pgxn.org/dist/kmeans/1.1.0/kmeans-1.1.0.zip
unzip kmeans-1.1.0.zip
cd kmeans-1.1.0/

Перед сборкой необходимо установить USE_PGXS переменную окружения (мой предыдущий пост проинструктировал удалить эту часть Makefile, что было не лучшим вариантом). Одна из этих двух команд должна работать для вашей оболочки Unix:

# bash
export USE_PGXS=1
# csh
setenv USE_PGXS 1

Теперь соберите и установите расширение:

make
make install
psql -f /usr/share/pgsql/contrib/kmeans.sql -U postgres -D postgis

(Примечание: я тоже пробовал это с Ubuntu 10.10, но не повезло, так как путь к нему pg_config --pgxsне существует! Это, вероятно, ошибка упаковки Ubuntu)

Использование / Пример: у вас должна быть где-то таблица точек (я нарисовал кучу псевдослучайных точек в QGIS). Вот пример того, что я сделал:

SELECT kmeans, count(*), ST_Centroid(ST_Collect(geom)) AS geom
FROM (
  SELECT kmeans(ARRAY[ST_X(geom), ST_Y(geom)], 5) OVER (), geom
  FROM rand_point
) AS ksub
GROUP BY kmeans
ORDER BY kmeans;

5я обеспечил во втором аргументе kmeansфункции окна является К целому числу , чтобы произвести пять кластеров. Вы можете изменить это на любое целое число, которое вы хотите.

Ниже приведены 31 псевдослучайная точка, которую я нарисовал, и пять центроидов с меткой, показывающей счет в каждом кластере. Это было создано с использованием вышеуказанного SQL-запроса.

Kmeans


Вы также можете попытаться проиллюстрировать, где находятся эти кластеры, с помощью ST_MinimumBoundingCircle :

SELECT kmeans, ST_MinimumBoundingCircle(ST_Collect(geom)) AS circle
FROM (
  SELECT kmeans(ARRAY[ST_X(geom), ST_Y(geom)], 5) OVER (), geom
  FROM rand_point
) AS ksub
GROUP BY kmeans
ORDER BY kmeans;

Kmeans2


2) Кластеризация в пределах порогового расстояния с ST_ClusterWithin

Эта агрегатная функция включена в PostGIS 2.2 и возвращает массив GeometryCollections, где все компоненты находятся на расстоянии друг от друга.

Вот пример использования, где расстояние 100,0 является порогом, который приводит к 5 различным кластерам:

SELECT row_number() over () AS id,
  ST_NumGeometries(gc),
  gc AS geom_collection,
  ST_Centroid(gc) AS centroid,
  ST_MinimumBoundingCircle(gc) AS circle,
  sqrt(ST_Area(ST_MinimumBoundingCircle(gc)) / pi()) AS radius
FROM (
  SELECT unnest(ST_ClusterWithin(geom, 100)) gc
  FROM rand_point
) f;

ClusterWithin100

Самый большой средний кластер имеет радиус окружающего круга 65,3 единиц или около 130, что больше, чем пороговое значение. Это связано с тем, что отдельные расстояния между геометриями элементов меньше порога, поэтому он связывает их вместе как один более крупный кластер.

Майк Т
источник
2
Отлично, эти модификации помогут при установке :-) Однако, боюсь, я не смогу действительно использовать это расширение в конце, потому что (если я правильно понял), ему нужно жестко закодированное магическое число кластеров, что хорошо для предотвращения статических данных. Вы можете настроить его заранее, но мне не подойдет для кластеризации произвольных (из-за различных фильтров) наборов данных, например, большого разрыва в кластере из 10 точек на последнем изображении. Однако это также поможет другим людям, потому что (afaik), это единственный существующий пример SQL (за исключением одного вкладыша на домашней странице расширения) для этого расширения.
Wildpeaks
(а вы ответили в то же время, что я удалил предыдущий комментарий, чтобы переформулировать его, извините)
wildpeaks
7
Для кластеризации kmeans вам необходимо заранее указать количество кластеров; Мне любопытно, есть ли альтернативные алгоритмы, где количество кластеров не требуется.
Djq
1
Версия 1.1.0 теперь доступна: api.pgxn.org/dist/kmeans/1.1.0/kmeans-1.1.0.zip
djq
1
@ maxd нет. Если A = πr², то r = √ (A / π).
Майк Т
27

Я написал функцию, которая вычисляет кластеры объектов на основе расстояния между ними и строит выпуклую оболочку над этими объектами:

CREATE OR REPLACE FUNCTION get_domains_n(lname varchar, geom varchar, gid varchar, radius numeric)
    RETURNS SETOF record AS
$$
DECLARE
    lid_new    integer;
    dmn_number integer := 1;
    outr       record;
    innr       record;
    r          record;
BEGIN

    DROP TABLE IF EXISTS tmp;
    EXECUTE 'CREATE TEMPORARY TABLE tmp AS SELECT '||gid||', '||geom||' FROM '||lname;
    ALTER TABLE tmp ADD COLUMN dmn integer;
    ALTER TABLE tmp ADD COLUMN chk boolean DEFAULT FALSE;
    EXECUTE 'UPDATE tmp SET dmn = '||dmn_number||', chk = FALSE WHERE '||gid||' = (SELECT MIN('||gid||') FROM tmp)';

    LOOP
        LOOP
            FOR outr IN EXECUTE 'SELECT '||gid||' AS gid, '||geom||' AS geom FROM tmp WHERE dmn = '||dmn_number||' AND NOT chk' LOOP
                FOR innr IN EXECUTE 'SELECT '||gid||' AS gid, '||geom||' AS geom FROM tmp WHERE dmn IS NULL' LOOP
                    IF ST_DWithin(ST_Transform(ST_SetSRID(outr.geom, 4326), 3785), ST_Transform(ST_SetSRID(innr.geom, 4326), 3785), radius) THEN
                    --IF ST_DWithin(outr.geom, innr.geom, radius) THEN
                        EXECUTE 'UPDATE tmp SET dmn = '||dmn_number||', chk = FALSE WHERE '||gid||' = '||innr.gid;
                    END IF;
                END LOOP;
                EXECUTE 'UPDATE tmp SET chk = TRUE WHERE '||gid||' = '||outr.gid;
            END LOOP;
            SELECT INTO r dmn FROM tmp WHERE dmn = dmn_number AND NOT chk LIMIT 1;
            EXIT WHEN NOT FOUND;
       END LOOP;
       SELECT INTO r dmn FROM tmp WHERE dmn IS NULL LIMIT 1;
       IF FOUND THEN
           dmn_number := dmn_number + 1;
           EXECUTE 'UPDATE tmp SET dmn = '||dmn_number||', chk = FALSE WHERE '||gid||' = (SELECT MIN('||gid||') FROM tmp WHERE dmn IS NULL LIMIT 1)';
       ELSE
           EXIT;
       END IF;
    END LOOP;

    RETURN QUERY EXECUTE 'SELECT ST_ConvexHull(ST_Collect('||geom||')) FROM tmp GROUP by dmn';

    RETURN;
END
$$
LANGUAGE plpgsql;

Пример использования этой функции:

SELECT * FROM get_domains_n('poi', 'wkb_geometry', 'ogc_fid', 14000) AS g(gm geometry)

'poi' - имя слоя, 'wkb_geometry' - имя столбца геометрии, 'ogc_fid' - первичный ключ таблицы, 14000 - расстояние кластера.

Результат использования этой функции:

введите описание изображения здесь

drnextgis
источник
Большой! Не могли бы вы добавить пример того, как использовать вашу функцию? Спасибо!
Подземье
1
Я немного изменил исходный код и добавил пример использования функции.
drnextgis
Только что попробовал использовать это на postgres 9.1 и строке "FOR innr IN EXECUTE 'SELECT' || gid || ' КАК гид, '|| geom ||' AS geom FROM tmp WHERE dmn IS NULL 'LOOP "выдает следующую ошибку. Есть идеи ? ОШИБКА: функция с множественным значением, вызываемая в контексте, которая не может принять набор
битбокс
Я не уверен, как использовать этот код в PG (PostGIS n00b) в моей таблице. с чего мне начать понимать этот синтаксис? У меня есть таблица с латами и латами, которую я хочу кластеризовать
mga
Прежде всего вы должны построить geometryстолбец в вашей таблице, чтобы не хранить lonlat отдельно и создать столбец с уникальными значениями (идентификаторами).
drnextgis
10

На сегодняшний день наиболее многообещающим я считаю это расширение для кластеризации K-средних в качестве оконной функции: http://pgxn.org/dist/kmeans/

Однако я еще не смог установить его успешно.


В противном случае для базовой кластеризации сетки вы можете использовать SnapToGrid .

SELECT
    array_agg(id) AS ids,
    COUNT( position ) AS count,
    ST_AsText( ST_Centroid(ST_Collect( position )) ) AS center,
FROM mytable
GROUP BY
    ST_SnapToGrid( ST_SetSRID(position, 4326), 22.25, 11.125)
ORDER BY
    count DESC
;
wildpeaks
источник
2

Дополняю ответ @MikeT ...

Для MS Windows:

Требования:

Что ты сделаешь:

  • Настройте исходный код, чтобы экспортировать функцию kmeans в DLL.
  • Скомпилируйте исходный код с помощью cl.exeкомпилятора, чтобы сгенерировать DLL с kmeansфункцией.
  • Поместите сгенерированную DLL в папку PostgreSQL \ lib.
  • Затем вы можете «создать» (связать) UDF в PostgreSQL с помощью команды SQL.

шаги:

  1. Скачать и установить / извлечь требования.
  2. Откройте kmeans.cв любом редакторе:

    1. После #includeстрок определите макрос DLLEXPORT с помощью:

      #if defined(_WIN32)
          #define DLLEXPORT __declspec(dllexport)
      #else
         #define DLLEXPORT
      #endif
    2. Поставьте DLLEXPORTперед каждой из этих строк:

      PG_FUNCTION_INFO_V1(kmeans_with_init);
      PG_FUNCTION_INFO_V1(kmeans);
      
      extern Datum kmeans_with_init(PG_FUNCTION_ARGS);
      extern Datum kmeans(PG_FUNCTION_ARGS);
  3. Откройте командную строку Visual C ++.

  4. В командной строке:

    1. Перейти к извлеченным kmeans-postgresql.
    2. Установите ваш POSTGRESPATH, мой, например: SET POSTGRESPATH=C:\Program Files\PostgreSQL\9.5
    3. Бегать

      cl.exe /I"%POSTGRESPATH%\include" /I"%POSTGRESPATH%\include\server" /I"%POSTGRESPATH%\include\server\port\win32" /I"%POSTGRESPATH%\include\server\port\win32_msvc" /I"C:\Program Files (x86)\Microsoft SDKs\Windows\v7.1A\Include" /LD kmeans.c "%POSTGRESPATH%\lib\postgres.lib"
  5. Скопируйте kmeans.dllв%POSTGRESPATH%\lib

  6. Теперь запустите команду SQL в вашей базе данных, чтобы «СОЗДАТЬ» функцию.

    CREATE FUNCTION kmeans(float[], int) RETURNS int
    AS '$libdir/kmeans'
    LANGUAGE c VOLATILE STRICT WINDOW;
    
    CREATE FUNCTION kmeans(float[], int, float[]) RETURNS int
    AS '$libdir/kmeans', 'kmeans_with_init'
    LANGUAGE C IMMUTABLE STRICT WINDOW;
caiohamamura
источник
2

Вот способ отображения в QGIS результата запроса PostGIS, данного в 2) в этом ответе

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

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

SELECT id,countfeature,circle FROM (SELECT row_number() over () AS id,
  ST_NumGeometries(gc) as countfeature,
  ST_MinimumBoundingCircle(gc) AS circle
FROM (
  SELECT unnest(ST_ClusterWithin(the_geom, 100)) gc
  FROM rand_point
) f) a WHERE ST_GeometryType(circle) = 'ST_Polygon'

Затем для кластеризованных точек необходимо преобразовать геометрические коллекции в многоточечные:

SELECT row_number() over () AS id,
  ST_NumGeometries(gc) as countfeature,
  ST_CollectionExtract(gc,1) AS multipoint
FROM (
  SELECT unnest(ST_ClusterWithin(the_geom, 100)) gc
  FROM rand_point
) f

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

Кластеризация в QGIS

Николя Буасто
источник
2

Вы можете более легко использовать решение Kmeans с помощью метода ST_ClusterKMeans, который доступен в postgis из 2.3. Пример:

SELECT kmean, count(*), ST_SetSRID(ST_Extent(geom), 4326) as bbox 
FROM
(
    SELECT ST_ClusterKMeans(geom, 20) OVER() AS kmean, ST_Centroid(geom) as geom
    FROM sls_product 
) tsub
GROUP BY kmean;

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

Оригинальные геометрии Функциональные кластеры

Volda
источник
1

Решение для кластеризации снизу вверх Получите один кластер из облака точек с максимальным диаметром в postgis, который не требует динамических запросов.

CREATE TYPE pt AS (
    gid character varying(32),
    the_geom geometry(Point))

и тип с идентификатором кластера

CREATE TYPE clustered_pt AS (
    gid character varying(32),
    the_geom geometry(Point)
    cluster_id int)

Следующая функция алгоритма

CREATE OR REPLACE FUNCTION buc(points pt[], radius integer)
RETURNS SETOF clustered_pt AS
$BODY$

DECLARE
    srid int;
    joined_clusters int[];

BEGIN

--If there's only 1 point, don't bother with the loop.
IF array_length(points,1)<2 THEN
    RETURN QUERY SELECT gid, the_geom, 1 FROM unnest(points);
    RETURN;
END IF;

CREATE TEMPORARY TABLE IF NOT EXISTS points2 (LIKE pt) ON COMMIT DROP;

BEGIN
    ALTER TABLE points2 ADD COLUMN cluster_id serial;
EXCEPTION
    WHEN duplicate_column THEN --do nothing. Exception comes up when using this function multiple times
END;

TRUNCATE points2;
    --inserting points in
INSERT INTO points2(gid, the_geom)
    (SELECT (unnest(points)).* ); 

--Store the srid to reconvert points after, assumes all points have the same SRID
srid := ST_SRID(the_geom) FROM points2 LIMIT 1;

UPDATE points2 --transforming points to a UTM coordinate system so distances will be calculated in meters.
SET the_geom =  ST_TRANSFORM(the_geom,26986);

--Adding spatial index
CREATE INDEX points_index
ON points2
USING gist
(the_geom);

ANALYZE points2;

LOOP
    --If the smallest maximum distance between two clusters is greater than 2x the desired cluster radius, then there are no more clusters to be formed
    IF (SELECT ST_MaxDistance(ST_Collect(a.the_geom),ST_Collect(b.the_geom))  FROM points2 a, points2 b
        WHERE a.cluster_id <> b.cluster_id
        GROUP BY a.cluster_id, b.cluster_id 
        ORDER BY ST_MaxDistance(ST_Collect(a.the_geom),ST_Collect(b.the_geom)) LIMIT 1)
        > 2 * radius
    THEN
        EXIT;
    END IF;

    joined_clusters := ARRAY[a.cluster_id,b.cluster_id]
        FROM points2 a, points2 b
        WHERE a.cluster_id <> b.cluster_id
        GROUP BY a.cluster_id, b.cluster_id
        ORDER BY ST_MaxDistance(ST_Collect(a.the_geom),ST_Collect(b.the_geom)) 
        LIMIT 1;

    UPDATE points2
    SET cluster_id = joined_clusters[1]
    WHERE cluster_id = joined_clusters[2];

    --If there's only 1 cluster left, exit loop
    IF (SELECT COUNT(DISTINCT cluster_id) FROM points2) < 2 THEN
        EXIT;

    END IF;

END LOOP;

RETURN QUERY SELECT gid, ST_TRANSFORM(the_geom, srid)::geometry(point), cluster_id FROM points2;
END;
$BODY$
LANGUAGE plpgsql

Использование:

WITH subq AS(
    SELECT ARRAY_AGG((gid, the_geom)::pt) AS points
    FROM data
    GROUP BY collection_id)
SELECT (clusters).* FROM 
    (SELECT buc(points, radius) AS clusters FROM subq
) y;
Рафаэль
источник