Szukam algorytmu klastrowania przestrzennego do używania go w bazie danych obsługującej PostGIS dla obiektów punktowych. Zamierzam napisać funkcję plpgsql, która pobiera odległość między punktami w tym samym klastrze co dane wejściowe. Na wyjściu funkcja zwraca tablicę klastrów. Najbardziej oczywistym rozwiązaniem jest zbudowanie określonych stref buforowych wokół obiektu i poszukiwanie funkcji do tego bufora. Jeśli takie funkcje istnieją, nadal buduj wokół nich bufor itp. Jeśli takie funkcje nie istnieją, oznacza to, że tworzenie klastrów zostało zakończone. Może są jakieś sprytne rozwiązania?

Istnieją co najmniej dwie dobre metody grupowania dla PostGIS: k- oznacza (poprzez kmeans-postgresqlrozszerzenie) lub geometrie grupowania w odległości progowej (PostGIS 2.2)

1) k- oznacza zkmeans-postgresql

Instalacja: Musisz mieć PostgreSQL 8.4 lub nowszy na systemie hosta POSIX (nie wiedziałbym od czego zacząć dla MS Windows). Jeśli masz to zainstalowane z pakietów, upewnij się, że masz również pakiety programistyczne (np. postgresql-develDla CentOS). Pobierz i rozpakuj:

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/

Przed budowaniem musisz ustawić USE_PGXS zmienną środowiskową (mój poprzedni post polecił usunąć tę część Makefile, która nie była najlepsza z opcji). Jedno z tych dwóch poleceń powinno działać dla twojej powłoki Unix:

# bash
export USE_PGXS=1
# csh
setenv USE_PGXS 1

Teraz skompiluj i zainstaluj rozszerzenie:

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

(Uwaga: próbowałem tego również z Ubuntu 10.10, ale bez powodzenia, ponieważ ścieżka pg_config --pgxsnie istnieje! To prawdopodobnie błąd pakowania Ubuntu)

Zastosowanie / przykład: Powinieneś gdzieś mieć tabelę punktów (narysowałem kilka pseudolosowych punktów w QGIS). Oto przykład tego, co zrobiłem:

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, że znajduje się w drugim argumencie kmeansfunkcji okna jest K całkowitą produkować pięciu klastrów. Możesz zmienić to na dowolną liczbę całkowitą, którą chcesz.

Poniżej znajduje się 31 pseudolosowych punktów, które narysowałem i pięć centroidów z etykietą pokazującą liczbę w każdej grupie. Zostało to utworzone przy użyciu powyższego zapytania SQL.

Możesz także spróbować zilustrować położenie tych klastrów za pomocą 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;

2) Grupowanie w odległości progowej za pomocą ST_ClusterWithin

Ta funkcja agregująca jest zawarta w PostGIS 2.2 i zwraca tablicę GeometryCollections, w której wszystkie komponenty znajdują się w pewnej odległości od siebie.

Oto przykład zastosowania, w którym odległość 100,0 to próg, który daje 5 różnych klastrów:

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;

Największa środkowa gromada ma promień otaczającego okręgu wynoszący 65,3 jednostki lub około 130, czyli więcej niż próg. Wynika to z faktu, że indywidualne odległości między geometriami prętów są mniejsze niż próg, więc wiąże je razem jako jedna większa gromada.

Napisałem funkcję, która oblicza skupienia cech na podstawie odległości między nimi i buduje wypukły kadłub nad tymi cechami:

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;

Przykład użycia tej funkcji:

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

„poi” - nazwa warstwy, „wkb_geometry” - nazwa kolumny geometrii, „ogc_fid” - klucz podstawowy tabeli, 14000 - odległość klastra.

Wynik użycia tej funkcji:

Jak dotąd najbardziej obiecujące jest to rozszerzenie dla K-oznacza klastrowanie jako funkcja okna: http://pgxn.org/dist/kmeans/

Nie udało mi się jednak zainstalować go pomyślnie.

W przeciwnym razie do podstawowego grupowania siatki można użyć 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
;

Uzupełnienie odpowiedzi @MikeT ...

W przypadku MS Windows:

Wymagania:

• Każda wersja Visual C ++ Express, taka jak ta
• Moduł kmeans-postgresql .

Co będziesz robić:

• Popraw kod źródłowy, aby wyeksportować funkcję kmeans do biblioteki DLL.
• Skompiluj kod źródłowy za pomocą cl.exekompilatora, aby wygenerować bibliotekę DLL z kmeansfunkcją.
• Umieść wygenerowaną bibliotekę DLL w folderze PostgreSQL \ lib.
• Następnie możesz „utworzyć” (link) UDF do PostgreSQL za pomocą polecenia SQL.

Kroki:

1. Pobierz i zainstaluj / rozpakuj wymagania.

2. Otwórz kmeans.cw dowolnym edytorze:

   1. Po #includewierszach zdefiniuj makro DLLEXPORT za pomocą:

      #if defined(_WIN32)
          #define DLLEXPORT __declspec(dllexport)
      #else
         #define DLLEXPORT
      #endif

   2. Umieść DLLEXPORTprzed każdą z tych linii:

      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. Otwórz wiersz poleceń Visual C ++.

4. W wierszu poleceń:

   1. Idź do wyodrębnionego kmeans-postgresql.
   2. Ustaw swoją POSTGRESPATH, moja na przykład to: SET POSTGRESPATH=C:\Program Files\PostgreSQL\9.5

   3. Biegać

      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. Skopiuj kmeans.dlldo%POSTGRESPATH%\lib

6. Teraz uruchom polecenie SQL w bazie danych, aby „UTWÓRZ” funkcję.

   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;

Oto sposób wyświetlenia w QGIS wyniku zapytania PostGIS podanego w 2) w tej odpowiedzi

Ponieważ QGIS nie obsługuje ani kolekcji geometrycznych, ani różnych typów danych w tej samej kolumnie geometrii, stworzyłem dwie warstwy, jedną dla klastrów i jedną dla punktów skupionych.

Po pierwsze dla klastrów potrzebujesz tylko wielokątów, inne wyniki to samotne punkty:

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'

Następnie dla punktów klastrowych musisz przekształcić geometrycollection w wielopunktowe:

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

Niektóre punkty mają te same współrzędne, więc etykieta może być myląca.

Możesz łatwiej używać rozwiązania Kmeans dzięki metodzie ST_ClusterKMeans dostępnej w Postgis od 2.3 Przykład:

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;

Obwiednia elementów jest używana jako geometria skupienia w powyższym przykładzie. Pierwszy obraz pokazuje oryginalne geometrie, a drugi jest wynikiem wyboru powyżej.

Oddolne rozwiązanie klastrowe od Uzyskaj pojedynczy klaster z chmury punktów o maksymalnej średnicy w postgis, który nie wymaga dynamicznych zapytań.

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

i typ z identyfikatorem klastra

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

Następnie funkcja algorytmu

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

Stosowanie:

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;