In [2]:
# ── Imports + Load point cloud ────────────────────────────
import laspy
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
import rasterio
from rasterio.transform import from_origin
plt.style.use('dark_background')
LAS_PATH = '../data/lidar/utah_state_capitol.las'
las = laspy.read(LAS_PATH)
print(f"✅ Point cloud loaded")
print(f" Total points: {len(las.points):,}")
print(f" X range: {las.x.min():.1f} to {las.x.max():.1f}")
print(f" Y range: {las.y.min():.1f} to {las.y.max():.1f}")
print(f" Z range: {las.z.min():.1f} to {las.z.max():.1f} m")
print(f" Point format: {las.point_format}")
print(f" Available dimensions: {list(las.point_format.dimension_names)}")
if 'classification' in las.point_format.dimension_names:
classes, counts = np.unique(las.classification, return_counts=True)
print(f"\n Classification codes present:")
for c, n in zip(classes, counts):
print(f" Class {c}: {n:,} points")
✅ Point cloud loaded
Total points: 2,515,216
X range: 424907.2 to 425199.1
Y range: 4514279.2 to 4514696.8
Z range: 1365.9 to 1458.6 m
Point format: <PointFormat(3, 0 bytes of extra dims)>
Available dimensions: ['X', 'Y', 'Z', 'intensity', 'return_number', 'number_of_returns', 'scan_direction_flag', 'edge_of_flight_line', 'classification', 'synthetic', 'key_point', 'withheld', 'scan_angle_rank', 'user_data', 'point_source_id', 'gps_time', 'red', 'green', 'blue']
Classification codes present:
Class 1: 1,675,602 points
Class 2: 836,431 points
Class 7: 3,183 points
In [3]:
# ── Filter noise + top-down visualization ─────────────────
mask_valid = las.classification != 7
x, y, z = las.x[mask_valid], las.y[mask_valid], las.z[mask_valid]
classification = las.classification[mask_valid]
ground_mask = classification == 2
print(f"✅ Filtered: {len(x):,} valid points ({ground_mask.sum():,} ground, {(~ground_mask).sum():,} non-ground)")
fig, ax = plt.subplots(figsize=(10, 12))
scatter = ax.scatter(x, y, c=z, cmap='turbo', s=0.5, alpha=0.6)
ax.set_title('Sprint 06 · Utah State Capitol · Point Cloud (Top-Down, colored by elevation)', color='white')
ax.set_xlabel('Easting (m)', color='white')
ax.set_ylabel('Northing (m)', color='white')
ax.set_aspect('equal')
ax.tick_params(colors='white')
plt.colorbar(scatter, label='Elevation (m)')
plt.tight_layout()
plt.savefig('../data/lidar_pointcloud_topdown.png', dpi=150, bbox_inches='tight')
plt.show()
print("✅ Point cloud visualization generado")
✅ Filtered: 2,512,033 valid points (836,431 ground, 1,675,602 non-ground)
✅ Point cloud visualization generado
In [4]:
# ── Generate DTM (ground) and DSM (surface) rasters ───────
from scipy.stats import binned_statistic_2d
RESOLUTION = 1.0 # meters per pixel
x_min, x_max = x.min(), x.max()
y_min, y_max = y.min(), y.max()
n_cols = int(np.ceil((x_max - x_min) / RESOLUTION))
n_rows = int(np.ceil((y_max - y_min) / RESOLUTION))
x_edges = np.linspace(x_min, x_max, n_cols + 1)
y_edges = np.linspace(y_min, y_max, n_rows + 1)
# DTM: mean elevation of GROUND points per cell
dtm, _, _, _ = binned_statistic_2d(
x[ground_mask], y[ground_mask], z[ground_mask],
statistic='mean', bins=[x_edges, y_edges]
)
# DSM: max elevation of ALL points per cell (highest surface)
dsm, _, _, _ = binned_statistic_2d(
x, y, z,
statistic='max', bins=[x_edges, y_edges]
)
# Transpose so rows=Y (north-up) and flip vertically for correct orientation
dtm = np.flipud(dtm.T)
dsm = np.flipud(dsm.T)
print(f"✅ Grid generated: {n_rows} rows x {n_cols} cols @ {RESOLUTION}m/px")
print(f" DTM valid cells: {np.sum(~np.isnan(dtm)):,} / {dtm.size:,}")
print(f" DSM valid cells: {np.sum(~np.isnan(dsm)):,} / {dsm.size:,}")
fig, axes = plt.subplots(1, 2, figsize=(16, 8))
fig.suptitle('Sprint 06 · DTM vs DSM', fontsize=13, color='white')
im0 = axes[0].imshow(dtm, cmap='terrain', extent=[x_min, x_max, y_min, y_max])
axes[0].set_title('DTM (Ground / Terrain)', color='white')
axes[0].tick_params(colors='white')
plt.colorbar(im0, ax=axes[0], label='Elevation (m)')
im1 = axes[1].imshow(dsm, cmap='terrain', extent=[x_min, x_max, y_min, y_max])
axes[1].set_title('DSM (Surface incl. structures)', color='white')
axes[1].tick_params(colors='white')
plt.colorbar(im1, ax=axes[1], label='Elevation (m)')
plt.tight_layout()
plt.savefig('../data/lidar_dtm_dsm.png', dpi=150, bbox_inches='tight')
plt.show()
print("✅ DTM/DSM generado")
✅ Grid generated: 418 rows x 292 cols @ 1.0m/px DTM valid cells: 102,251 / 122,056 DSM valid cells: 121,992 / 122,056
✅ DTM/DSM generado
In [5]:
# ── Fill DTM gaps + compute height model ──────────────────
from scipy.interpolate import griddata
rows, cols = np.indices(dtm.shape)
valid_mask = ~np.isnan(dtm)
dtm_filled = griddata(
(rows[valid_mask], cols[valid_mask]), dtm[valid_mask],
(rows, cols), method='nearest'
)
height_model = dsm - dtm_filled
height_model = np.where(np.isnan(dsm), np.nan, height_model)
height_model = np.clip(height_model, 0, None) # no negative heights
print(f"✅ Height model generado")
print(f" Max height above ground: {np.nanmax(height_model):.1f} m")
print(f" Mean height above ground: {np.nanmean(height_model):.2f} m")
fig, ax = plt.subplots(figsize=(10, 12))
im = ax.imshow(height_model, cmap='inferno', extent=[x_min, x_max, y_min, y_max])
ax.set_title('Sprint 06 · Height Above Ground (DSM - DTM)', color='white')
ax.set_xlabel('Easting (m)', color='white')
ax.set_ylabel('Northing (m)', color='white')
ax.tick_params(colors='white')
plt.colorbar(im, label='Height (m)')
plt.tight_layout()
plt.savefig('../data/lidar_height_model.png', dpi=150, bbox_inches='tight')
plt.show()
✅ Height model generado Max height above ground: 75.9 m Mean height above ground: 3.51 m
In [6]:
# ── Structure detection via height threshold ──────────────
THRESHOLD_LOW = 2.0 # meters - cualquier estructura sobre el suelo (vegetación + edificio)
THRESHOLD_HIGH = 10.0 # meters - solo estructuras de escala edificio
structure_mask_low = height_model > THRESHOLD_LOW
structure_mask_high = height_model > THRESHOLD_HIGH
area_low = np.sum(structure_mask_low) * (RESOLUTION ** 2)
area_high = np.sum(structure_mask_high) * (RESOLUTION ** 2)
print(f"✅ Structure detection")
print(f" Footprint > {THRESHOLD_LOW}m: {area_low:,.0f} m² ({area_low/10000:.2f} ha)")
print(f" Footprint > {THRESHOLD_HIGH}m (building-scale): {area_high:,.0f} m² ({area_high/10000:.2f} ha)")
fig, axes = plt.subplots(1, 2, figsize=(16, 8))
fig.suptitle('Sprint 06 · Structure Detection by Height Threshold', fontsize=13, color='white')
axes[0].imshow(dsm, cmap='gray', extent=[x_min, x_max, y_min, y_max])
overlay_low = np.ma.masked_where(~structure_mask_low, structure_mask_low)
axes[0].imshow(overlay_low, cmap='Reds', alpha=0.6, extent=[x_min, x_max, y_min, y_max])
axes[0].set_title(f'Threshold > {THRESHOLD_LOW}m (veg + structures)', color='white')
axes[0].tick_params(colors='white')
axes[1].imshow(dsm, cmap='gray', extent=[x_min, x_max, y_min, y_max])
overlay_high = np.ma.masked_where(~structure_mask_high, structure_mask_high)
axes[1].imshow(overlay_high, cmap='Reds', alpha=0.6, extent=[x_min, x_max, y_min, y_max])
axes[1].set_title(f'Threshold > {THRESHOLD_HIGH}m (building-scale)', color='white')
axes[1].tick_params(colors='white')
plt.tight_layout()
plt.savefig('../data/lidar_structure_detection.png', dpi=150, bbox_inches='tight')
plt.show()
print("✅ Structure detection map generado")
✅ Structure detection Footprint > 2.0m: 30,838 m² (3.08 ha) Footprint > 10.0m (building-scale): 16,387 m² (1.64 ha)
✅ Structure detection map generado
In [7]:
# ── Export to GeoTIFF + final summary ─────────────────────
from rasterio.transform import from_origin
from rasterio.crs import CRS
try:
crs = las.header.parse_crs()
epsg_code = crs.to_epsg()
print(f"CRS found in file: EPSG:{epsg_code}")
except Exception:
epsg_code = 26912
print(f"No CRS embedded — using EPSG:{epsg_code} (NAD83 UTM Zone 12N, Salt Lake City)")
transform = from_origin(x_min, y_max, RESOLUTION, RESOLUTION)
def save_raster(array, filename, dtype='float32'):
with rasterio.open(
filename, 'w', driver='GTiff',
height=array.shape[0], width=array.shape[1],
count=1, dtype=dtype,
crs=CRS.from_epsg(epsg_code), transform=transform,
nodata=np.nan
) as dst:
dst.write(array.astype(dtype), 1)
save_raster(dtm_filled, '../data/sprint06_dtm.tif')
save_raster(dsm, '../data/sprint06_dsm.tif')
save_raster(height_model, '../data/sprint06_height_model.tif')
print(f"\n✅ GeoTIFFs exportados: DTM, DSM, height model")
summary = f"""
╔══════════════════════════════════════════╗
║ SPRINT 06 · LIDAR STRUCTURE DETECTION ║
╠══════════════════════════════════════════╣
║ Site : Utah State Capitol
║ Total points : {len(las.points):,}
║ Ground points : {ground_mask.sum():,}
║ Grid resolution : {RESOLUTION}m/px
║ Max height : {np.nanmax(height_model):.1f}m
║ Footprint >2m : {area_low/10000:.2f} ha
║ Footprint >10m : {area_high/10000:.2f} ha
║ CRS : EPSG:{epsg_code}
║ Source : USGS / OpenTopography
╚══════════════════════════════════════════╝
"""
print(summary)
CRS found in file: EPSG:26912 ✅ GeoTIFFs exportados: DTM, DSM, height model ╔══════════════════════════════════════════╗ ║ SPRINT 06 · LIDAR STRUCTURE DETECTION ║ ╠══════════════════════════════════════════╣ ║ Site : Utah State Capitol ║ Total points : 2,515,216 ║ Ground points : 836,431 ║ Grid resolution : 1.0m/px ║ Max height : 75.9m ║ Footprint >2m : 3.08 ha ║ Footprint >10m : 1.64 ha ║ CRS : EPSG:26912 ║ Source : USGS / OpenTopography ╚══════════════════════════════════════════╝
In [ ]: