Ich suche nach einem räumlichen Clustering-Algorithmus, um ihn in einer PostGIS-fähigen Datenbank für Punkt-Features zu verwenden. Ich werde die Funktion plpgsql schreiben, die den Abstand zwischen Punkten innerhalb desselben Clusters als Eingabe misst. Bei der Ausgabe gibt die Funktion ein Array von Clustern zurück. Die naheliegendste Lösung besteht darin, um das Feature herum Pufferzonen zu erstellen und in diesem Puffer nach Features zu suchen. Wenn solche Features vorhanden sind, erstellen Sie weiterhin einen Puffer um sie herum usw. Wenn solche Features nicht vorhanden sind, bedeutet dies, dass die Clustererstellung abgeschlossen ist. Vielleicht gibt es einige clevere Lösungen?

Es gibt mindestens zwei gute Clustering-Methoden für PostGIS: k- Mittel (über kmeans-postgresqlErweiterung) oder Clustering-Geometrien innerhalb eines Schwellenabstands (PostGIS 2.2).

1) k bedeutet mitkmeans-postgresql

Installation: Auf einem POSIX-Hostsystem muss PostgreSQL 8.4 oder höher installiert sein (ich würde nicht wissen, wo ich mit MS Windows anfangen soll). Wenn Sie dies aus Paketen installiert haben, stellen Sie sicher, dass Sie auch die Entwicklungspakete haben (z. B. postgresql-develfür CentOS). Herunterladen und extrahieren:

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/

Vor dem Erstellen müssen Sie die USE_PGXS Umgebungsvariable festlegen (in meinem vorherigen Beitrag wurde angewiesen, diesen Teil der zu löschen Makefile, was nicht die beste Option war). Einer dieser beiden Befehle sollte für Ihre Unix-Shell funktionieren:

# bash
export USE_PGXS=1
# csh
setenv USE_PGXS 1

Erstellen und installieren Sie nun die Erweiterung:

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

(Hinweis: Ich habe dies auch mit Ubuntu 10.10 versucht, aber kein Glück, da der Pfad in pg_config --pgxsnicht existiert! Dies ist wahrscheinlich ein Ubuntu-Paketierungsfehler.)

Verwendung / Beispiel: Sie sollten irgendwo eine Punktetabelle haben (ich habe eine Reihe von Pseudozufallspunkten in QGIS gezeichnet). Hier ist ein Beispiel für das, was ich getan habe:

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;

Das 5im zweiten Argument der kmeansFensterfunktion angegebene I ist die K- Ganzzahl, um fünf Cluster zu erzeugen. Sie können dies in eine beliebige Ganzzahl ändern.

Unten sind die 31 Pseudozufallspunkte, die ich gezeichnet habe, und die fünf Zentroide mit der Bezeichnung, die die Anzahl in jedem Cluster angibt. Dies wurde mit der obigen SQL-Abfrage erstellt.

Sie können auch versuchen, mit ST_MinimumBoundingCircle zu veranschaulichen, wo sich diese Cluster befinden :

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) Clustering innerhalb eines Schwellenabstandes mit ST_ClusterWithin

Diese Aggregatfunktion ist in PostGIS 2.2 enthalten und gibt ein Array von GeometryCollections zurück, in dem sich alle Komponenten in einem Abstand voneinander befinden.

Hier ein Anwendungsbeispiel, bei dem ein Abstand von 100,0 der Schwellenwert ist, der zu 5 verschiedenen Clustern führt:

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;

Der größte mittlere Cluster hat einen umschließenden Kreisradius von 65,3 Einheiten oder ungefähr 130, was größer als der Schwellenwert ist. Dies liegt daran, dass die einzelnen Abstände zwischen den Elementgeometrien kleiner als der Schwellenwert sind, sodass sie zu einem größeren Cluster zusammengefasst werden.

Ich habe eine Funktion geschrieben, die Cluster von Features basierend auf dem Abstand zwischen ihnen berechnet und einen konvexen Rumpf über diesen Features erstellt:

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;

Beispiel für die Verwendung dieser Funktion:

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

'poi' - Name des Layers, 'wkb_geometry' - Name der Geometriespalte, 'ogc_fid' - Primärschlüssel der Tabelle, 14000 - Clusterabstand.

Das Ergebnis der Verwendung dieser Funktion:

Das bisher vielversprechendste, was ich gefunden habe, ist diese Erweiterung für K-Means-Clustering als Fensterfunktion: http://pgxn.org/dist/kmeans/

Ich konnte es jedoch noch nicht erfolgreich installieren.

Andernfalls können Sie für das grundlegende Grid-Clustering SnapToGrid verwenden .

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
;

Komplementäre Antwort von @MikeT ...

Für MS Windows:

Bedarf:

• Beliebige Visual C ++ Express-Version wie diese
• Das kmeans-postgresql- Modul.

Was du tun wirst:

• Optimieren Sie den Quellcode, um die kmeans-Funktion in eine DLL zu exportieren.
• Kompilieren Sie den Quellcode mit dem cl.exeCompiler, um eine DLL mit kmeansFunktion zu generieren .
• Legen Sie die generierte DLL im Ordner PostgreSQL \ lib ab.
• Anschließend können Sie die UDF über den SQL-Befehl in PostgreSQL "erstellen" (verknüpfen).

Schritte:

1. Anforderungen herunterladen und installieren / extrahieren.

2. Öffnen Sie das kmeans.cin einem beliebigen Editor:

   1. Nach den #includeZeilen definieren Sie das DLLEXPORT-Makro mit:

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

   2. Stellen Sie DLLEXPORTvor jede dieser Zeilen:

      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. Öffnen Sie die Visual C ++ - Befehlszeile.

4. In der Kommandozeile:

   1. Gehe zum Extrahierten kmeans-postgresql.
   2. Stellen Sie Ihren POSTGRESPATH ein, mein Beispiel ist: SET POSTGRESPATH=C:\Program Files\PostgreSQL\9.5

   3. Lauf

      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. Kopieren Sie das kmeans.dllan%POSTGRESPATH%\lib

6. Führen Sie nun den SQL-Befehl in Ihrer Datenbank aus, um die Funktion zu "ERSTELLEN".

   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;

Hier ist eine Möglichkeit, das Ergebnis der in 2) angegebenen PostGIS-Abfrage in dieser Antwort in QGIS anzuzeigen

Da QGIS weder Geometriesammlungen noch andere Datentypen in derselben Geometriespalte verarbeitet, habe ich zwei Ebenen erstellt, eine für Cluster und eine für Clusterpunkte.

Erstens brauchen Sie für Cluster nur Polygone, andere Ergebnisse sind Einsamkeitspunkte:

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'

Dann müssen Sie für Clusterpunkte Geometriesammlungen in Mehrpunkt transformieren:

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

Einige Punkte haben die gleichen Koordinaten, sodass die Beschriftung verwirrend sein kann.

Sie können die Kmeans-Lösung einfacher mit der ST_ClusterKMeans- Methode verwenden, die in postgis ab 2.3 verfügbar ist. Beispiel:

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;

Der Begrenzungsrahmen für Features wird im obigen Beispiel als Clustergeometrie verwendet. Das erste Bild zeigt die Originalgeometrien und das zweite ist das Ergebnis der obigen Auswahl.

Bottom-up-Clustering-Lösung von Holen Sie sich einen einzelnen Cluster aus Punktewolken mit maximalem Durchmesser in Postgis, für den keine dynamischen Abfragen erforderlich sind .

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

und ein Typ mit Cluster-ID

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

Als nächstes die Algorithmusfunktion

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

Verwendungszweck:

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;