0

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

aogino
  • 1
  • Please provide sample data or at least what you chose from I assume https://apps.nationalmap.gov/downloader/ – Serge de Gosson de Varennes May 11 '23 at 13:23
  • I used the 25DN2_20 tile from the AHN4 pointcloud that you can find in geotiles.nl . The building footprints are extracted for the same extent of the point cloud from https://extract.bbbike.org/ . – aogino May 11 '23 at 14:05

0 Answers0