Pengelompokan spasial dengan PostGIS?

97

Saya mencari algoritme pengelompokan spasial untuk menggunakannya dalam database yang mendukung PostGIS untuk fitur titik. Saya akan menulis fungsi plpgsql yang membutuhkan jarak antara titik-titik dalam cluster yang sama sebagai input. Pada fungsi output mengembalikan array cluster. Solusi yang paling jelas adalah membangun zona buffer yang ditentukan jarak di sekitar fitur dan mencari fitur ke buffer ini. Jika fitur tersebut ada maka teruslah membangun buffer di sekitar mereka, dll. Jika fitur tersebut tidak ada itu berarti pembangunan cluster selesai. Mungkin ada beberapa solusi pintar?

drnextgis
sumber
4
Ada berbagai macam metode pengelompokan karena sifat data yang berbeda dan tujuan pengelompokan yang berbeda. Untuk tinjauan umum tentang apa yang ada di luar sana dan untuk beberapa bacaan mudah tentang apa yang dilakukan orang lain untuk mengelompokkan matriks jarak, cari situs CV @ SE . Bahkan, "memilih metode pengelompokan" hampir merupakan duplikat yang tepat dari Anda dan memiliki jawaban yang baik.
whuber
8
Memberi +1 pada pertanyaan karena menemukan contoh SQL PostGIS aktual alih-alih tautan ke algoritme adalah misi yang mustahil untuk apa pun selain pengelompokan kisi dasar, terutama untuk pengelompokan yang lebih eksotis seperti MCL
wildpeaks

Jawaban:

112

Setidaknya ada dua metode pengelompokan yang baik untuk PostGIS: k- berarti (melalui kmeans-postgresqlekstensi) atau geometri pengelompokan dalam jarak ambang batas (PostGIS 2.2)


1) k- artinya dengankmeans-postgresql

Instalasi: Anda harus memiliki PostgreSQL 8.4 atau lebih tinggi pada sistem host POSIX (saya tidak tahu harus mulai dari mana untuk MS Windows). Jika Anda menginstal ini dari paket, pastikan Anda juga memiliki paket pengembangan (misalnya, postgresql-develuntuk CentOS). Unduh dan ekstrak:

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/

Sebelum membangun, Anda perlu mengatur USE_PGXS variabel lingkungan (posting saya sebelumnya diinstruksikan untuk menghapus bagian ini Makefile, yang bukan pilihan terbaik). Salah satu dari dua perintah ini harus bekerja untuk shell Unix Anda:

# bash
export USE_PGXS=1
# csh
setenv USE_PGXS 1

Sekarang bangun dan instal ekstensi:

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

(Catatan: Saya juga mencoba ini dengan Ubuntu 10.10, tetapi tidak berhasil, karena jalurnya pg_config --pgxstidak ada! Ini mungkin merupakan bug pengemasan Ubuntu)

Penggunaan / Contoh: Anda harus memiliki tabel poin di suatu tempat (saya menggambar banyak poin pseudo acak di QGIS). Berikut ini adalah contoh dengan apa yang saya lakukan:

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;

yang 5saya disediakan dalam argumen kedua dari kmeansfungsi jendela adalah K integer untuk menghasilkan lima cluster. Anda bisa mengubahnya ke bilangan bulat apa pun yang Anda inginkan.

Di bawah ini adalah 31 poin acak semu yang saya gambar dan lima centroid dengan label yang menunjukkan jumlah di setiap cluster. Ini dibuat menggunakan query SQL di atas.

Kmeans


Anda juga dapat mencoba mengilustrasikan di mana cluster ini berada bersama 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) Clustering dalam jarak ambang batas dengan ST_ClusterWithin

Fungsi agregat ini disertakan dengan PostGIS 2.2, dan mengembalikan array GeometryCollections di mana semua komponen berada dalam jarak satu sama lain.

Berikut ini adalah contoh penggunaan, di mana jarak 100.0 adalah ambang batas yang menghasilkan 5 kelompok berbeda:

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

Cluster tengah terbesar memiliki jari-jari lingkaran melingkar dari 65,3 unit atau sekitar 130, yang lebih besar dari ambang batas. Ini karena jarak individual antara geometri anggota kurang dari ambang, sehingga mengikatnya bersama sebagai satu cluster yang lebih besar.

Mike T
sumber
2
Hebat, modifikasi ini akan membantu untuk instalasi :-) Namun saya khawatir saya tidak bisa menggunakan ekstensi itu pada akhirnya karena (jika saya mengerti dengan benar), diperlukan sejumlah cluster hardcoded magic, yang tidak masalah dengan data statis. Anda dapat menyetelnya terlebih dahulu tetapi tidak cocok untuk mengelompokkan set data yang sewenang-wenang (karena berbagai filter), misalnya kesenjangan besar dalam kluster 10-poin pada gambar terakhir. Namun ini akan membantu orang lain juga karena (afaik), ini adalah satu-satunya contoh SQL yang ada (kecuali satu liners di beranda ekstensi) untuk ekstensi itu.
wildpeaks
(ah Anda menjawab pada saat yang sama saya menghapus komentar sebelumnya untuk memformulasikannya, maaf)
wildpeaks
7
Untuk pengelompokan kmeans Anda perlu menentukan jumlah cluster di muka; Saya ingin tahu apakah ada algoritma alternatif di mana jumlah cluster tidak diperlukan.
djq
1
Versi 1.1.0 sekarang tersedia: api.pgxn.org/dist/kmeans/1.1.0/kmeans-1.1.0.zip
djq
1
@maxd no. Diberikan A = πr², lalu r = √ (A / π).
Mike T
27

Saya telah menulis fungsi yang menghitung kelompok fitur berdasarkan jarak di antara mereka dan membangun lambung cembung atas fitur ini:

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;

Contoh menggunakan fungsi ini:

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

'poi' - nama layer, 'wkb_geometry' - nama kolom geometri, 'ogc_fid' - kunci utama tabel, jarak 14000-cluster.

Hasil menggunakan fungsi ini:

masukkan deskripsi gambar di sini

drnextgis
sumber
Bagus! Bisakah Anda menambahkan contoh cara menggunakan fungsi Anda juga? Terima kasih!
underdark
1
Saya telah memodifikasi sedikit kode sumber dan telah menambahkan contoh penggunaan fungsi.
drnextgis
Baru saja mencoba menggunakan ini pada postgres 9.1 dan baris "FOR innr IN EXECUTE 'SELECT' || gid || ' SEBAGAI gid, '|| geom ||' SEBAGAI geom DARI tmp DI MANA dmn ADALAH NULL 'LOOP "menghasilkan kesalahan berikut. Ada ide? GALAT: fungsi set-dihargai yang dipanggil dalam konteks yang tidak dapat menerima set
bitbox
Saya tidak yakin bagaimana cara menggunakan kode ini dalam PG (PostGIS n00b) di tabel saya. di mana saya bisa mulai memahami sintaks ini? Saya punya meja dengan lat dan lons yang ingin saya klaster
mga
Pertama-tama Anda harus membuat geometrykolom di dalam tabel Anda, bukan untuk menyimpan lonlat secara terpisah dan membuat kolom dengan nilai unik (ID).
drnextgis
10

Sejauh ini, yang paling menjanjikan yang saya temukan adalah ekstensi untuk pengelompokan K-means sebagai fungsi jendela: http://pgxn.org/dist/kmeans/

Namun saya belum berhasil menginstalnya.


Jika tidak, untuk pengelompokan kisi dasar, Anda bisa menggunakan 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
sumber
2

Melengkapi jawaban @MikeT ...

Untuk MS Windows:

Persyaratan:

Apa yang akan kamu lakukan:

  • Tweak kode sumber untuk mengekspor fungsi kmeans ke DLL.
  • Kompilasi kode sumber dengan cl.exekompiler untuk menghasilkan DLL dengan kmeansfungsi.
  • Masukkan DLL yang dihasilkan ke folder PostgreSQL \ lib.
  • Kemudian Anda dapat "membuat" (tautan) UDF ke PostgreSQL melalui perintah SQL.

Langkah:

  1. Unduh & pasang / ekstrak persyaratan.
  2. Buka kmeans.cdi editor mana saja:

    1. Setelah #includebaris mendefinisikan makro DLLEXPORT dengan:

      #if defined(_WIN32)
          #define DLLEXPORT __declspec(dllexport)
      #else
         #define DLLEXPORT
      #endif
    2. Letakkan di DLLEXPORTdepan setiap baris ini:

      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. Buka Visual C ++ Command Line.

  4. Di baris perintah:

    1. Pergi ke yang diekstraksi kmeans-postgresql.
    2. Atur POSTGRESPATH Anda, milik saya misalnya adalah: SET POSTGRESPATH=C:\Program Files\PostgreSQL\9.5
    3. Lari

      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. Salin kmeans.dllke%POSTGRESPATH%\lib

  6. Sekarang jalankan perintah SQL di database Anda untuk "MENCIPTAKAN" fungsi.

    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
sumber
2

Berikut adalah cara untuk menampilkan di QGIS hasil dari permintaan PostGIS yang diberikan dalam 2) di server ini

Karena QGIS tidak menangani pengumpulan geometri atau tipe data yang berbeda dalam kolom geometri yang sama, saya telah membuat dua lapisan, satu untuk cluster dan satu untuk poin-poin yang dikelompokkan.

Pertama untuk cluster, Anda hanya perlu poligon, hasil lainnya adalah poin kesepian:

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'

Kemudian untuk titik berkerumun, Anda perlu mengubah koleksi geometri dalam multipoint:

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

Beberapa titik berada pada koordinat yang sama sehingga label dapat membingungkan.

Pengelompokan dalam QGIS

Nicolas Boisteault
sumber
2

Anda dapat menggunakan solusi Kmeans lebih mudah dengan metode ST_ClusterKMeans yang tersedia di postgis dari 2.3 Contoh:

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;

Kotak pembatas fitur digunakan sebagai geometri klaster pada contoh di atas. Gambar pertama menunjukkan geometri asli dan yang kedua adalah hasil pilih di atas.

Geometri asli Cluster fitur

Volda
sumber
1

Solusi pengelompokan bawah ke atas dari Dapatkan satu kluster dari awan titik dengan diameter maksimum pada postgis yang tidak melibatkan kueri dinamis.

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

dan tipe dengan id cluster

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

Selanjutnya fungsi algoritma

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

Pemakaian:

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;
raphael
sumber