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)?