Add LiDAR point cloud support via OptiX sphere primitives

examples/nyc_lidar.py

@@ -0,0 +1,107 @@
+"""
+NYC Real LiDAR Point Cloud Demo
+================================
+
+Downloads real USGS 3DEP LiDAR LAZ tiles from the TNM API for
+Manhattan and visualizes them over the NYC DEM using OptiX sphere
+primitives with ASPRS classification coloring.
+
+First run downloads ~2.7 GB of LAZ data; subsequent runs use cache.
+LAZ decompression is parallelized across CPU cores.
+
+LiDAR point clouds and buildings are cached in the zarr store after
+the first run. Subsequent runs skip all LAZ decompression, filtering,
+thinning, and coordinate transforms — loading directly from zarr.
+
+Usage:
+    python examples/nyc_lidar.py
+"""
+
+from pathlib import Path
+
+import numpy as np
+import xarray as xr
+import zarr
+
+import rtxpy
+
+
+def _has_lidar_cache(zarr_path):
+    """Check if the zarr store has cached LiDAR sphere geometries."""
+    try:
+        store = zarr.open(str(zarr_path), mode='r', use_consolidated=False)
+        if 'meshes' not in store:
+            return False
+        mg = store['meshes']
+        for gid in mg:
+            if mg[gid].attrs.get('type', '') == 'sphere':
+                return True
+    except Exception:
+        pass
+    return False
+
+
+def main():
+    # Resolve paths relative to this script's directory
+    _here = Path(__file__).resolve().parent
+
+    # Full NYC DEM extent in WGS84
+    dem_bounds = (-74.26, 40.49, -73.70, 40.92)
+    # Manhattan LiDAR subset
+    lidar_bounds = (-74.02, 40.70, -73.97, 40.88)
+    cache_dir = _here / 'cache' / 'lidar'
+    zarr_path = str(_here / 'nyc_dem.zarr')
+
+    # Load NYC DEM (write CRS back — rioxarray doesn't auto-detect from zarr)
+    ds = xr.open_zarr(zarr_path)
+    terrain = ds['elevation'].load().astype(np.float32)
+    terrain = terrain.rio.write_crs('EPSG:32618')
+    terrain = terrain.rtx.to_cupy()
+
+    # Build terrain mesh
+    terrain.rtx.triangulate()
+
+    # Satellite tiles
+    terrain.rtx.place_tiles('satellite')
+
+    # Check for cached lidar + buildings in zarr
+    if _has_lidar_cache(zarr_path):
+        print("Loading cached meshes from zarr...")
+        terrain.rtx.load_meshes(zarr_path)
+    else:
+        # Use cached LAZ files if available, otherwise download
+        cached = sorted(cache_dir.glob('*.laz')) if cache_dir.exists() else []
+        if cached:
+            laz_paths = cached
+            print(f"Using {len(laz_paths)} cached LAZ file(s)")
+        else:
+            laz_paths = rtxpy.fetch_lidar(lidar_bounds, cache_dir=str(cache_dir))
+        print(f"Got {len(laz_paths)} LAZ file(s)")
+
+        # Buildings from Overture Maps (full DEM extent)
+        buildings = rtxpy.fetch_buildings(
+            bounds=dem_bounds, source='overture',
+            cache_path=str(_here / 'cache' / 'nyc_lidar_buildings.geojson'))
+        terrain.rtx.place_buildings(buildings)
+
+        # Place all LAZ tiles in parallel (threaded LAZ decompression)
+        # ASPRS classes: 1=unclassified, 2=ground, 10=rail, 17=bridge
+        # Exclude: 7=low noise, 9=water, 18=high noise
+        # thin=0.5 removes flight-line overlap striations (1 point per 0.5m cell)
+        terrain.rtx.place_pointclouds(
+            laz_paths,
+            point_size=0.25,
+            color='classification',
+            classification=[1, 2, 10, 17],
+            thin=0.5,
+        )
+
+        # Save all meshes (buildings + lidar) to zarr for next run
+        terrain.rtx.save_meshes(zarr_path)
+
+    # Launch interactive viewer
+    terrain.rtx.explore(width=1600, height=1200)
+
+
+if __name__ == '__main__':
+    main()

rtxpy/__init__.py

@@ -1,15 +1,15 @@
 from .rtx import (
     RTX,
     has_cupy,
     get_device_count,
     get_device_properties,
     list_devices,
     get_current_device,
     get_capabilities,
 )
 from .mesh import (
     triangulate_terrain,
     voxelate_terrain,
     write_stl,
     load_glb,
     load_mesh,
@@ -28,7 +28,7 @@
 
 # Optional convenience — network helpers with lazy dependency checks
 try:
-    from .remote_data import fetch_dem, fetch_osm, fetch_buildings, fetch_roads, fetch_water, fetch_wind, fetch_firms, fetch_places, fetch_infrastructure, fetch_land_use, fetch_restaurant_grades, fetch_gtfs
+    from .remote_data import fetch_dem, fetch_lidar, fetch_osm, fetch_buildings, fetch_roads, fetch_water, fetch_wind, fetch_firms, fetch_places, fetch_infrastructure, fetch_land_use, fetch_restaurant_grades, fetch_gtfs
 except ImportError:
     pass