Joining attributes after selecting one polygon which intersects another using geopandas?

Question

To select the polygons which overlap other layer's polygons, I came up with this code based on what I found on this site:

import geopandas as gpd
import fiona
import os
import sys

# Get the current working directory
path = os.getcwd() 

# Parameters to fix
input1 = 'building.shp'    #polygons to keep from
input2 = 'cadastre.shp'    #polygons which must be intersected
outputname = 'intersect_selected.shp'
intersect = 0.75


# Read data
building = gpd.read_file(os.path.join(path,input1))
cadastre = gpd.read_file(os.path.join(path,input2))

# Check crs
if building.crs != cadastre.crs:
    print('CRS not match \nEnding script.')
    sys.exit()

#List crs
crs3 = building.crs

#Processing iteration
data = gpd.GeoDataFrame()
for index1, cad in cadastre.iterrows():
    for index2, buil in building.iterrows(): 
        if buil['geometry'].intersects(cad['geometry']):
            area_int = buil['geometry'].intersection(cad['geometry']).area
            area_buil = buil['geometry'].area
            area_cad = cad['geometry'].area
            crit1 = area_int/area_buil
            crit2 = area_int/area_cad
            if crit1 > intersect or crit2 == 1.0:
                geobuil = gpd.GeoDataFrame([buil], crs=crs3)
                geocad = gpd.GeoDataFrame([cad], crs=crs3)
                intersec = gpd.sjoin(geobuil, geocad, op='intersects')
                data = data.append(intersec)

#Export data
data.to_file('intersection.shp')

Example: I am trying to select the polygons from building.shp which intersect more than 75% the ones from cadastre.shp or that contain polygons from cadastre. Once selected, I would like to add all the attributs from cadastre features.

Is the .sjoin() method not too much time consuming? Is there a better way to join the attributs of concerned polygons (cad and buil)?

Answer 1

Use GeoPandas Overlay

# areas of buildings
building['area_buil'] =building.geometry.apply(lambda x: x.area)
# first row
print(building.head(1))
   id               geometry                                area_buil
0   1  POLYGON ((682.0664490861618 -56.9528720626632,...  14234.471441
# areas of cadastre parcels
cadastre['area_cad'] =cadastre.geometry.apply(lambda x: x.area)
# first row
cadastre.head()
    id            geometry                                 area_cad
0   1  POLYGON ((522.2243230638396 -72.21015374517569...  859.121096

Intersection of the two GeoDataFrames

inter = gpd.overlay(building, cadastre, how='intersection')
# areas of intersection polygons
inter['area_int'] =inter.geometry.apply(lambda x: x.area)
# first row
inter.head(1)
    id   area_buil    id_2   area_cad   geometry    area_int
0   1  14234.471441     3  2569.372692    POLYGON ((686.3676624325899 -106.6574169572678   2569.372692 

Now you can compute crit1 and crit2

inter['crit1'] =inter.apply(lambda row: row.area_int/row.area_buil,axis=1)
inter['crit2'] =inter.apply(lambda row: row.area_int/row.area_cad,axis=1)

The result is a GeoDataFrame with everything you need.

  inter[['area_int','area_buil','crit1','area_cad','crit2']]
       area_int     area_buil     crit1     area_cad     crit2
  0  2569.372692  14234.471441  0.180504  2569.372692  1.000000
  1   776.791694    913.700570  0.850160  1791.819989  0.433521
  2   247.790095    545.270793  0.454435   859.121096  0.288423
  .....