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

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

1) k- significa conkmeans-postgresql

Instalación: debe tener PostgreSQL 8.4 o superior en un sistema host POSIX (no sabría dónde comenzar para MS Windows). Si tiene esto instalado desde paquetes, asegúrese de tener también los paquetes de desarrollo (por ejemplo, postgresql-develpara 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 compilar, debe establecer la USE_PGXS variable de entorno (mi publicación anterior me indicó que eliminara esta parte de 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 compila e instala 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 pg_config --pgxsno existe! Probablemente sea un error de empaquetado de Ubuntu)

Uso / Ejemplo: debe tener una tabla de puntos en alguna parte (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 5que proporcioné en el segundo argumento de la kmeansfunción de ventana es el entero K para producir cinco grupos. Puede cambiar esto a cualquier número entero que desee.

A continuación se muestran los 31 puntos seudoaleatorios que dibujé y los cinco centroides con la etiqueta que muestra el conteo en cada grupo. Esto fue creado usando la consulta SQL anterior.

También puede intentar ilustrar dónde están estos grupos con 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ón dentro de una distancia umbral con ST_ClusterWithin

Esta función agregada se incluye con PostGIS 2.2 y devuelve una matriz de GeometryCollections donde todos los componentes están a una distancia entre sí.

Aquí hay 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 envolvente 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 lo une como un grupo más grande.

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 primaria de la tabla, 14000 - distancia del clúster.

El resultado de usar esta función:

Hasta ahora, lo más prometedor que encontré es esta extensión para la agrupación de K-means como una función de ventana: http://pgxn.org/dist/kmeans/

Sin embargo, aún no he podido instalarlo con éxito.

De lo contrario, para la agrupación en cuadrícula básica, puede 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
;

Complementando la respuesta @MikeT ...

Para MS Windows:

Requisitos:

• Cualquier versión de Visual C ++ Express como esta
• El módulo kmeans-postgresql .

Qué harás:

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

Pasos:

1. Descargar e instalar / extraer requisitos.

2. Abra el kmeans.cen cualquier editor:

   1. Después de las #includelíneas, defina la macro DLLEXPORT con:

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

   2. Poner DLLEXPORTantes 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. Abra la línea de comandos de Visual C ++.

4. En la linea de comando:

   1. Ir a lo extraído kmeans-postgresql.
   2. Establezca su POSTGRESPATH, el mío, por ejemplo, es: SET POSTGRESPATH=C:\Program Files\PostgreSQL\9.5

   3. correr

      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.dlla%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;

Aquí hay una manera de mostrar en QGIS el resultado de la consulta PostGIS dada en 2) en esta respuesta

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

Primero para los clústeres, solo necesita 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, debe 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 puede ser confusa.

Puede usar la solución Kmeans más fácilmente con el método ST_ClusterKMeans que está disponible en postgis desde 2.3 Ejemplo:

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;

El cuadro delimitador de características se utiliza como geometría de clúster en el ejemplo anterior. La primera imagen muestra las geometrías originales y la segunda es el resultado de seleccionar arriba.

Solución de agrupamiento ascendente de Get a single cluster desde la nube de puntos con un 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)

Luego la función del 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;