I am trying to intersect two MultiPolygon tables, which are in the same UTM projection, and calculate the area of their intersections. I am getting some surprising results when I examine the intersections, which I retain in a new table. Specifically, each table contains either 2000 or 2010 Census block group boundaries. Here, I find the intersections with a
WHERE clause, calculate the weight associated with the area of that intersection, and put the results into a new table.
SELECT c2000.bkgpidfp00 AS fips2000, c2010.geoid10 AS fips2010, ST_Area(ST_Intersection(c2000.geom, c2010.geom)) / ST_Area(c2000.geom) AS weight, ST_Intersection(c2010.geom, c2000.geom) AS geom INTO xwalk FROM census2000_block_groups c2000, census2010_block_groups c2010 WHERE ST_Intersects(c2000.geom, c2010.geom) AND NOT ST_Touches(c2000.geom, c2010.geom);
Problem is, there are gaps in the final result. Not all of the intersections are present! See images below.
Here is the intersection--gaps are in white.
Here are the 2000 (blue) and 2010 (red) block group boundaries--you can see there are other intersections that should have been found!
Finally, just to prove that these intersections make sense, I used QGIS to calculate the intersections as well--and it worked! **Below is the correct intersection from QGIS ("Vector" > "Geoprocessing" > "Intersect").
Any ideas why my PostGIS
ST_Intersection example is failing to produce all of the possible intersections?
In the interest of reproducibility, here's a shell script to download and subset the data from the U.S. Census:
wget ftp://ftp2.census.gov/geo/pvs/tiger2010st/26_Michigan/26/tl_2010_26_bg00.zip wget ftp://ftp2.census.gov/geo/pvs/tiger2010st/26_Michigan/26/tl_2010_26_bg10.zip unzip ~/Downloads/tl_2010_26_bg00.zip ogr2ogr -f "ESRI Shapefile" -where "COUNTYFP00 IN ('099', '125', '163')" -s_srs "EPSG:4269" -t_srs "EPSG:32617" census_block_groups_2000_utm17n.shp tl_2010_26_bg00.shp unzip ~/Downloads/tl_2010_26_bg10.zip -d $WD ogr2ogr -f "ESRI Shapefile" -where "COUNTYFP10 IN ('099', '125', '163')" -s_srs "EPSG:4269" -t_srs "EPSG:32617" $WD/census_block_groups_2010_utm17n.shp tl_2010_26_bg10.shp
To insert into PostGIS:
shp2pgsql census_block_groups_2000_utm17n.shp public.census2000_block_groups | psql -h localhost -d detroit shp2pgsql census_block_groups_2010_utm17n.shp public.census2010_block_groups | psql -h localhost -d detroit
Software versions are:
- QGIS: 2.8.1-Wien
- GDAL/OGR: 1.11.2
- GEOS: 3.4.2-CAPI-1.8.2
- Proj.4: 480
- PostgreSQL: 9.3.6
- PostGIS: 2.1.0