¿Clúster espacial con PostGIS?

91

Estoy buscando un algoritmo de agrupamiento espacial para usarlo dentro de la base de datos habilitada con PostGIS para las características de los puntos. Voy a escribir la función plpgsql que toma distancia entre los puntos dentro del mismo grupo que la entrada. En la función de salida devuelve la matriz de los clústeres. La solución más obvia es crear una distancia especificada de zonas de búfer alrededor de la característica y buscar características en este búfer. Si existen tales características, continúe construyendo un búfer alrededor de ellas, etc. Si no existen tales características, eso significa que se completó la construcción del clúster. Tal vez hay algunas soluciones inteligentes?

    
pregunta drnextgis 28.06.2011 - 17:39

6 respuestas

102

Hay al menos dos buenos métodos de agrupación para PostGIS: k -means (a través de la extensión kmeans-postgresql ) o geometrías de agrupación dentro de una distancia de umbral (PostGIS 2.2)

1) k -significa con kmeans-postgresql

Instalación: necesita tener PostgreSQL 8.4 o superior en un sistema host POSIX (no sabría por dónde empezar para MS Windows). Si tiene esto instalado desde los paquetes, asegúrese de tener también los paquetes de desarrollo (por ejemplo, postgresql-devel para CentOS). Descargar y extraer:

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/

Antes de construir, debes configurar la variable de entorno USE_PGXS (mi publicación anterior me indicó que eliminara esta parte de el Makefile , que no era la mejor de las opciones). Uno de estos dos comandos debería funcionar para su shell de Unix:

# bash
export USE_PGXS=1
# csh
setenv USE_PGXS 1

Ahora compile e instale la extensión:

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

(Nota: También probé esto con Ubuntu 10.10, pero no tuve suerte, ya que la ruta en pg_config --pgxs no existe. Probablemente sea un error de empaquetado de Ubuntu)

Uso / Ejemplo: Deberías tener una tabla de puntos en algún lugar (dibujé un montón de puntos pseudoaleatorios en QGIS). Aquí hay un ejemplo con lo que hice:

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;

el 5 que proporcioné en el segundo argumento de la función de ventana kmeans es el entero K para producir cinco grupos. Puedes cambiar esto a cualquier entero que quieras.

Debajo están los 31 puntos pseudoaleatorios que dibujé y los cinco centroides con la etiqueta que muestra el conteo en cada grupo. Esto se creó utilizando la consulta SQL anterior.

Tambiénpuedeintentarilustrardóndeseencuentranestosgruposcon 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)AgrupacióndentrodeunadistanciadeumbralconST_ClusterWithin

Esta función agregada se incluye con PostGIS 2.2 y devuelve una variedad de GeometryCollections donde todos los componentes están a una distancia el uno del otro.

Este es un ejemplo de uso, donde una distancia de 100.0 es el umbral que da como resultado 5 grupos diferentes:

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;

El grupo central más grande tiene un radio de círculo circundante de 65.3 unidades o aproximadamente 130, que es más grande que el umbral. Esto se debe a que las distancias individuales entre las geometrías de los miembros son menores que el umbral, por lo que las une como un grupo más grande.

    
respondido por el Mike T 04.07.2011 - 11:36
26

He escrito una función que calcula grupos de características basadas en la distancia entre ellas y construye un casco convexo sobre estas características:

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;

Ejemplo de uso de esta función:

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

'poi' - nombre de la capa, 'wkb_geometry' - nombre de la columna de geometría, 'ogc_fid' - clave principal de la tabla, 14000 - distancia de agrupación.

El resultado de usar esta función:

    
respondido por el drnextgis 05.07.2011 - 06:40
10

Hasta ahora, lo más prometedor que encontré es esta extensión para K-means clustering como una función de ventana: enlace

Sin embargo, aún no he podido instalarlo correctamente.

De lo contrario, para la agrupación en cuadrícula básica, podría usar 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
;
    
respondido por el wildpeaks 29.06.2011 - 20:58
2

Complementando la respuesta de @MikeT ...

Para MS Windows:

Requisitos:

Lo que harás:

  • Ajuste el código fuente para exportar la función kmeans a una DLL.
  • Compile el código fuente con el compilador cl.exe para generar una DLL con la función kmeans .
  • Coloque la DLL generada en la carpeta PostgreSQL \ lib.
  • Luego, puede "crear" (vincular) la UDF en PostgreSQL a través del comando SQL.

Pasos:

  1. descargar & requisitos de instalación / extracción.
  2. Abra el kmeans.c en cualquier editor:

    1. Después de #include , las líneas definen la macro DLLEXPORT con:

      #if defined(_WIN32)
          #define DLLEXPORT __declspec(dllexport)
      #else
         #define DLLEXPORT
      #endif
      
    2. Coloque DLLEXPORT antes de cada una de estas líneas:

      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. Abre la línea de comandos de Visual C ++.

  4. En la línea de comando:

    1. Vaya al kmeans-postgresql extraído.
    2. Establezca su POSTGRESPATH, el mío, por ejemplo, es: SET POSTGRESPATH=C:\Program Files\PostgreSQL.5
    3. Ejecutar

      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. Copie el kmeans.dll a %POSTGRESPATH%\lib

  6. Ahora ejecute el comando SQL en su base de datos para "CREAR" la función.

    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;
    
respondido por el caiohamamura 23.01.2016 - 01:09
2

Esta es una forma de mostrar en QGIS el resultado de la consulta de PostGIS dada en 2) en esta respuesta

Como QGIS no maneja ni las colecciones de geometría ni los tipos de datos diferentes en la misma columna de geometría, he creado dos capas, una para grupos y otra para puntos agrupados.

Primero para los grupos, solo necesitas polígonos, otros resultados son puntos solitarios:

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'

Luego, para los puntos agrupados, necesita transformar las colecciones de geometría en multipunto:

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

Algunos puntos están en las mismas coordenadas, por lo que la etiqueta podría ser confusa.

    
respondido por el Nicolas Boisteault 03.03.2017 - 16:30
1

Solución de agrupamiento ascendente desde Obtenga un solo clúster desde una nube de puntos con diámetro máximo en postgis que no implica consultas dinámicas.

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

y un tipo con ID de clúster

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

A continuación la función de algoritmo

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

Uso:

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;
    
respondido por el raphael 30.09.2014 - 21:22

Lea otras preguntas en las etiquetas