201 lines
6.6 KiB
Python
201 lines
6.6 KiB
Python
#!/usr/bin/env python3
|
|
"""
|
|
Convert LAS lidar files to .mound binary format for Three.js rendering.
|
|
|
|
Usage:
|
|
python las2mound.py input.las output.mound
|
|
|
|
.mound Binary Format Specification
|
|
==================================
|
|
|
|
Header (64 bytes):
|
|
Offset Size Type Description
|
|
------ ---- ---------- -----------
|
|
0 4 char[4] Magic number: 'LIDR'
|
|
4 4 uint32 Version number (currently 1)
|
|
8 4 uint32 Point count (number of vertices)
|
|
12 4 uint32 Triangle count (number of triangles)
|
|
16 4 float32 Min X coordinate (Web Mercator meters)
|
|
20 4 float32 Min Y coordinate (Web Mercator meters)
|
|
24 4 float32 Min Z coordinate (elevation in meters)
|
|
28 4 float32 Max X coordinate (Web Mercator meters)
|
|
32 4 float32 Max Y coordinate (Web Mercator meters)
|
|
36 4 float32 Max Z coordinate (elevation in meters)
|
|
40 24 bytes Reserved (padding to 64 bytes)
|
|
|
|
Vertex Data (point_count * 12 bytes):
|
|
Series of vertices in XYZ float32 triplets.
|
|
X, Y are in Web Mercator meters (EPSG:3857)
|
|
Z is elevation in meters
|
|
Total size: point_count * 3 * 4 bytes
|
|
|
|
Index Data (triangle_count * 12 bytes):
|
|
Series of triangle indices (3 uint32 per triangle).
|
|
Each index references a vertex in the vertex data.
|
|
Total size: triangle_count * 3 * 4 bytes
|
|
|
|
Example layout for 100 points and 50 triangles:
|
|
Bytes 0-63: Header
|
|
Bytes 64-1263: Vertex data (100 * 12)
|
|
Bytes 1264-1863: Index data (50 * 12)
|
|
Total file size: 1864 bytes
|
|
|
|
Coordinate System:
|
|
Input: Ohio State Plane North (EPSG:3734) in US Survey Feet
|
|
Output: Web Mercator (EPSG:3857) in meters
|
|
This ensures compatibility with MapLibre GL JS and web mapping standards.
|
|
"""
|
|
|
|
import sys
|
|
import struct
|
|
import numpy as np
|
|
from pathlib import Path
|
|
|
|
try:
|
|
import laspy
|
|
except ImportError:
|
|
print("Error: laspy not installed. Run: pip install laspy")
|
|
sys.exit(1)
|
|
|
|
try:
|
|
from scipy.spatial import Delaunay
|
|
except ImportError:
|
|
print("Error: scipy not installed. Run: pip install scipy")
|
|
sys.exit(1)
|
|
|
|
try:
|
|
from pyproj import Transformer
|
|
except ImportError:
|
|
print("Error: pyproj not installed. Run: pip install pyproj")
|
|
sys.exit(1)
|
|
|
|
|
|
def transform_to_webmercator(x, y, z):
|
|
"""Transform Ohio State Plane coordinates to Web Mercator (EPSG:3857)."""
|
|
print("Transforming coordinates to Web Mercator...")
|
|
|
|
# Ohio State Plane North (EPSG:3734) in US Survey Feet to Web Mercator (EPSG:3857) in meters
|
|
# Newark is in the North zone
|
|
transformer = Transformer.from_crs("EPSG:3734", "EPSG:3857", always_xy=True)
|
|
|
|
# Transform x,y (easting, northing) to Web Mercator meters
|
|
merc_x, merc_y = transformer.transform(x, y)
|
|
|
|
# Convert elevation from US Survey Feet to meters
|
|
z_meters = z * 0.3048006096012192
|
|
|
|
print(f"Transformed to Web Mercator bounds:")
|
|
print(f" X (meters): [{merc_x.min():.2f}, {merc_x.max():.2f}] (span: {merc_x.max() - merc_x.min():.2f}m)")
|
|
print(f" Y (meters): [{merc_y.min():.2f}, {merc_y.max():.2f}] (span: {merc_y.max() - merc_y.min():.2f}m)")
|
|
print(f" Z (meters): [{z_meters.min():.2f}, {z_meters.max():.2f}] (span: {z_meters.max() - z_meters.min():.2f}m)")
|
|
|
|
return merc_x, merc_y, z_meters
|
|
|
|
|
|
def read_las(filepath):
|
|
"""Read LAS file and extract ground points."""
|
|
print(f"Reading {filepath}...")
|
|
las = laspy.read(filepath)
|
|
|
|
# Extract coordinates
|
|
x = las.x
|
|
y = las.y
|
|
z = las.z
|
|
|
|
# Filter for ground points (classification 2) if available
|
|
if hasattr(las, 'classification'):
|
|
ground_mask = las.classification == 2
|
|
if ground_mask.any():
|
|
print(f"Filtering {ground_mask.sum()} ground points from {len(x)} total points")
|
|
x = x[ground_mask]
|
|
y = y[ground_mask]
|
|
z = z[ground_mask]
|
|
|
|
print(f"Loaded {len(x)} points")
|
|
return x, y, z
|
|
|
|
|
|
def triangulate_points(x, y, z):
|
|
"""Perform Delaunay triangulation on XY coordinates."""
|
|
print("Performing Delaunay triangulation...")
|
|
|
|
# Create 2D point array for triangulation
|
|
points_2d = np.column_stack([x, y])
|
|
|
|
# Triangulate
|
|
tri = Delaunay(points_2d)
|
|
|
|
print(f"Generated {len(tri.simplices)} triangles")
|
|
|
|
return tri.simplices
|
|
|
|
|
|
def write_mound(filepath, x, y, z, indices):
|
|
"""Write .mound binary format."""
|
|
print(f"Writing {filepath}...")
|
|
|
|
# Convert to float32
|
|
positions = np.column_stack([x, y, z]).astype(np.float32)
|
|
indices = indices.astype(np.uint32)
|
|
|
|
# Calculate bounds
|
|
min_x, min_y, min_z = positions.min(axis=0)
|
|
max_x, max_y, max_z = positions.max(axis=0)
|
|
|
|
point_count = len(positions)
|
|
triangle_count = len(indices)
|
|
|
|
with open(filepath, 'wb') as f:
|
|
# Header (64 bytes)
|
|
f.write(b'LIDR') # Magic number (4 bytes)
|
|
f.write(struct.pack('I', 1)) # Version (4 bytes)
|
|
f.write(struct.pack('I', point_count)) # Point count (4 bytes)
|
|
f.write(struct.pack('I', triangle_count)) # Triangle count (4 bytes)
|
|
f.write(struct.pack('f', min_x)) # Min X (4 bytes)
|
|
f.write(struct.pack('f', min_y)) # Min Y (4 bytes)
|
|
f.write(struct.pack('f', min_z)) # Min Z (4 bytes)
|
|
f.write(struct.pack('f', max_x)) # Max X (4 bytes)
|
|
f.write(struct.pack('f', max_y)) # Max Y (4 bytes)
|
|
f.write(struct.pack('f', max_z)) # Max Z (4 bytes)
|
|
f.write(b'\x00' * 24) # Reserved (24 bytes)
|
|
|
|
# Position data
|
|
f.write(positions.tobytes())
|
|
|
|
# Index data
|
|
f.write(indices.tobytes())
|
|
|
|
file_size = Path(filepath).stat().st_size / (1024 * 1024)
|
|
print(f"Wrote {point_count} points, {triangle_count} triangles ({file_size:.2f} MB)")
|
|
print(f"Bounds: X[{min_x:.2f}, {max_x:.2f}] Y[{min_y:.2f}, {max_y:.2f}] Z[{min_z:.2f}, {max_z:.2f}]")
|
|
|
|
|
|
def main():
|
|
if len(sys.argv) != 3:
|
|
print("Usage: python las2mound.py input.las output.mound")
|
|
sys.exit(1)
|
|
|
|
input_file = sys.argv[1]
|
|
output_file = sys.argv[2]
|
|
|
|
if not Path(input_file).exists():
|
|
print(f"Error: Input file '{input_file}' not found")
|
|
sys.exit(1)
|
|
|
|
# Read LAS
|
|
x, y, z = read_las(input_file)
|
|
|
|
# Transform to Web Mercator
|
|
merc_x, merc_y, z_meters = transform_to_webmercator(x, y, z)
|
|
|
|
# Triangulate (using Web Mercator coordinates)
|
|
indices = triangulate_points(merc_x, merc_y, z_meters)
|
|
|
|
# Write output (Web Mercator X/Y in meters, elevation in meters)
|
|
write_mound(output_file, merc_x, merc_y, z_meters, indices)
|
|
|
|
print("Done!")
|
|
|
|
|
|
if __name__ == '__main__':
|
|
main() |