1

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

enter image description here

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)

enter image description here

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())

enter image description here

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

What am I doing wrong?

  • 1
    You're going to encounter this over and over again if you're moving from QGIS to POSTGIS. PostGIS is much more powerful, but also much more picky that QGIS. When you run an intersection in POSTGIS, it's going to return the intersection points, lines, and polygons in a geometry collection. Look at the documentation for st_collectionextract - you can extract just the polygons from your output. – jbalk Sep 20 at 21:52
  • 2
    Also read the documentation for st_intersection and st_intersects. It will help to understand what is being returned by your query. – jbalk Sep 20 at 21:54
  • Understood about the st_intersection return formats etc. Is there not a way to achieve a similar output to the QGIS function then? – deepsky Sep 20 at 22:38
  • 1
    Could it be that in QGIS the layer is defined to show only polygons (and not multipolygons)? – JGH Sep 21 at 1:32
  • 2
    Could you reduce your example to 1 or 2 polygons in each layer in order to reproduce and simplify the problem and then provide their wkt representation? I made a similar example here and it behave as expected. There must be something with your polygons... – Pierre Racine Sep 23 at 14:56
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.

  • worked like a charm - thank you so much! – deepsky Sep 30 at 6:19

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.