So replizieren Sie ArcGIS Intersect in PostGIS

Ich versuche, diesen ArcGIS-Prozess in PostGIS zu replizieren: http://blogs.esri.com/esri/arcgis/2012/11/13/spaghetti_and_meatballs/ . Hier wird beschrieben, wie gepufferte Punkte anhand ihrer Schnittpunkte in Polygone zerlegt werden, die Anzahl der Ebenen gezählt und diese den Polygonen zugeordnet werden, um sie zu klassifizieren. Ich verwende es, um eine grobe Punktdichtekarte mit Vektoren zu erstellen, und die Ergebnisse waren für meinen Datensatz in ArcGIS überraschend gut. Ich habe jedoch Probleme, in PostGIS etwas zu finden, das für die Erstellung dynamischer Punktdichteebenen für eine Webkarte erforderlich ist.

In ArcGIS habe ich einfach das Schnittwerkzeug auf meiner Ebene mit gepufferten Punkten ausgeführt, um die benötigten Formen zu erstellen.

In PostGIS habe ich diese Abfrage ausgeführt:

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;

Die Ausgabe sieht ziemlich identisch mit der ArcGIS-Ausgabe aus, außer dass die Polygone nicht in dem Maße zerlegt werden, wie es für eine aussagekräftige Dichtekarte erforderlich ist. Hier sind Screenshots von dem, was ich meine:

ArcGIS befindet sich links und PostGIS rechts. Es ist etwas schwer zu sagen, aber das ArcGIS-Bild zeigt das 'innere' Polygon, das dort erstellt wurde, wo sich alle 3 Puffer schneiden. Die PostGIS-Ausgabe hingegen erstellt dieses innere Polygon nicht und behält stattdessen seine Komponenten bei. Dies macht es unmöglich, eine Klassifizierung für nur diesen Innenbereich mit 3 übereinanderliegenden Schichten im Vergleich zu nur 1 für die Außenteile bereitzustellen.

Kennt jemand eine PostGIS-Funktion, um das Polygon so weit aufzuschlüsseln, wie ich es brauche? Kennt jemand einen besseren Weg, um eine Punktdichtekarte mit Vektoren in PostGIS zu erstellen?

Sie können dies alles in einem Schritt tun, indem Sie die CTEs miteinander verketten. Ich habe es jedoch in mehreren Schritten getan, damit ich mir die Ergebnisse in QGIS im weiteren Verlauf ansehen kann.

Generieren Sie zunächst eine Reihe von zufälligen Punkten, mit denen Sie arbeiten möchten, und verwenden Sie dabei eine Gaußsche Verteilung, damit sich in der Mitte mehr Überlappungen ergeben.

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;

Puffern Sie nun die Punkte in Kreise, damit wir eine gewisse Überlappung erhalten.

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

Extrahieren Sie nun nur die Grenzen aus den Kreisen. Wenn Sie Polygone mit Löchern haben, müssen Sie ST_DumpRings () verwenden und hier mehr Fantasie bekommen. Ich habe einfache Polygone, also betrüge ich. Sobald Sie die Grenzen haben, vereinigen Sie sie gegen sich selbst (tatsächlich reicht jedes kleine Stück zufälliger Linienführung aus), um sie zum Nicken und Deduplizieren zu zwingen. (Das ist Magie.)

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

Bauen Sie nun Bereiche mit den geknoteten Linien wieder auf. Dies sind die aufgeschlüsselten Bereiche mit nur einem Polygon pro Bereich. Speichern Sie nach der Polygonisierung die einzelnen Polygone aus der Multipolygonausgabe.

create sequence polyseq;

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

Fügen Sie nun einen Platz für die Polygonzahl hinzu und füllen Sie ihn aus, indem Sie die Schwerpunkte der kleinen zerschnittenen Polygone mit den ursprünglichen Kreisen verbinden und für jedes kleine Stück zusammenfassen. Für größere Datenmengen ist mindestens ein Index für die Kreistabelle erforderlich, um die Dinge nicht unmöglich langsam zu machen.

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;

Das war's, Sie haben jetzt keine überlappenden Polygone, aber jedes resultierende Polygon hat eine Zählung, die angibt, für wie viele Überlappungen es steht.

Die Methode, die ich letztendlich verwendete, bestand darin, in meinem Interessengebiet ein Netzgitter mit einer ausreichend hohen "Auflösung" zu erstellen, um die Daten in angemessenem Maße zu formatieren und wiederzugeben. Informationen zur Fischnetzfunktion finden Sie hier: Wie erstelle ich ein reguläres Polygonraster in PostGIS?

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

Dadurch entsteht das Fischnetz mit 800 Zeilen und 850 Spalten mit einer Höhe und Länge von 0,0005 Radiant (unter Verwendung der WGS84-Projektion in Lat / Long und einer ausreichend kleinen geografischen Ausdehnung, dass die Verzerrung vernachlässigbar ist - dh alle sind mehr oder weniger gleichmäßig verzerrt ) und dann die Koordinaten unten links im Raster.

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

Da dadurch eine große Anzahl von Polygonen erstellt wurde, auf denen Abfragen ausgeführt werden, habe ich einen Index erstellt und die Statistiken aktualisiert. Dies reduzierte meine typischen Abfragen von 50+ Sekunden auf 4-5 Sekunden.

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;

Die Unterabfrage zählt hier die Anzahl der Vorfälle innerhalb von 0,002 Radian (ca. 220 Meter) jedes Fischnetzgitterpolygons und gruppiert sie nach dem Netzgitter. Dies zählt effektiv die Anzahl überlappender Kreise zur Auflösung des Gitters.

Die äußere Abfrage, mit der ich den Zählwert jedes Polygons vereinigt und die Anzahl auf 3 oder mehr beschränkt habe. Während die Vereinigung nicht unbedingt erforderlich ist und den ressourcenintensivsten Teil der Abfrage darstellt, ist sie für das Web-Mapping von entscheidender Bedeutung, da sie Zehntausende von Rasterpolygonen effektiv umwandelt, was nicht direkt funktioniert, wenn sie direkt für Openlayer bereitgestellt werden Multipolygone mit so vielen verschiedenen Zählwerten (normalerweise ein paar Dutzend für meine Daten).

Das Einschränken des Zählwerts ist eine wichtige Funktion für Heatmaps, damit sie nicht zu viele Daten darstellen, bis sie nicht mehr interpretiert werden können. Außerdem bietet es den zusätzlichen Nutzen, die Abfrage erheblich zu beschleunigen.

Endergebnis: