2

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;

bounding box

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_contains

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_within


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_intersects

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);

st_within and st_centroid


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

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;
    
  • This ultimately lead me to what I wanted to achieve! Rather than cut/splice the geom like I thought I wanted to, I ended up selecting blocks with an area of more than 100 square feet in common with the zip in addition to the bounding box. You the man! – Tyler Sep 9 '14 at 19:05
0

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)

  • Thanks for the idea, I updated the question and added 2 more screenshots. It was helpful and I learned something but still not exactly what I'm looking for. – Tyler Aug 7 '14 at 18:10

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.