I have two shapefiles:

  • catchments: polygons for 18 catchment areas (tend to be 'regional')
  • DA: polygons for ~2300 Dissemination Areas from Statistics Canada (tend to be 'neighborhood sized')

I want to determine how many Dissemination Areas as in each catchment. For example, catchment A spans 3 DAs, while catchment B spans 5 DAs.

import geopandas as gpd

DA = gpd.read_file("path to shapefile")
catchments = gpd.read_file("path to shapefile")

DA = DA.to_crs(epsg = 26917)
catchments = catchments.to_crs(epsg = 26917)

catchment_names = catchments["CatchName"].unique()

selection = catchments.loc[catchments["CatchName"]==catchment_names[0]]
result = selection.intersection(DA)

The results GeoSeries is as long as DA, but the results are a few "Polygon Empty" and the rest of the items are 'None'.

1 Answer 1


I would use overlay for this, e.g.

result_gdf = geopandas.overlay( DA , catchments , how='intersection')

The resulting GDF will have a row for the intersection of each DA and catchment shape. catchments that spread across multiple DAs will result in multiple output rows. The columns will contain all the columns from the 2 input GDFs.

  • This worked - thank you very much!
    – pubb
    Commented Apr 1, 2021 at 13:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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