I'm working with two datasets: a set of precincts and a set of census block groups. The census block groups include water; the precincts don't. I would like to trim the census block groups to get only the land area of the census block groups (for a later calculation).

To do this, I find all the precincts that intersect the census block group, and calculate their intersections. I then union all of the precinct intersections together. The end result would, ideally, be the land area of the census block group. However, this is not the case. Rather, the intersections return a GeometryCollection of assorted LineStrings, MultiLineStrings, and Polygons. Attempts to union this collection are unsuccessful (the end goal would be to union it all into a polygon). Here's a rough example: example map

The translucent red/pink area would be the census block group, the light purple off to the right would be the water area of the census block group, and the blue, orange, and red areas would be precincts (land only). The aim is to get the areas of the translucent figure that intersect with the blue, orange, and red polygons, and then union those areas together to get the land area of the census block group.

I'm currently using GeoTools, a Java library built on the JTS implementation the OpenGIS standard, but I can use/translate a PostGIS query. I would also be open to any suggestions on better ways to do this (I'm sure there are a few).

The actual shapefiles can be obtained from the Harvard Dataverse Database (MD precinct map), and the census website (MD statewide census block group map).

  • Could you add the actual query you ran, and table/column info? Aren't the block groups corresponding to water bodies denoted as such or are that true off-shore regions only? – ThingumaBob Aug 13 '19 at 21:25
  • @ThingumaBob Block groups can cover both land and water. For example, if there is a large bay separating two block groups, the extents of the block groups will go to the middle of the bay. There are a few water-only precincts, but, as you said, they only encompass offshore areas. I've been using Java for this, which has resulted in a lot of verbose code, but I'll add a translated PostGIS query that has the same problem. – Alex Quilliam Aug 13 '19 at 21:50

I've solved my problem. I had to move away from GeoTools and write a PostGIS query, because GeoTools doesn't appear have the features needed to do this successfully.

Here's the query:

SELECT intersection_query.census_block_group_gid, ST_Multi(ST_Union(intersection_query.geom)) AS geom
 SELECT md_census_block_groups.gid AS census_block_group_gid, md_precincts.gid AS precinct_gid, ST_Intersection(md_census_block_groups.geom, md_precincts.geom) AS geom
 FROM md_census_block_groups, md_precincts
 WHERE ST_Intersects(md_census_block_groups.geom, md_precincts.geom)
) AS intersection_query
GROUP BY intersection_query.census_block_group_gid;

md_census_block_groups are the census block groups, md_precincts are the precincts.

  • Generally speaking unioning a set of geometries is fairly efficient. So It might be faster to take the union of all precincts which intersect a census block, and then compute the intersection of the union with the census block. The end result should be the same. – dr_jts Sep 7 '19 at 19:15

Not the answer you're looking for? Browse other questions tagged or ask your own question.