4

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

California State Leg Lower

Upper Chamber Districts

California State Leg Upper

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

QGIS Intersections

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.

partial intersections

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

6
  • ST_Intersects works in both directions, so ST_Intersects(a.geom, b.geom) will give same answers as ST_Intersects(b.geom, a.geom). You often get data that claims to be a multipolygon, but whish is in fact a polygon, so you can probably ignore that. If you run ST_NumGeometries you will get the answer to that. If greater than 1, use ST_Dump to get just the polygons. Jul 16, 2018 at 14:17
  • So, Select sll.name, slu.name, ST_Intersections(sll.geom, slu.geom) FROM state_leg_lower_single sll, state_leg_upper_single slu WHERE ST_Intersects(sll.geom, slu.geom) should be all you need. Jul 16, 2018 at 14:19
  • @JohnPowellakaBarça Hi, That was the first version of the query that I ran, but I got the same results. Does this work for you? Jul 16, 2018 at 18:24
  • @ThingumaBob : I can’t figure out how to add a bounty to this but if you show me that query, I will happily hire you to write it. Jul 17, 2018 at 2:24
  • I deleted all my traces here for now...my line of thoughts was plain wrong ,)
    – geozelot
    Jul 18, 2018 at 10:16

2 Answers 2

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.

2
  • How do I mark this as correct? Aug 28, 2018 at 18:28
  • 1
    according to this, there should be a green checkmark to the left of my answer
    – yorkie
    Aug 28, 2018 at 19:06
0

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.

California legislative overlaps image

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.