[Python] Is there a way to intersect 2d polygons with a point cloud and filtering them out...

I'm trying to intersect some building footprints in a shapefile with a point cloud. My goal is to filter out the polygons(i.e. the footprints) that intersect with a certain number of points (a threshold that I'm willing to fix) and output them onto a new shapefile.

This is the code i'm currently using :

import geopandas as gpd
import laspy
from shapely.geometry import Polygon
import fiona

# read in the las file using laspy
las_file = laspy.read(
"results/C2C_using_sampled_model/25DN2_20_AHN4_appeared_building_sampled_model.las")

# create a GeoDataFrame from the x,y coordinates of the point cloud
points = gpd.GeoDataFrame(geometry=gpd.points_from_xy(
las_file.x, las_file.y), crs="EPSG:7415")

# load the shapefile containing building footprints into a geopandas dataframe
buildings = gpd.read_file("buildings_shapefile/buildings_clipped.shp")

# Project the points of the point cloud onto the footprints
joined = gpd.sjoin(buildings, points, predicate='intersects')

# Filter the geodataframe to only keep the polygons that intersect with a certain number of points
threshold = 100
filtered = joined.groupby('index_right').filter(lambda x: len(x) >= threshold)

# Save the filtered geodataframe as a new shapefile
filtered.to_file("filtered_footprints.shp")



The code outputs a result (after like 40 minutes or so of processing), but the shapefile is empty.

I don't really know where the problem came from. Also I'm not very familiar with geopandas functions.

Any help would be very much appreciated !!!

Continue reading...