Comment répliquer ArcGIS Intersect dans PostGIS

J'essaie de reproduire ce processus ArcGIS dans PostGIS: http://blogs.esri.com/esri/arcgis/2012/11/13/spaghetti_and_meatballs/ . Il décrit comment diviser les points tamponnés en polygones en fonction de leurs intersections, en comptant le nombre de couches et en l'attribuant aux polygones afin de les classer. Je l'utilise pour créer une carte de densité de points approximative avec des vecteurs, et les résultats étaient étonnamment bons pour mon ensemble de données dans ArcGIS. Cependant, j'ai du mal à trouver quelque chose de réalisable dans PostGIS où j'en ai besoin pour produire des couches de densité de points dynamiques pour une carte Web.

Dans ArcGIS, j'ai simplement exécuté l'outil Intersection sur ma couche de points tamponnés pour créer les formes dont j'avais besoin.

Dans PostGIS, j'ai exécuté cette requête:

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 sortie semble à peu près identique à la sortie ArcGIS, sauf qu'elle ne décompose pas les polygones dans la même mesure que celle requise pour une carte de densité significative. Voici des captures d'écran de ce que je veux dire:

ArcGIS est à gauche et PostGIS à droite. C'est un peu difficile à dire, mais l'image ArcGIS montre le polygone «intérieur» créé à l'intersection des 3 tampons. En revanche, la sortie PostGIS ne crée pas ce polygone intérieur et conserve ses composants intacts. Il est donc impossible de fournir une classification uniquement pour cette zone intérieure avec 3 couches les unes sur les autres, contre 1 pour les parties extérieures.

Quelqu'un connaît-il une fonction PostGIS pour décomposer le polygone autant que je le souhaite? Sinon, quelqu'un connaît-il une meilleure façon de produire une carte de densité de points avec des vecteurs dans PostGIS?

Vous pouvez faire tout cela en une seule étape en enchaînant les CTE ensemble, mais je l'ai fait en plusieurs afin de pouvoir consulter les résultats dans QGIS au fur et à mesure de ma progression.

Tout d'abord, générez un tas de points aléatoires avec lesquels travailler, en utilisant une distribution gaussienne afin d'obtenir plus de chevauchement au milieu.

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;

Maintenant, tamponnez les points dans des cercles pour obtenir un certain chevauchement.

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

Maintenant, extrayez uniquement les limites des cercles. Si vous avez des polygones avec des trous, vous devrez utiliser ST_DumpRings () et obtenir plus de fantaisie ici. J'ai des polygones simples donc je triche. Une fois que vous avez les limites, réunissez-les contre elles-mêmes (en fait, n'importe quel petit morceau de dessin coïncident fera l'affaire) pour les forcer à hocher la tête et à dédupliquer. (C'est magique.)

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

Reconstruisez maintenant les zones en utilisant le dessin au trait incliné. Il s'agit des zones décomposées, avec un seul polygone par zone. Après la polygonisation, videz les polygones individuels de la sortie multipolygone.

create sequence polyseq;

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

Maintenant, ajoutez une place pour le nombre de polygones et remplissez-la en joignant les centroïdes des petits polygones découpés aux cercles d'origine et en résumant pour chaque petit morceau. Pour les ensembles de données plus volumineux, un index sur la table des cercles au moins sera nécessaire pour rendre les choses pas incroyablement lentes.

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;

C'est tout, vous n'avez plus de polygones qui se chevauchent, mais chaque polygone résultant a un décompte qui indique le nombre de chevauchements pour lesquels il représente.

La méthode que j'ai finalement utilisée était de créer une grille en résille dans ma zone d'intérêt avec une "résolution" suffisamment élevée pour styliser et refléter les données à un degré raisonnable. Vous pouvez lire sur la fonction résille ici: Comment créer une grille polygonale régulière dans PostGIS?

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

Cela crée le résille avec 800 lignes, 850 colonnes, qui sont de 0,0005 radians en hauteur et en longueur (en utilisant la projection WGS84 en lat / long et sa suffisamment petite étendue géographique pour que la distorsion soit négligeable - c'est-à-dire qu'elles sont toutes déformées plus ou moins également ), puis les coordonnées en bas à gauche de la grille.

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

Parce que cela a créé une énorme quantité de polygones sur lesquels des requêtes seront exécutées, j'ai créé un index et mis à jour les statistiques. Cela a réduit mes requêtes typiques de 50+ secondes à 4-5 secondes.

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 sous-requête compte ici le nombre d'incidents à moins de 0,002 radians (environ 220 mètres) de chaque polygone de la grille de résille, et les regroupe par la grille de résille. Cela compte effectivement le nombre de cercles qui se chevauchent jusqu'à la résolution de la grille.

La requête externe que j'ai utilisée pour réunir la valeur de comptage de chaque polygone et restreindre le nombre à 3 ou plus. Bien que l'union ne soit pas strictement nécessaire et constitue la partie la plus gourmande en ressources de la requête, elle est essentielle pour la cartographie Web car elle transforme efficacement des dizaines de milliers de polygones de grille, ce qui ne fonctionne pas trop bien en servant directement aux couches ouvertes, en multipolygones de toutes sortes de valeurs de comptage différentes (généralement quelques dizaines pour mes données).

La restriction de la valeur de comptage est une capacité importante pour les cartes de chaleur afin qu'elles ne représentent pas trop de données au point de ne pas pouvoir l'interpréter - elle a également l'utilité supplémentaire d'accélérer considérablement la requête.

Résultat final: