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;