Selecting Census Block Groups inside Zip code using PostGIS?

Question

I already have my data imported, problem is my SELECT statements aren't yielding the results I'm looking for. I specifically need all the Block Group multipolygons that are inside a Zip Code multipolygon so that all of the zip code is covered and no Block Group leaves the boundaries of the zip code.

Blue = Zip Code

Red = Block Group

&& -- bounding box

select blockgroups.the_geom 
from blockgroups, zips 
where zips.zip = '92656' AND zips.the_geom && blockgroups.the_geom;

ST_Contains

select blockgroups.the_geom 
from blockgroups, zips 
where zips.zip = '92656' AND zips.the_geom && blockgroups.the_geom 
    AND ST_Contains(zips.the_geom,blockgroups.the_geom);

ST_Within

select blockgroups.the_geom 
from blockgroups, zips 
where zips.zip = '92656' AND blockgroups.the_geom && zips.the_geom 
    AND ST_Within(blockgroups.the_geom, zips.the_geom);

---

ST_Intersects

select blockgroups.the_geom 
from blockgroups, zips 
where zips.zip = '92656' AND blockgroups.the_geom && zips.the_geom 
    AND ST_Intersects(blockgroups.the_geom,zips.the_geom);

ST_Within + ST_Centroid

select blockgroups.the_geom 
from blockgroups, zips 
where zips.zip = '92656' AND blockgroups.the_geom && zips.the_geom 
    AND ST_Within(ST_Centroid(blockgroups.the_geom),zips.the_geom);

---

If block group multipolygon does leave the boundary of the zip code, I would like it to be cut by the zip code boundary, resulting in a fully covered zip code by block groups.

---

@mapBaker had a good idea but still not exactly what I want. Using the centroid preserves the block group shapes, but still doesn't cover the entire zip boundary. Added screenshots of ST_Intersects and ST_Within + ST_Centroid

Answer 1

You're really most of the way there. While && or ST_Intersects will find you which block groups overlap which ZIP codes, you are returning the geometry unmodified in your SELECT list. If you want the block groups clipped to the ZIP codes, you have to use a function which actually modifies the geometry, viz. ST_Intersection. In the PostGIS Reference this is the difference between the "Spatial Relationships and Measurements" functions, which generally return a number or a boolean (e.g., does it touch, how far away is it?) and the "Geometry Processing" functions which actually return a new geometry.

SELECT ST_Multi(ST_Intersection(b.the_geom, z.the_geom)) AS geom
FROM blockgroups b, zips z
WHERE z.zip = '92656' AND ST_Intersects(b.the_geom, z.the_geom);

or, as I prefer to write these queries:

SELECT ST_Multi(ST_Intersection(b.the_geom, z.the_geom)) AS geom
FROM blockgroups b JOIN zips z ON (ST_Intersects(b.the_geom, z.the_geom))
WHERE z.zip = '92656';

Some notes:

1. ST_Multi is necessary to coerce the POLYGONs into MULTIPOLYGONs or your resultset will include a mix. Some GIS applications will choke on this.

2. You don't need both the bounding box test and ST_Intersects. From the docs:

   This function call will automatically include a bounding box comparison that will make use of any indexes that are available on the geometries.

   This applies to pretty much all of the spatial relationship functions (ST_Overlaps, ST_Touches, etc.).

3. You may run into a problem with slivers if your data are from different sources. You can eliminate them in the query by filtering out geometries with areas below a reasonable threshold. In data I was testing that were in a foot-based projection, I eliminated all resulting geometries below 100 square feet:

   SELECT ST_Multi(ST_Intersection(b.the_geom, z.the_geom) AS geom
   FROM blockgroups b JOIN zips z ON (ST_Intersects(b.the_geom, z.the_geom))
   WHERE z.zip = '92656'
       AND ST_Area(ST_Intersection(b.the_geom, z.the_geom)) < 100;
   

Answer 2

Using SQL Server, I've found luck converting the target layer to points (centroids) then selecting them:

select * from polygonA, polygonB
where
    polygonA.shape.STIntersects.STCentroid().shape(polygonB.shape) = 1)

and here's the syntax I use for STWithin, which also does a great job:

select * from polygonA, polygonB
where
    polygonA.shape.STWithin(polygonB.shape) = 1)

The syntax for PostGIS isn't much different...

(I posted this as an answer on this thread)