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.
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?