I am having a table of polygons with some 15 Mio. records. The polygons hold an attribute called "dn" with integer values between 0 and 6. My goal is to simplify that table so that polygons below a certain size get merged with neighboring ones. First step would be to merge dn=0 polygons of 100 square meters or below with neighboring polygons with dn=1.

The polygons are simple and angular shaped. Resulting from gdal_polygonize.py a raster representation of the data.

My naive attempt was to :

create table temp_union as  select 1 as dn, st_union(a_geom,b_geom) from
   (select geom as a_geom from big_poly_table where dn=1) as foo,
   (select geom as b_geom from big_poly_table where dn=0 AND st_area(geom)<=100) as bar
   WHERE st_touches(a_geom,b_geom);

This for testing purposes on a subset of some 200000 members.

Two problems I am running into:

  • The query runs some 3 hours on a very decent server. Speeding up is important.
  • The resulting table temp_union seems to get really huge as the respective table dump is five times bigger than the original table.

Is this to be expected?

select postgis_full_version();
POSTGIS="2.0.3 r11128" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.1, released 2013/08/2
6" LIBXML="2.7.6" RASTER
(1 Zeile)

EDIT: As pointed out in my comment to @MickyT his solution solved my problem. A minor change is referring to the attribute filter on the bigger polygons the smaller ones should be merged to. They have to hold an specific value for the atrribute "dn".

Another - for me - important finding is that the test on "connectedness" of the polygons should not be performed with "ST_Touches" but with "ST_Intersects" as this runs by an order of magnitude faster. In my case both yield the same results.

In this blogpost by Paul Ramsey on "Performance: Prepared Geometry", "ST_Intersects" is in the list of PostGis functions that take advantage of that, "ST_Touches" is not in the list. So here the slightly revamped SQL script that is now in production use:

 SELECT ogc_fid, geom,
  (SELECT ogc_fid from w001001_1_1 l
   WHERE dn = 1 and
    ST_Intersects(s.geom, l.geom)
    ORDER  BY l.the_area DESC
    LIMIT 1 ) as d_id
into temp_dissolvelist 
 FROM w001001_1_1 s
WHERE dn = 0 and the_area <= 100;
  • The reason that your temp_union is larger than the original, is that you are most likely doing a many to many join. For example a large polygon that has 12 small polygons touching it will be returned 12 times unioned with each of the individual smaller polygons. Those smaller polygons could also be unioned with another larger polygon in turn. – MickyT Dec 6 '15 at 20:51
  • @Micky, that sounds quite on the point. There should be a filter introduced (on the fids, I would think) . Small polygons should be discarded from further unioning once the have been merged – waterfront Dec 6 '15 at 22:12

Here's something that may help get you started and is based on a couple of assumptions that I have made.
Unfortunately I don't have a lot of time that I can spend on this at the moment, so I haven't tested this. I try to mock up some test data for testing later.

The general idea with this process is:

  • Foreach small polygon (dn = 0 and area <= 100) find a single large polygon to merge with. In my example they a ranked by dn increasing and area decreasing.
  • Foreach large polygon referred to by the previous step, aggregate and union with all the small polygons referencing it.

I have assumed that there is some IDs in place.

/* Create a set of small polygons with the larger polygon it should dissolve with */
SELECT id, geom,
    (SELECT id 
    FROM big_poly_table l 
    WHERE dn <> 0 and 
        ST_Touches(s.geom, l.geom) 
    ORDER BY dn, ST_Area(l.geom) DESC 
    LIMIT 1) d_id
INTO temp_dissolvelist
FROM big_poly_table s
WHERE dn = 0 and ST_Area(geom) <= 100

/* Create a dissolved table */
SELECT id, MAX(dn) dn, ST_Union(geom) geom
INTO temp_union
    /* large polygons to use in union */
    SELECT id, dn, geom
    FROM big_poly_table a
    WHERE exists (SELECT 1 FROM temp_dissolvelist b WHERE b.d_id = a.id)
    SELECT d_id, 0, geom
    FROM temp_dissolvelist
    ) d
  • Thanks, MickyT, for taking the time. Part one of your solution takes about 1/3 of the processing time with my test table, part two uses a few seconds. The result exactly is what's intended. I added a "the_area" column (int) to the table beforehand, I would think that may speed up the area comparison. Question: do I really need to put the geometry (geom) into table "temp_dissolvelist"? Couldn't I just have pairs of small-poly-fids for dn=0, touching poly-fids for dn=1 and union their geoms from the original big_poly_table? Would it make any difference in terms of processing time? – waterfront Dec 7 '15 at 18:05
  • 1
    @waterfront you don't really need to put the geometry into the temp_dissolvelist. I put it there to keep the query as simple and clear as I could, also I was a bit lazy. In terms of processing time, it may have an impact, but I wouldn't expect it to be too large. – MickyT Dec 7 '15 at 18:15

Are you sure you have built spatial index beforehand?

CREATE INDEX ON big_poly_table USING gist(geom);

I'm not very sure how PostGIS deals with subqueries you use, hence I propose you use a simple JOIN clause.

CREATE TABLE temp_union AS
    ST_Union(a.geom, b.geom)
FROM big_poly_table a
JOIN big_poly_table b ON ST_Touches(a.geom ,b.geom)
WHERE a.dn = 1
    AND b.dn = 0
    AND ST_Area(b.geom) <= 100;
  • Thanks @Michal. Yes, the index is there. Your suggested JOIN is running as I write. – waterfront Dec 6 '15 at 22:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.