1

I am using pgAdmin (v.1.22) for my postgres(9.5). I am still fairly new to postgres. I am trying to perform a union(intersection) of multiple layers where there is both complete and no overlap between layers. I need all linework and attributes to come through to the final product. I am playing around with st_union and st_intersection but can't quite get to the final layer I need.

DROP TABLE IF EXISTS TEST;
CREATE TABLE TEST AS 
SELECT
 (st_dump(st_union(a.shape, b.shape))).geom geom, a.status, b.code
FROM 
 p760_passive_lb a, p760_fmu b 
WHERE
 st_intersects(a.shape, b.shape)
;
 ALTER TABLE TEST ADD COLUMN id SERIAL;
ALTER TABLE TEST ADD PRIMARY KEY (id);
 CREATE INDEX idx_TEST_geom ON TEST USING gist(geom);

Where am I going wrong?

This is similar to: Is there a union function for multiple layers comparable to ArcGIS in open source GIS like QGIS?

So, to help illustrate, I have attached an image of what is going on. The top two layers are my inputs. The bottom left is the output of the above code. Two things of note: (1) the features in the western portion are not visible, and (2) when I use the identify tool (bottom right) in QGIS on any ONE area in the result, the whole area gets chosen and multiple features are selected. Strangely, the linework in the western section is viewable with the identify tool. This leads me to think that the layers are just stacked on top of each other and not forming one geometry.

enter image description here

  • How is it not working? – John Powell Feb 5 '16 at 5:47
  • 2
    You say you want the layers merged even when there is no overlap -- so having ST_Intersects(a.shape, b.shape) means you only merge when there is an overlap. Is this your problem? – John Powell Feb 5 '16 at 6:57
  • If I do not have the intersect I do not get any linework shared between the boundaries (bottom left image). – D_C Feb 5 '16 at 17:13
1

I found the answer to be a two part solution. Get the parts that overlap and then add the pieces that don't:

CREATE TABLE TEST AS 
SELECT fix_geometry(CASE WHEN l.id IS NULL THEN f.geom ELSE st_intersection(f.geom, l.geom) END) geom, f.FIELD, l.FIELD
FROM X f 
LEFT JOIN Y l ON st_intersects(l.geom, f.geom);

INSERT INTO TEST (geom, f_FIELD, L_FIELD)
SELECT fix_geometry(st_difference(st_union(f.geom), st_union(l.geom))) geom, f.FIELD, l.FIELD
FROM X f
JOIN Y l ON st_intersects(l.geom, f.geom)
GROUP BY f.objectid, l.FIELD;
ALTER TABLE TEST ADD COLUMN id SERIAL;
ALTER TABLE TEST ADD PRIMARY KEY (id);
CREATE INDEX idx_TEST_geom ON TEST USING gist(geom);

Fix_geometry is a function that makes the returning geometry only be allowed to be polygonal. it buffers by 0 makes valid, and snaps to a small grid

  • What is fix_geometry? – jpmc26 Jan 18 at 18:39

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.