Visualize Johson Highway surface

roadsurface.py

import laspy
from sklearn.cluster import OPTICS
import numpy as np
from scipy.spatial import KDTree
import pyvista as pv
import geopandas as gpd
import pandas as pd
from shapely.geometry import box
from shapely.geometry import Point
import matplotlib.pyplot as plt

# Load the LAS file
las_file = "/home/dylanblee/USGS_LPC_TN_Eastern_TN_LiDAR_2016_B16_2640629NE.las"
las_data = laspy.read(las_file)
points = np.vstack((las_data.x, las_data.y, las_data.z)).transpose()

# Load the road area
geo_df = gpd.read_file('/home/dylanblee/JohnsonHighway.geojson')

# clip road area to LAS bounding box
bbox = box(-83.7028063143902, 36.05100475507778, -83.6788119843394, 36.06244358817336)
cropped_geometries = geo_df.intersection(bbox)
geo_df['geometry'] = cropped_geometries

# clip LAS points to just the road
gdf_points = gpd.GeoDataFrame(
    pd.DataFrame(points, columns=['x', 'y', 'z']),
    geometry=gpd.points_from_xy(points[:, 0], points[:, 1], points[:, 2]),
    crs="EPSG:6576"
)

#reproject points
gdf_points = gdf_points.to_crs("EPSG:3857")

gdf_geometries = gpd.GeoDataFrame(geometry=cropped_geometries,crs=4326)
gdf_geometries = gdf_geometries.to_crs("EPSG:3857")

buffer_distance = 20 #m since using a projected CRS
gdf_geometries['geometry'] = gdf_geometries.buffer(buffer_distance)

joined = gpd.sjoin(gdf_points, gdf_geometries, how="inner", op="within")

clipped_points = joined[['x', 'y', 'z']].values

#visualize
cloud = pv.PolyData(clipped_points)

# Generate the mesh using Delaunay triangulation
mesh = cloud.delaunay_2d()

plotter = pv.Plotter()
plotter.add_mesh(mesh, show_edges=True, color="lightgrey")
plotter.enable_terrain_style(mouse_wheel_zooms=True, shift_pans=True)
plotter.show()