LiDAR and point clouds

"""A synthetic airborne LiDAR survey of 120 m by 120 m on the eastern flank of the module's hill, in metres of UTM zone
43N with heights above the geoid: about two pulses a square metre, each giving one to three returns.

Every point is (x, y, z, intensity, return number, number of returns, true class), with the true class an ASPRS code:
2 ground, 5 high vegetation, 6 building, 7 low point (noise), 18 high noise. The truth lets a classifier be measured.
"""
import math
import random

WEST, SOUTH, SIZE = 374_700.0, 2_042_940.0, 120.0
BUILDINGS = [(374_712.0, 2_043_010.0, 30.0, 24.0, 9.0), (374_760.0, 2_042_950.0, 48.0, 40.0, 7.0)]  # x, y, width, depth, height
rng0 = random.Random(5)
TREES = [(WEST + rng0.uniform(5, 115), SOUTH + rng0.uniform(70, 115), rng0.uniform(3, 5), rng0.uniform(8, 15)) for _ in range(18)]


def ground(x, y):
    """The hill's bare ground (m above the geoid), as in the elevation models lesson."""
    u = 1 - ((x - 374_000.0) / 1_000.0) ** 2 - ((y - 2_043_000.0) / 800.0) ** 2
    return 520.0 + 150.0 * max(0.0, u) ** 2


def roof(x, y):
    """The height of a flat roof over (x, y), or None: each roof is level, its height above the ground's highest corner."""
    for bx, by, w, d, h in BUILDINGS:
        if bx <= x < bx + w and by <= y < by + d:
            return max(ground(bx, by), ground(bx + w, by), ground(bx, by + d), ground(bx + w, by + d)) + h
    return None


def canopy(x, y):
    """The height of a tree crown over (x, y), or None (crowns are domes)."""
    best = None
    for tx, ty, r, h in TREES:
        d = math.hypot(x - tx, y - ty)
        if d < r:
            top = ground(tx, ty) + h - (h * 0.4) * (d / r) ** 2
            best = top if best is None else max(best, top)
    return best


def survey(seed=3, spacing=0.7):
    """All points of the survey: pulses on a jittered grid, returns from roofs, crowns and the ground, and noise."""
    rng = random.Random(seed)
    points = []
    steps = int(SIZE / spacing)
    for i in range(steps):
        for j in range(steps):
            x = WEST + (i + rng.random()) * spacing
            y = SOUTH + (j + rng.random()) * spacing
            g = ground(x, y) + rng.gauss(0, 0.05)
            if (r := roof(x, y)) is not None:
                points.append((x, y, r + rng.gauss(0, 0.05), 820 + rng.randint(-40, 40), 1, 1, 6))
            elif (c := canopy(x, y)) is not None:
                n = 3 if rng.random() < 0.3 else 2 if rng.random() < 0.7 else 1  # the pulse may not reach the ground
                points.append((x, y, c + rng.gauss(0, 0.1), 160 + rng.randint(-30, 30), 1, n, 5))
                if n == 3:
                    points.append((x, y, (c + g) / 2 + rng.gauss(0, 0.5), 90 + rng.randint(-20, 20), 2, 3, 5))
                if n > 1:
                    points.append((x, y, g, 120 + rng.randint(-30, 30), n, n, 2))
            else:
                points.append((x, y, g, 300 + rng.randint(-50, 50), 1, 1, 2))
    for _ in range(40):  # low noise: echoes that come back from below the ground
        x, y = WEST + rng.uniform(0, SIZE), SOUTH + rng.uniform(0, SIZE)
        points.append((x, y, ground(x, y) - rng.uniform(3, 10), rng.randint(5, 40), 1, 1, 7))
    for _ in range(25):  # high noise: birds and haze far above the surface
        x, y = WEST + rng.uniform(0, SIZE), SOUTH + rng.uniform(0, SIZE)
        points.append((x, y, ground(x, y) + rng.uniform(40, 80), rng.randint(5, 40), 1, 1, 18))
    return points

"""One pulse through a tree crown, its three returns as LAS point records of format 6 (30 bytes each), and back."""
import struct

from cloud import survey

SCALE, OFFSET = (0.01, 0.01, 0.01), (374_700.0, 2_042_900.0, 500.0)  # from the file's header
FORMAT6 = "<iiiHBBBBhHd"  # X Y Z intensity returns flags class user-data scan-angle source gps-time


def encode(p, gps_time):
    x, y, z, intensity, ret, n, cls = p
    ix, iy, iz = (round((v - o) / s) for v, o, s in zip((x, y, z), OFFSET, SCALE))  # coordinates as scaled integers
    return struct.pack(FORMAT6, ix, iy, iz, intensity, ret | n << 4, 0, cls, 0, 0, 1, gps_time)


def decode(raw):
    ix, iy, iz, intensity, returns, _, cls, *_rest, gps_time = struct.unpack(FORMAT6, raw)
    x, y, z = (i * s + o for i, s, o in zip((ix, iy, iz), SCALE, OFFSET))
    return (ix, iy, iz), (x, y, z), intensity, returns & 0x0F, returns >> 4, cls


points = survey()
k = next(i for i, p in enumerate(points) if p[5] == 3)  # the first pulse that came back three times
pulse = points[k : k + 3]
records = [encode(p, 1_000.0 + 0.000_01 * k) for p in pulse]
print(f"record size: {struct.calcsize(FORMAT6)} bytes; the first return as stored:\n  {records[0].hex(' ')}")
for raw in records:
    ints, (x, y, z), intensity, ret, n, cls = decode(raw)
    print(f"return {ret} of {n}: stored {ints}, so x {x:,.2f}, y {y:,.2f}, z {z:.2f} m; intensity {intensity}, class {cls}")
print(f"crown to ground along this pulse: {pulse[0][2] - pulse[2][2]:.2f} m")

record size: 30 bytes; the first return as stored:
  fe 00 00 00 41 36 00 00 24 19 00 00 ae 00 31 00 05 00 00 00 01 00 b2 4b 54 6f 0d 40 8f 40
return 1 of 3: stored (254, 13889, 6436), so x 374,702.54, y 2,043,038.89, z 564.36 m; intensity 174, class 5
return 2 of 3: stored (254, 13889, 6094), so x 374,702.54, y 2,043,038.89, z 560.94 m; intensity 99, class 5
return 3 of 3: stored (254, 13889, 5809), so x 374,702.54, y 2,043,038.89, z 558.09 m; intensity 96, class 2
crown to ground along this pulse: 6.27 m

"""The survey's points by class and by return: what a LAS file of them holds, and what it would weigh."""
from cloud import BUILDINGS, TREES, survey

CODES = {2: "ground", 5: "high vegetation", 6: "building", 7: "low point (noise)", 18: "high noise"}
points = survey()
print(f"{len(points):,} points over 120 m by 120 m: {len(points) / 14_400:.2f} a square metre")
for code, name in CODES.items():
    n = sum(p[6] == code for p in points)
    print(f"  class {code:>2} {name:<18} {n:>6,}")
first = sum(p[4] == 1 for p in points)
last_ground = sum(p[6] == 2 and p[5] > 1 for p in points)
print(f"first returns: {first:,}; ground points that were the last of several returns: {last_ground:,}")
print(f"{len(BUILDINGS)} buildings and {len(TREES)} trees")
print(f"as LAS 1.5, point format 6: a 393-byte header and {len(points):,} records of 30 bytes, {393 + 30 * len(points):,} bytes before the VLRs")

30,952 points over 120 m by 120 m: 2.15 a square metre
  class  2 ground             23,595
  class  5 high vegetation     1,922
  class  6 building            5,370
  class  7 low point (noise)      40
  class 18 high noise             25
first returns: 29,306; ground points that were the last of several returns: 1,200
2 buildings and 18 trees
as LAS 1.5, point format 6: a 393-byte header and 30,952 records of 30 bytes, 928,953 bytes before the VLRs