1

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 cadastreenter image description here. 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)?

  • I don't think the sjoin call is doing much in your case, as geobuil and geocad (which you are joining) both only consist of a single geometry without any attributes. A small reproducible example would help. – joris Oct 5 '18 at 8:03
  • No, indeed, but it is all I found to achieve the merged attributs. I tried a pd.dataframe().merge() method, but i got an error saying the selected polygon with merged dataframe had a wrong geometry. – francois Oct 5 '18 at 12:13
  • I recently had to solve this problem; I ended up using shapely's intersection() and area(). Geopandas stores the geometries as shapely geometry types, and you can access the geometries easily via geoms = geocad.geometry.values. – Jon Oct 5 '18 at 13:51
2

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

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.