drone-site-survey

Drone Site Survey Processing

Overview

This skill implements drone data processing for construction site monitoring. Process aerial imagery to generate maps, measure volumes, track progress, and compare with design models.

Capabilities:

• Orthomosaic generation
• Digital Elevation Model (DEM) creation
• Point cloud processing
• Volume calculations
• Progress monitoring
• BIM comparison
• Stockpile measurement

Quick Start

from dataclasses import dataclass
from typing import List, Dict, Tuple, Optional
from datetime import datetime
import numpy as np

@dataclass
class DroneImage:
    filename: str
    timestamp: datetime
    latitude: float
    longitude: float
    altitude: float
    heading: float
    pitch: float
    roll: float
    camera_model: str

@dataclass
class PointCloud:
    points: np.ndarray  # Nx3 array
    colors: Optional[np.ndarray] = None  # Nx3 RGB
    normals: Optional[np.ndarray] = None  # Nx3

@dataclass
class VolumeResult:
    volume_m3: float
    area_m2: float
    method: str
    reference_plane: str
    confidence: float

def calculate_volume_simple(point_cloud: PointCloud,
                           reference_z: float = None) -> VolumeResult:
    """Simple volume calculation from point cloud"""
    points = point_cloud.points

    if reference_z is None:
        reference_z = np.min(points[:, 2])

    # Grid-based volume calculation
    x_min, x_max = np.min(points[:, 0]), np.max(points[:, 0])
    y_min, y_max = np.min(points[:, 1]), np.max(points[:, 1])

    grid_size = 0.5  # 50cm grid
    x_bins = np.arange(x_min, x_max + grid_size, grid_size)
    y_bins = np.arange(y_min, y_max + grid_size, grid_size)

    volume = 0
    cell_area = grid_size ** 2

    for i in range(len(x_bins) - 1):
        for j in range(len(y_bins) - 1):
            mask = (
                (points[:, 0] >= x_bins[i]) & (points[:, 0] < x_bins[i + 1]) &
                (points[:, 1] >= y_bins[j]) & (points[:, 1] < y_bins[j + 1])
            )
            cell_points = points[mask]
            if len(cell_points) > 0:
                max_z = np.max(cell_points[:, 2])
                height = max_z - reference_z
                if height > 0:
                    volume += height * cell_area

    area = (x_max - x_min) * (y_max - y_min)

    return VolumeResult(
        volume_m3=volume,
        area_m2=area,
        method='grid_based',
        reference_plane=f'z={reference_z:.2f}',
        confidence=0.9
    )

# Example usage
sample_points = np.random.rand(10000, 3) * [100, 100, 10]  # 100x100m, 10m height
point_cloud = PointCloud(points=sample_points)
result = calculate_volume_simple(point_cloud)
print(f"Volume: {result.volume_m3:.2f} m³, Area: {result.area_m2:.2f} m²")

Comprehensive Drone Survey System

Image Processing Pipeline

from dataclasses import dataclass, field
from typing import List, Dict, Tuple, Optional
from datetime import datetime
import numpy as np
from pathlib import Path
import json

@dataclass
class CameraParameters:
    focal_length_mm: float
    sensor_width_mm: float
    sensor_height_mm: float
    image_width_px: int
    image_height_px: int

@dataclass
class GeoReference:
    crs: str  # Coordinate Reference System (e.g., "EPSG:4326")
    origin: Tuple[float, float, float]  # lat, lon, alt
    rotation: Tuple[float, float, float]  # heading, pitch, roll

@dataclass
class SurveyFlight:
    flight_id: str
    date: datetime
    site_name: str
    images: List[DroneImage]
    camera: CameraParameters
    geo_reference: GeoReference
    flight_altitude: float
    overlap_forward: float = 0.8
    overlap_side: float = 0.7
    gsd: float = 0  # Ground Sample Distance (cm/pixel)

    def __post_init__(self):
        if self.gsd == 0 and self.camera:
            # Calculate GSD
            sensor_width = self.camera.sensor_width_mm
            focal_length = self.camera.focal_length_mm
            image_width = self.camera.image_width_px
            altitude = self.flight_altitude

            self.gsd = (altitude * sensor_width) / (focal_length * image_width) * 100  # cm

@dataclass
class ProcessingResult:
    orthomosaic_path: Optional[str] = None
    dem_path: Optional[str] = None
    dsm_path: Optional[str] = None
    point_cloud_path: Optional[str] = None
    report_path: Optional[str] = None
    statistics: Dict = field(default_factory=dict)

class DroneDataProcessor:
    """Process drone survey data"""

    def __init__(self, output_dir: str):
        self.output_dir = Path(output_dir)
        self.output_dir.mkdir(parents=True, exist_ok=True)

    def process_survey(self, flight: SurveyFlight,
                      generate_ortho: bool = True,
                      generate_dem: bool = True,
                      generate_pointcloud: bool = True) -> ProcessingResult:
        """Process drone survey data"""
        result = ProcessingResult()
        result.statistics['flight_id'] = flight.flight_id
        result.statistics['image_count'] = len(flight.images)
        result.statistics['gsd_cm'] = flight.gsd
        result.statistics['flight_date'] = flight.date.isoformat()

        # In production, these would call actual photogrammetry libraries
        # like OpenDroneMap, Pix4D API, or custom SfM pipeline

        if generate_ortho:
            result.orthomosaic_path = str(self.output_dir / f"{flight.flight_id}_ortho.tif")
            result.statistics['ortho_resolution'] = flight.gsd

        if generate_dem:
            result.dem_path = str(self.output_dir / f"{flight.flight_id}_dem.tif")
            result.dsm_path = str(self.output_dir / f"{flight.flight_id}_dsm.tif")

        if generate_pointcloud:
            result.point_cloud_path = str(self.output_dir / f"{flight.flight_id}_pointcloud.las")

        # Generate report
        result.report_path = str(self.output_dir / f"{flight.flight_id}_report.json")
        with open(result.report_path, 'w') as f:
            json.dump(result.statistics, f, indent=2)

        return result

    def extract_point_cloud(self, las_path: str) -> PointCloud:
        """Extract point cloud from LAS file"""
        # In production, use laspy or similar
        # Simulated point cloud for demonstration
        n_points = 100000
        points = np.random.rand(n_points, 3) * [100, 100, 20]
        colors = np.random.randint(0, 255, (n_points, 3), dtype=np.uint8)

        return PointCloud(points=points, colors=colors)

    def compare_surveys(self, survey1: ProcessingResult,
                       survey2: ProcessingResult) -> Dict:
        """Compare two surveys for change detection"""
        # Load point clouds
        pc1 = self.extract_point_cloud(survey1.point_cloud_path)
        pc2 = self.extract_point_cloud(survey2.point_cloud_path)

        # Calculate elevation differences
        # In production, use proper point cloud registration and comparison

        comparison = {
            'survey1_date': survey1.statistics.get('flight_date'),
            'survey2_date': survey2.statistics.get('flight_date'),
            'point_count_diff': len(pc2.points) - len(pc1.points),
            'changes_detected': []
        }

        return comparison

Volume Calculation Engine

from scipy.spatial import Delaunay
from scipy.interpolate import griddata
import numpy as np

class VolumeCalculator:
    """Advanced volume calculations from drone data"""

    def __init__(self, point_cloud: PointCloud):
        self.points = point_cloud.points
        self.colors = point_cloud.colors

    def calculate_cut_fill(self, design_surface: np.ndarray,
                          grid_size: float = 0.5) -> Dict:
        """Calculate cut and fill volumes compared to design surface"""
        # Create grid
        x_min, x_max = np.min(self.points[:, 0]), np.max(self.points[:, 0])
        y_min, y_max = np.min(self.points[:, 1]), np.max(self.points[:, 1])

        x_grid = np.arange(x_min, x_max, grid_size)
        y_grid = np.arange(y_