Programmiererforum

Code: Select all

# === STAGE: Imports and Environment ===
from whitebox_workflows import WbEnvironment
import whitebox_tools as wbt
import os, time, sys, shutil, glob

wbe = WbEnvironment()
wbe.verbose = True

# === STAGE: Define Directories and Parameters ===
base_dir = '/home/rdc/data/las_files/'
las_dir = f'{base_dir}LAS_orig/'
dtm_dir = f'{base_dir}ground/'
dsm_dir = f'{base_dir}first/'
big_dir = f'{base_dir}bigfiles/'
clean_dir = f'{base_dir}clean/'

max_tin_edge = 250  # Important parameter for TIN gridding

# === STAGE: Filter 'Ground' Points from LAS Files ===
wbe.working_directory = las_dir

all_files = os.listdir(las_dir)
las_files = [f for f in all_files if f.endswith('.las')]

for las in las_files:
lidar = wbe.read_lidar(las)
out_file = f'{dtm_dir}/{las[:-4]}_ground.las'
ground = wbe.filter_lidar_classes(lidar, exclusion_classes=[1, *range(3, 19)])
wbe.write_lidar(ground, out_file)

# === STAGE: Merge Filtered LAS Files ===
merged_las = 'lidar_ground_join.las'
las_files = glob.glob(rf'{dtm_dir}/*las')
lidars = wbe.read_lidars(*las_files)  # Removed 3 corrupted files manually
joined = wbe.lidar_join(lidars)
wbe.write_lidar(joined, merged_las)
lidars = None  # Clear memory

# === STAGE: Remove Overlaps ===
las_ground_final = f'{big_dir}lidar_ground_join.las'
las_no_overlap = os.path.join(clean_dir, 'lidar_no_overlap.las')

if not os.path.exists(las_no_overlap):
print(f"[Removing overlaps] {las_no_overlap}")
wbt.classify_overlap_points(
i=las_ground_final,
output=las_no_overlap,
resolution=1.0,
criterion="not max point source ID",
filter=True
)
else:
print(f"[SKIP] Overlap-free file already exists: {las_no_overlap}")

# === STAGE: Thinning LAS Data ===
las_thinned = os.path.join(clean_dir, 'lidar_thinned.las')

if not os.path.exists(las_thinned):
print(f"[Thinning] {las_thinned}")
wbt.lidar_thin(
i=las_no_overlap,
output=las_thinned,
resolution=2.0,
method="lowest",
save_filtered=False
)
else:
print(f"[SKIP] Already thinned file: {las_thinned}")

# === STAGE: Generate DEM via TIN Gridding ===
lidar = wbe.read_lidar(las_thinned)
max_tin_edge = 500  # Updated for TIN gridding

dem = wbe.lidar_tin_gridding(
input_lidar=lidar,
interpolation_parameter='elevation',
returns_included='all',  # Include single returns
excluded_classes=[1, *range(3, 10), 17],  # Keep class 2 ('ground')
cell_size=1.0,
max_triangle_edge_length=max_tin_edge
)

dem_out = f'{las_thinned[:-4]}_tin{max_tin_edge}.tif'
wbe.write_raster(dem, dem_out, compress=True)

# === STAGE: Fill DEM Voids using IDW ===
dem_fill = wbe.fill_missing_data(
dem=dem,
filter_size=51,
weight=2.0,
exclude_edge_nodata=True
)

dem_fill_out = f'{las_thinned[:-4]}_tin{max_tin_edge}_fill51.tif'
wbe.write_raster(dem_fill, dem_fill_out, compress=True)

# === STAGE: Generate Hillshade from Filled DEM ===
shade = f'{dem_fill_out[:-4]}_shade360.tif'
hillshade = wbe.multidirectional_hillshade(dem, full_360_mode=False)
wbe.write_raster(hillshade, shade, compress=True)