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?

  • 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. – John Powell Jul 16 '18 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. – John Powell Jul 16 '18 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? – Dheeraj Chand Jul 16 '18 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. – Dheeraj Chand Jul 17 '18 at 2:24
  • I deleted all my traces here for now...my line of thoughts was plain wrong ,) – ThingumaBob Jul 18 '18 at 10:16
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.

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

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.