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