Replicating"Intersection" function in QGIS using PostGIS queries

Question

I have two multipolygon shapefiles that I've imported into my PostGIS database - one of congressional districts (cds) and one of US counties (ctys). Below shows the counties in blue and the districts overlain with orange border, no fill.

INPUT DATA

I would like to perform a simple intersection using a PostGIS/PostgreSQL query that would produce the same output as the "Intersection" tool in QGIS, which creates a multipolygon output that consists of polygons for all intersected areas between the two files, as well as all of the attributes of each file. This desired output is shown here:

DESIRED OUTPUT (QGIS/Intersection)

I tried to achieve this in PostGIS using ST_Intersection() along with ST_Intersect() like this:

CREATE TABLE test AS
SELECT ST_Intersection(a.geom, b.geom) AS geom
FROM cds AS a, ctys AS b
WHERE ST_Intersects(a.geom, b.geom);

but the resulting output is a messy, multiple geometry-type output with polygons, lines, points and gaps in seemingly random places.

BAD OUTPUT (PostGIS/ST_Intersection())

The geometries of each input layer (checked with the QGIS/Check Validity function) appear to be valid.

What am I doing wrong?

Answer 1

I don't think you have a processing problem, I think you probably have a visualization problem, related to the fact that your output will be quite diverse in terms of geometry type, while your viewer might expect only to get one geometry type. So filter out things, like so.

CREATE TABLE test AS
SELECT ST_Multi(ST_CollectionExtract(ST_Intersection(a.geom, b.geom),3))::geometry(MultiPolygon) AS geom
FROM cds AS a, ctys AS b
WHERE ST_Intersects(a.geom, b.geom);

If the processing was erroring out, the whole query would have stopped and you've have no output table.