Finding combinations of polygons in PostGIS

Question

all,

I have found something that has made me a bit curious about how to correctly use ST_INTERSECTION(A,B). I am looking for every combination of two political district boundary sets in California, the state senate and the state assembly. These are two different sets of districts in California, both of which cover the entire state. There are fewer polygons in the upper chamber (senate) than in the lower chamber (assembly), which means that each upper polygon will be bigger and intersect with multiple lower polygons, possibly containing them.

The shapefiles that I am using are from the US Census, but I have preserved the specific state of the files that I am using in data.world if you want to try them.

The first thing to note is that loading these into PSQL via shp2pgsql results in them coming in as MULTIPOLYGON, so I recast them to POLYGON using these queries:

DROP TABLE IF EXISTS state_leg_lower_single;
CREATE TABLE state_leg_lower_single AS
  (
    SELECT
      sll.name,
      sll.statefp,
      sll.sldlst,
      sll.lsad,
      sll.lsy,
      (ST_DUMP(sll.geom)).geom :: geometry(Polygon, 4326) AS geom
    FROM
      state_legislative_lower sll
  );

DROP TABLE IF EXISTS state_leg_upper_single;
CREATE TABLE state_leg_upper_single AS
  (
    SELECT
      slu.name,
      slu.statefp,
      slu.sldust,
      slu.lsad,
      slu.lsy,
      (ST_DUMP(slu.geom)).geom :: geometry(Polygon, 4326) AS geom
    FROM
      state_legislative_upper slu
  );

Both of these seem to come through correctly in QGIS.

Lower Chamber Districts

Upper Chamber Districts

Now, here is what Quantum GIS does when I use Vector -> Geoprocessing -> Intersection, which is what I expected.

To get this result in PostGIS, I have tried so many variations on the following query:

DROP TABLE IF EXISTS state_legislative_intersections;
CREATE TABLE state_legislative_intersections AS
  (

    SELECT
      'SLL-' || sll.name || ': SLU-' || slu.name AS composing_districts,
      sll.name                                   AS state_leg_lower_name,
      slu.name                                   AS state_leg_upper_name,
      COALESCE(
          ST_INTERSECTION(sll.geom, slu.geom),
          ST_INTERSECTION(slu.geom, sll.geom)
      )                                          as intersecting_geom
    FROM state_leg_lower_single sll, state_leg_upper_single slu
    WHERE
      (ST_INTERSECTS(sll.geom, slu.geom) OR ST_INTERSECTS(slu.geom, sll.geom))
      AND
      GeometryType(st_intersection(sll.geom, slu.geom)) IN ('POLYGON', 'MULTIPOLYGON')

  );

What I get is a subset of what I imagine should be the result, only 71 polygons, with massive lacunae.

This has become a matter of curiosity to me, so I would like to see what I am doing incorrectly. Can anyone help?

Answer 1

I ignored the conversion from polygons to multipolygons, and your data seems to be in EPSG:4269 rather than EPSG:4326. Having a close look at the documentation of ST_Intersection gives an example on how to handle these situations:

---Clip all lines (trails) by country (here we assume country geom are POLYGON or MULTIPOLYGONS)
-- NOTE: we are only keeping intersections that result in a LINESTRING or MULTILINESTRING because we don't
-- care about trails that just share a point
-- the dump is needed to expand a geometry collection into individual single MULT* parts
-- the below is fairly generic and will work for polys, etc. by just changing the where clause
SELECT clipped.gid, clipped.f_name, clipped_geom
FROM (SELECT trails.gid, trails.f_name, (ST_Dump(ST_Intersection(country.the_geom, trails.the_geom))).geom As clipped_geom
FROM country
    INNER JOIN trails
    ON ST_Intersects(country.the_geom, trails.the_geom))  As clipped
    WHERE ST_Dimension(clipped.clipped_geom) = 1 ;

--For polys e.g. polygon landmarks, you can also use the sometimes faster hack that buffering anything by 0.0
-- except a polygon results in an empty geometry collection
--(so a geometry collection containing polys, lines and points)
-- buffered by 0.0 would only leave the polygons and dissolve the collection shell
SELECT poly.gid,  ST_Multi(ST_Buffer(
                ST_Intersection(country.the_geom, poly.the_geom),
                0.0)
                ) As clipped_geom
FROM country
    INNER JOIN poly
    ON ST_Intersects(country.the_geom, poly.the_geom)
    WHERE Not ST_IsEmpty(ST_Buffer(ST_Intersection(country.the_geom, poly.the_geom),0.0));

Meaning, in your context the query would look something like this:

create table state_legislative_intersections as(
select
    'SLL-' || sll.name || ': SLU-' || slu.name AS composing_districts,
    sll.name                                   AS state_leg_lower_name,
    slu.name                                   AS state_leg_upper_name,
    st_multi(st_buffer(st_intersection(slu.geom,sll.geom),0.0)) as geom
from state_leg_lower_single as sll inner join state_leg_upper_single as slu 
    on (st_intersects(sll.geom,slu.geom) )
where not st_isempty(st_buffer(st_intersection(slu.geom,sll.geom),0.0)));

So, st_buffer combined with the where clause removes points and lines, and st_multi transforms it to multipolygons, I suspect just to be sure not to have more than one line with the same content in the result. This query gave me 198 multipolygons that does not seem to overlap and look equal to the screenshot of the QGIS intersection you posted.

Answer 2

I think your only problem was the filter on POLYGON and MULTIPOLYGON since there were probably some geometry collections that got returned. I dumped those, then applied the filter, and then regrouped. I tested against my 2020 copy of Tiger data.

WITH a AS ( 
    SELECT 
          'SLL-' || sll.namelsad  || ': SLU-' || slu.namelsad AS composing_districts,
          sll.namelsad                                   AS state_leg_lower_name,
          slu.namelsad                                   AS state_leg_upper_name,
              (st_dump(ST_INTERSECTION(sll.geom, slu.geom))).* 
    FROM tiger2020.sldl sll 
    JOIN tiger2020.sldu slu ON st_intersects(sll.geom, slu.geom)
    WHERE sll.statefp ='06' 
) 
SELECT composing_districts, state_leg_lower_name, state_leg_upper_name 
    , st_union(geom) AS intersecting_geom 
FROM a 
WHERE GeometryType(geom) IN ('POLYGON', 'MULTIPOLYGON')
GROUP BY 1,2,3 

Here's the result that I get.