1

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.

Result of intersection; note the gaps!

Here are the 2000 (blue) and 2010 (red) block group boundaries--you can see there are other intersections that should have been found!

enter image description here

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").

enter image description here

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
  • Odd, it has to be said. Perhaps you could provide your raw data, as it is hard to see anything wrong with your query. – John Powell Feb 10 '16 at 10:24
  • Sure, I added a short shell script to download and subset the data using GDAL/OGR and shp2pgsql. – Arthur Feb 10 '16 at 16:07
  • Arthur, I'm a bit surprised that the output of this query was loadable in qgis, because ST_Intersects usually gives some GeometryCollections which are unreadable for qgis. I checked with your data and query and this was indeed the case. Also I am missing the fips columns in the data. Are you sure you this is the only query you used? In the end I had a complete coverage without gaps on my ST_Intersection, so it should be possible. – tilt Feb 11 '16 at 20:54
  • Hi @tilt thanks for your feedback. With regards to GeometryCollections in QGIS, there are three geom features associated with the layer when I open it in QGIS but two of them grayed out; I choose the one that is not and which is defined as a Polygon feature and it displays just fine (gaps notwithstanding). With regards to the FIPS code, in my original post I changed the field names in the SQL script for readability. I have now changed it back so it matches the attribute table. – Arthur Feb 11 '16 at 23:44
3

I should have seen this earlier: Since ST_Intersection gives geometrycollections in cases where the result gives more than one type of geometry (this happens especially when intersecting near-equal geometries), you have to extract the geometries from the collection first. I always put a (ST_Dump(ST_Intersection(geoma, geomb))).geom around the intersection, regardless of the expected output (note the extra brackets around dump) That should solve your issue.

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.