#!/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 South (EPSG:3735) 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 South (EPSG:3735) in US Survey Feet to Web Mercator (EPSG:3857) in meters # The lidar data uses FIPS 3402 (Ohio South) per metadata transformer = Transformer.from_crs("EPSG:3735", "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()