Cómo replicar ArcGIS intersección en PostGIS

Question

Estoy tratando de replicar este proceso ArcGIS en PostGIS: http://blogs.esri.com/esri/arcgis/2012/11/13/spaghetti_and_meatballs/. Esto describe cómo romper búfer de puntos en polígonos basándose en sus intersecciones, contando el número de capas, y la atribución de que a los polígonos con el fin de clasificar. Lo estoy utilizando para crear un punto áspero mapa de densidad con los vectores, y los resultados fueron sorprendentemente bueno para mi conjunto de datos en ArcGIS. Sin embargo, estoy luchando para llegar a algo viable en PostGIS que me falta para la producción de la dinámica de la densidad de puntos capas de un mapa web.

En ArcGIS, yo simplemente corrió la herramienta Intersecar en mi búfer puntos de capa para crear las formas que yo necesitaba.

En PostGIS, me encontré con esta consulta:

CREATE TABLE buffer_table AS SELECT a.gid AS gid, ST_Buffer(a.geo,.003) AS geo FROM public.pointTable a;

CREATE TABLE intersections AS SELECT a.gid AS gid_a, b.gid AS gid_b, ST_Intersection(a.geo,b.geo) AS geo FROM public.pointTable a, public.pointTable b WHERE ST_Intersects(a.geo, b.geo) AND a.gid < b.gid;

DELETE FROM intersections WHERE id_a = id_b;

La salida se ve casi idéntica a la de ArcGIS salida, excepto que no es romper los polígonos abajo en la misma medida que se requiere de una significativa mapa de densidad. Aquí están las capturas de pantalla de lo que quiero decir:

ArcGIS es en la izquierda, y PostGIS está a la derecha. Es un poco difícil de decir, pero el ArcGIS imagen muestra el interior del polígono creado donde todos 3 los búferes se cruzan. El PostGIS salida, por otro lado, no cree que el interior del polígono y en su lugar mantiene sus componentes intactos. Esto hace que sea imposible proporcionar una clasificación sólo para esa zona interior con 3 capas en la parte superior de uno al otro en comparación a sólo el 1 por la parte exterior.

¿Alguien sabe de alguna PostGIS función para romper el polígono hasta el punto de que necesito? Alternativamente, ¿alguien sabe de una mejor manera para producir una densidad de puntos de mapa con vectores en PostGIS?

Answer 1

Usted puede hacer todo esto en un solo paso por el encadenamiento de las Cte juntos, pero yo lo hice en varias para que yo pudiera mirar los resultados en QGIS a medida que avanzaba.

En primer lugar, generar un montón de puntos aleatorios para trabajar, utilizando una distribución de gauss, de modo que más se superponen en el medio.

create table pts as with 
    rands as (
        select generate_series as id, random() as u1, random() as u2 
        from generate_series(1,100))
select
      id,
      st_setsrid(st_makepoint(
            50 * sqrt(-2 * ln(u1)) * cos(2*pi()*u2), 
            50 * sqrt(-2 * ln(u1)) * sin(2*pi()*u2)),4326) as geom
from rands;

Ahora búfer de los puntos en círculos así que tenemos algunas coincidencias.

create table circles as
    select id, st_buffer(geom, 10) as geom from pts;

Ahora, la extracción de los límites de los círculos. Si usted tiene los polígonos con agujeros, tendrás que utilizar ST_DumpRings() y obtener más elegante, aquí. He polígonos simples así que me estafan. Una vez que usted tiene de los límites, de la unión de ellos en contra de ellos (en realidad, cualquier pequeño trozo de coincidencia de la línea de trabajo va a hacer) para obligarlos a ser noded y desduplicados. (Esta es la magia.)

create table boundaries as
    select st_union(st_exteriorring(geom)) as geom from circles;

Ahora la reconstrucción de áreas utilizando el noded de tendido de redes. Este es el roto áreas, con sólo un polígono por la zona. Después de polygonizing, el volcado de los polígonos individuales de la multipolígono de salida.

create sequence polyseq;

create table polys as
    select 
        nextval('polyseq') as id, 
        (st_dump(st_polygonize(geom))).geom as geom 
    from boundaries;

Ahora, agregue un lugar para el conteo de polígonos y la llenan de uniéndose a los centroides de los pequeños trozos de polígonos a los círculos originales, y resumiendo para cada pieza pequeña. Para grandes conjuntos de datos de un índice en los círculos de la tabla, al menos, estarán obligados a hacer cosas que no extremadamente lento.

create index circles_gix on circles using gist (geom);

alter table polys add column count integer default 0;

update polys set count = p.count 
from (
    select count(*) as count, 
           p.id as id 
    from polys p 
    join circles c 
    on st_contains(c.geom, st_pointonsurface(p.geom)) 
    group by p.id
) as p
where p.id = polys.id;

Eso es todo, que ahora no tienen la superposición de polígonos, pero cada polígono resultante tiene un recuento sobre lo que dice de cómo muchos se superpone que está de pie en la.

Answer 2

El método terminé usando era crear una red de cuadrícula en mi área de interés con un nivel suficientemente alto de "resolución" del estilo y reflejan los datos a un grado razonable. Usted puede leer acerca de la rejilla de la función aquí: Cómo crear un polígono regular de la cuadrícula en PostGIS?

CREATE TABLE fishnet AS
SELECT * FROM ST_CreateFishnet(800,850,.0005,.0005,-104.9190,38.7588);

Esto crea la red con 800 filas, 850 columnas, que son 0.0005 radianes en altura y en longitud (utilizando WGS84 proyección en lat/long y es lo suficientemente pequeña extensión geográfica que la distorsión es insignificante - es decir, todos están distorsionados más o menos igual) y, a continuación, las coordenadas de la esquina inferior izquierda de la cuadrícula.

UPDATE fishnet SET geom = ST_SetSRID(geom,4326);
CREATE INDEX fishnet_geom ON fishnet USING gist (geom);
ANALYZE fishnet;

Porque esto ha creado una gran cantidad de polígonos que se han consultas corrió en ellos, he creado un índice y actualizado de las estadísticas. Este reducido mis consultas más típicas de 50 segundos para 4-5 segundos.

SELECT ST_Union(a.geom), a.count
FROM (SELECT count(*) as count, fishnet.geom as geom
    FROM fishnet, incidents
    WHERE ST_DWithin(incidents.geo,fishnet.geom,.002) AND (incidents.incidenttype = 'Burglary')
    GROUP BY fishnet.geom) a
WHERE a.count >= 3
GROUP BY a.count;

La subconsulta aquí cuenta el número de incidentes dentro de .002 radianes (aprox 220 metros) de cada red de cuadrícula poligonal, y los grupos por la red grid. De esta forma, se cuenta el número de círculos que se superponen a la resolución de la cuadrícula.

La consulta externa he utilizado para la Unión de cada polígono del valor de recuento, y de restringir el conteo de 3 o mayor. Mientras que la unión no es estrictamente necesario y es el más intensivo en recursos, parte de la consulta, es fundamental para el mapeo web como se convierte a decenas de miles de cuadrícula de polígonos, que no funciona demasiado bien al momento de servir directamente a openlayers, en multipolígonos de sin embargo muchos diferentes valores de recuento de hay (por lo general un par de docenas de mis datos).

La restricción de que el valor de recuento es una habilidad importante para los mapas de calor para que no representan demasiado de datos hasta el punto de ser incapaz de interpretar - también tiene la mayor utilidad de la aceleración de la consulta de manera significativa.

Resultado Final: