Como replicar o ArcGIS Intersect no PostGIS

Estou tentando replicar esse processo do ArcGIS no PostGIS: http://blogs.esri.com/esri/arcgis/2012/11/13/spaghetti_and_meatballs/ . Isso descreve como dividir os pontos em buffer em polígonos com base em suas interseções, contando o número de camadas e atribuindo isso aos polígonos para classificá-los. Estou usando-o para criar um mapa de densidade de pontos aproximados com vetores, e os resultados foram surpreendentemente bons para meus dados no ArcGIS. No entanto, estou lutando para encontrar algo viável no PostGIS, onde preciso para produzir camadas dinâmicas de densidade de pontos para um mapa da web.

No ArcGIS, eu simplesmente executei a ferramenta Intersect na minha camada de pontos em buffer para criar as formas necessárias.

No PostGIS, executei 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;

A saída parece praticamente idêntica à saída do ArcGIS, exceto que não está dividindo os polígonos na mesma extensão que é necessária para um mapa de densidade significativo. Aqui estão as capturas de tela do que quero dizer:

O ArcGIS está à esquerda e o PostGIS está à direita. É um pouco difícil dizer, mas a imagem do ArcGIS mostra o polígono 'interior' criado onde os três buffers se cruzam. A saída do PostGIS, por outro lado, não cria esse polígono interno e mantém seus componentes intactos. Isso torna impossível fornecer uma classificação apenas para a área interna com 3 camadas uma sobre a outra, em comparação com apenas 1 para as partes externas.

Alguém sabe de alguma função PostGIS para quebrar o polígono na medida em que eu preciso? Como alternativa, alguém conhece uma maneira melhor de produzir um mapa de densidade de pontos com vetores no PostGIS?

Você pode fazer tudo isso de uma só vez, encadeando as CTEs, mas eu fiz isso em várias para que eu pudesse ver os resultados no QGIS à medida que progredia.

Primeiro, gere um monte de pontos aleatórios para trabalhar, usando uma distribuição gaussiana, para obter mais sobreposições no meio.

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;

Agora, coloque os pontos em círculos, para que se sobreponham.

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

Agora, extraia apenas os limites dos círculos. Se você tiver polígonos com orifícios, precisará usar ST_DumpRings () e obter mais sofisticação aqui. Eu tenho polígonos simples, então eu trapaceio. Depois de ter os limites, una-os contra si mesmos (na verdade, qualquer pequeno pedaço de linha coincidente servirá) para forçá-los a serem acenados e deduplicados. (Isto é mágico.)

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

Agora reconstrua as áreas usando o trabalho de linha assentido. Essas são as áreas divididas, com apenas um polígono por área. Após a poligonalização, despeje os polígonos individuais da saída do multipolígono.

create sequence polyseq;

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

Agora, adicione um local para a contagem de polígonos e preencha juntando os centróides dos pequenos polígonos de recorte aos círculos originais e resumindo cada pedaço pequeno. Para conjuntos de dados maiores, é necessário pelo menos um índice na tabela de círculos para tornar as coisas não impossivelmente lentas.

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;

É isso aí, agora você não tem polígonos sobrepostos, mas cada polígono resultante tem uma contagem que indica quantas sobreposições ele está representando.

O método que acabei usando foi criar uma grade de rede de pesca na minha área de interesse com uma "resolução" alta o suficiente para estilizar e refletir os dados em um grau razoável. Você pode ler sobre a função arrastão aqui: Como criar uma grade regular de polígonos no PostGIS?

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

Isso cria a rede de pesca com 800 linhas e 850 colunas, com 0,0005 radianos de altura e comprimento (usando a projeção WGS84 em lat / long e é uma extensão geográfica suficientemente pequena para que a distorção seja desprezível - ou seja, todas elas são distorcidas mais ou menos igualmente ) e as coordenadas da parte inferior esquerda da grade.

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

Como isso criou uma enorme quantidade de polígonos que terão consultas executadas neles, criei um índice e atualizei as estatísticas. Isso reduziu minhas consultas típicas de mais 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;

A subconsulta aqui conta o número de incidentes dentro de 0,002 radianos (aproximadamente 220 metros) de cada polígono da grade da rede de pesca e os agrupa pela grade da rede de pesca. Isso conta efetivamente o número de círculos sobrepostos à resolução da grade.

A consulta externa que usei para unir o valor de contagem de cada polígono e restringir a contagem a 3 ou mais. Embora a união não seja estritamente necessária e seja a parte da consulta que consome mais recursos, ela é essencial para o mapeamento da Web, pois efetivamente transforma dezenas de milhares de polígonos de grade, o que não funciona muito bem ao servir diretamente a camadas abertas, em multipolígonos com quantos valores de contagem diferentes existem (geralmente algumas dúzias para meus dados).

Restringir o valor da contagem é uma habilidade importante para os mapas de calor, para que eles não representem muitos dados a ponto de serem incapazes de interpretá-los - ele também possui o utilitário adicional de acelerar significativamente a consulta.

Resultado final: