How to replicate ArcGIS Intersect in PostGIS

Question

I'm trying to replicate this ArcGIS process in PostGIS: http://blogs.esri.com/esri/arcgis/2012/11/13/spaghetti_and_meatballs/. This describes how to break buffered points down into polygons based on their intersections, counting the number of layers, and attributing that to the polygons in order to classify them. I'm using it to create a rough point density map with vectors, and the results were surprisingly nice for my data set in ArcGIS. However, I am struggling to come up with something workable in PostGIS where I need it for producing dynamic point density layers for a web map.

In ArcGIS, I simply ran the Intersect tool on my buffered points layer to create the shapes I needed.

In PostGIS, I ran this query:

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;

The output looks pretty much identical to the ArcGIS output, except that it is not breaking the polygons down to the same extent that is required for a meaningful density map. Here are screenshots of what I mean:

ArcGIS is on the left, and PostGIS is on the right. It is slightly difficult to tell, but the ArcGIS image shows the 'interior' polygon created where all 3 buffers intersect. The PostGIS output, on the other hand, does not create that interior polygon and instead it keeps its components intact. This makes it impossible to provide a classification for just that interior area with 3 layers on top of each other compared to just 1 for the outer parts.

Does anyone know of any PostGIS function to break the polygon down to the extent I need? Alternatively, does anyone know of a better way to produce a point density map with vectors in PostGIS?

Answer 1

You can do this all in one step by chaining the CTEs together, but I did it in several so I could look at the results in QGIS as I progressed.

First, generate a bunch of random points to work with, using a gaussian distribution so we get more overlap in the middle.

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;

Now buffer the points into circles so we get some overlap.

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

Now, extract just the boundaries from the circles. If you have polygons with holes, you'll have to use ST_DumpRings() and get more fancy here. I have simple polygons so I cheat. Once you have the boundaries, union them against themselves (actually any small piece of coincident linework will do) to force them to be noded and deduplicated. (This is magic.)

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

Now rebuild areas using the noded linework. This is the broken down areas, with only one polygon per area. After polygonizing, dump the individual polygons out of the multipolygon output.

create sequence polyseq;

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

Now, add a place for the polygon count and fill it up by joining the centroids of the small cut-up polygons to the original circles, and summarizing for each small piece. For larger data sets an index on the circles table at least will be required to make things not impossibly slow.

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;

That's it, you now have no overlapping polygons, but each resultant polygon has a count on it that says how many overlaps it is standing in for.

Answer 2

The method I ended up using was to create a fishnet grid in my area of interest with a high enough "resolution" to style and reflect the data to reasonable degree. You can read about the fishnet function here: How to create a regular polygon grid in PostGIS?

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

This creates the fishnet with 800 rows, 850 columns, that are 0.0005 radians in height and length (using WGS84 projection in lat/long and its a small enough geographic extent that the distortion is negligible - ie they're all distorted more or less equally), and then the coordinates for the bottom left of the grid.

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

Because this created a huge amount of polygons which will have queries ran on them, I created an index and updated the statistics. This reduced my typical queries from 50+ seconds to 4-5 seconds.

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;

The subquery here counts the number of incidents within .002 radians (approx 220 meters) of each fishnet grid polygon, and groups them by the fishnet grid. This effectively counts the number of overlapping circles to the resolution of the grid.

The outer query I used to Union each polygon's count value, and restrict the count to 3 or greater. While the union isn't strictly necessary and is the most resource intensive part of the query, it is critical for web mapping as it effectively turns tens of thousands of grid polygons, which doesn't work too well when serving directly to openlayers, into multipolygons of however many different count values there are (usually a few dozen for my data).

Restricting the count value is an important ability for heat maps so they don't depict too much data to the point of being unable to interpret it - it also has the added utility of speeding the query up significantly.

Final result: