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