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_min, y_max, grid_size)
        xx, yy = np.meshgrid(x_grid, y_grid)

        # Interpolate actual surface
        actual_z = griddata(
            self.points[:, :2],
            self.points[:, 2],
            (xx, yy),
            method='linear'
        )

        # Interpolate design surface
        design_z = griddata(
            design_surface[:, :2],
            design_surface[:, 2],
            (xx, yy),
            method='linear'
        )

        # Calculate differences
        diff = actual_z - design_z
        cell_area = grid_size ** 2

        cut_volume = np.nansum(diff[diff > 0]) * cell_area
        fill_volume = np.nansum(np.abs(diff[diff < 0])) * cell_area
        net_volume = cut_volume - fill_volume

        return {
            'cut_volume_m3': float(cut_volume),
            'fill_volume_m3': float(fill_volume),
            'net_volume_m3': float(net_volume),
            'balance': 'cut' if net_volume > 0 else 'fill',
            'grid_size_m': grid_size,
            'area_m2': float((x_max - x_min) * (y_max - y_min))
        }

    def calculate_stockpile_volume(self, base_method: str = 'lowest_perimeter') -> VolumeResult:
        """Calculate stockpile volume"""
        # Find boundary points
        from scipy.spatial import ConvexHull
        hull = ConvexHull(self.points[:, :2])
        boundary_points = self.points[hull.vertices]

        if base_method == 'lowest_perimeter':
            reference_z = np.min(boundary_points[:, 2])
        elif base_method == 'average_perimeter':
            reference_z = np.mean(boundary_points[:, 2])
        elif base_method == 'triangulated':
            # Create base triangulation from boundary
            reference_z = np.min(self.points[:, 2])
        else:
            reference_z = np.min(self.points[:, 2])

        # Calculate volume using triangulated surface
        volume = self._triangulated_volume(reference_z)

        hull_area = self._calculate_hull_area(boundary_points[:, :2])

        return VolumeResult(
            volume_m3=volume,
            area_m2=hull_area,
            method='triangulated_' + base_method,
            reference_plane=f'z={reference_z:.2f}',
            confidence=0.95
        )

    def _triangulated_volume(self, reference_z: float) -> float:
        """Calculate volume using Delaunay triangulation"""
        # Create 2D triangulation
        points_2d = self.points[:, :2]
        tri = Delaunay(points_2d)

        volume = 0
        for simplex in tri.simplices:
            # Get triangle vertices
            p1 = self.points[simplex[0]]
            p2 = self.points[simplex[1]]
            p3 = self.points[simplex[2]]

            # Calculate prism volume
            avg_height = (p1[2] + p2[2] + p3[2]) / 3 - reference_z
            if avg_height > 0:
                # Triangle area
                area = 0.5 * abs(
                    (p2[0] - p1[0]) * (p3[1] - p1[1]) -
                    (p3[0] - p1[0]) * (p2[1] - p1[1])
                )
                volume += area * avg_height

        return volume

    def _calculate_hull_area(self, points_2d: np.ndarray) -> float:
        """Calculate area of convex hull"""
        hull = ConvexHull(points_2d)
        return hull.volume  # In 2D, volume is actually area

    def generate_contours(self, interval: float = 1.0) -> List[Dict]:
        """Generate elevation contours"""
        z_min = np.min(self.points[:, 2])
        z_max = np.max(self.points[:, 2])

        contour_levels = np.arange(
            np.floor(z_min / interval) * interval,
            np.ceil(z_max / interval) * interval + interval,
            interval
        )

        contours = []
        for level in contour_levels:
            # In production, use proper contouring algorithm
            contours.append({
                'elevation': float(level),
                'type': 'major' if level % 5 == 0 else 'minor'
            })

        return contours

Progress Monitoring

from datetime import date
from typing import List, Dict

@dataclass
class ProgressPoint:
    date: date
    point_cloud: PointCloud
    orthomosaic_path: str
    annotations: Dict = field(default_factory=dict)

class ConstructionProgressMonitor:
    """Monitor construction progress from drone surveys"""

    def __init__(self, project_name: str):
        self.project_name = project_name
        self.surveys: List[ProgressPoint] = []
        self.volume_calculator = None

    def add_survey(self, survey_date: date, point_cloud: PointCloud,
                  orthomosaic_path: str):
        """Add survey for progress tracking"""
        self.surveys.append(ProgressPoint(
            date=survey_date,
            point_cloud=point_cloud,
            orthomosaic_path=orthomosaic_path
        ))
        self.surveys.sort(key=lambda x: x.date)

    def calculate_earthwork_progress(self, design_surface: np.ndarray) -> List[Dict]:
        """Calculate earthwork progress over time"""
        progress = []

        for i, survey in enumerate(self.surveys):
            calc = VolumeCalculator(survey.point_cloud)
            cut_fill = calc.calculate_cut_fill(design_surface)

            progress.append({
                'date': survey.date.isoformat(),
                'survey_index': i + 1,
                'cut_volume_m3': cut_fill['cut_volume_m3'],
                'fill_volume_m3': cut_fill['fill_volume_m3'],
                'net_volume_m3': cut_fill['net_volume_m3']
            })

        # Calculate progress percentages
        if len(progress) > 1:
            initial = progress[0]
            for p in progress[1:]:
                if initial['net_volume_m3'] != 0:
                    p['progress_pct'] = (
                        (initial['net_volume_m3'] - p['net_volume_m3']) /
                        abs(initial['net_volume_m3']) * 100
                    )
                else:
                    p['progress_pct'] = 0

        return progress

    def detect_changes(self, survey1_idx: int, survey2_idx: int,
                      threshold_m: float = 0.5) -> Dict:
        """Detect changes between two surveys"""
        pc1 = self.surveys[survey1_idx].point_cloud
        pc2 = self.surveys[survey2_idx].point_cloud

        # Create grids for comparison
        grid_size = 1.0  # 1m grid

        x_min = min(np.min(pc1.points[:, 0]), np.min(pc2.points[:, 0]))
        x_max = max(np.max(pc1.points[:, 0]), np.max(pc2.points[:, 0]))
        y_min = min(np.min(pc1.points[:, 1]), np.min(pc2.points[:, 1]))
        y_max = max(np.max(pc1.points[:, 1]), np.max(pc2.points[:, 1]))

        x_bins = np.arange(x_min, x_max + grid_size, grid_size)
        y_bins = np.arange(y_min, y_max + grid_size, grid_size)

        changes = {
            'increased': [],  # Areas where elevation increased
            'decreased': [],  # Areas where elevation decreased
            'unchanged': 0
        }

        for i in range(len(x_bins) - 1):
            for j in range(len(y_bins) - 1):
                # Get points in cell for both surveys
                cell_x = (x_bins[i] + x_bins[i + 1]) / 2
                cell_y = (y_bins[j] + y_bins[j + 1]) / 2

                mask1 = (
                    (pc1.points[:, 0] >= x_bins[i]) & (pc1.points[:, 0] < x_bins[i + 1]) &
                    (pc1.points[:, 1] >= y_bins[j]) & (pc1.points[:, 1] < y_bins[j + 1])
                )
                mask2 = (
                    (pc2.points[:, 0] >= x_bins[i]) & (pc2.points[:, 0] < x_bins[i + 1]) &
                    (pc2.points[:, 1] >= y_bins[j]) & (pc2.points[:, 1] < y_bins[j + 1])
                )

                if np.sum(mask1) > 0 and np.sum(mask2) > 0:
                    z1 = np.mean(pc1.points[mask1, 2])
                    z2 = np.mean(pc2.points[mask2, 2])
                    diff = z2 - z1

                    if diff > threshold_m:
                        changes['increased'].append({
                            'x': cell_x, 'y': cell_y, 'change_m': diff
                        })
                    elif diff < -threshold_m:
                        changes['decreased'].append({
                            'x': cell_x, 'y': cell_y, 'change_m': diff
                        })
                    else:
                        changes['unchanged'] += 1

        changes['summary'] = {
            'increased_cells': len(changes['increased']),
            'decreased_cells': len(changes['decreased']),
            'unchanged_cells': changes['unchanged'],
            'date_from': self.surveys[survey1_idx].date.isoformat(),
            'date_to': self.surveys[survey2_idx].date.isoformat()
        }

        return changes

    def generate_progress_report(self, output_path: str,
                                design_surface: np.ndarray = None) -> str:
        """Generate progress report"""
        import pandas as pd

        report_data = {
            'project': self.project_name,
            'surveys': len(self.surveys),
            'date_range': {
                'start': self.surveys[0].date.isoformat() if self.surveys else None,
                'end': self.surveys[-1].date.isoformat() if self.surveys else None
            }
        }

        if design_surface is not None:
            report_data['earthwork_progress'] = self.calculate_earthwork_progress(design_surface)

        with pd.ExcelWriter(output_path, engine='openpyxl') as writer:
            # Summary
            pd.DataFrame([report_data]).to_excel(writer, sheet_name='Summary', index=False)

            # Survey list
            survey_data = [{
                'Date': s.date,
                'Orthomosaic': s.orthomosaic_path
            } for s in self.surveys]
            pd.DataFrame(survey_data).to_excel(writer, sheet_name='Surveys', index=False)

            # Progress (if available)
            if 'earthwork_progress' in report_data:
                pd.DataFrame(report_data['earthwork_progress']).to_excel(
                    writer, sheet_name='Progress', index=False
                )

        return output_path

BIM Comparison

class BIMDroneComparator:
    """Compare drone survey with BIM model"""

    def __init__(self, bim_surface: np.ndarray, drone_pointcloud: PointCloud):
        self.bim = bim_surface
        self.drone = drone_pointcloud

    def compare_elevations(self, grid_size: float = 1.0) -> Dict:
        """Compare drone elevations with BIM design"""
        # Create comparison grid
        x_min = max(np.min(self.bim[:, 0]), np.min(self.drone.points[:, 0]))
        x_max = min(np.max(self.bim[:, 0]), np.max(self.drone.points[:, 0]))
        y_min = max(np.min(self.bim[:, 1]), np.min(self.drone.points[:, 1]))
        y_max = min(np.max(self.bim[:, 1]), np.max(self.drone.points[:, 1]))

        x_grid = np.arange(x_min, x_max, grid_size)
        y_grid = np.arange(y_min, y_max, grid_size)
        xx, yy = np.meshgrid(x_grid, y_grid)

        # Interpolate both surfaces
        bim_z = griddata(self.bim[:, :2], self.bim[:, 2], (xx, yy), method='linear')
        drone_z = griddata(
            self.drone.points[:, :2],
            self.drone.points[:, 2],
            (xx, yy),
            method='linear'
        )

        # Calculate differences
        diff = drone_z - bim_z

        valid_mask = ~np.isnan(diff)
        valid_diff = diff[valid_mask]

        return {
            'mean_diff_m': float(np.mean(valid_diff)),
            'std_diff_m': float(np.std(valid_diff)),
            'max_diff_m': float(np.max(valid_diff)),
            'min_diff_m': float(np.min(valid_diff)),
            'rmse_m': float(np.sqrt(np.mean(valid_diff ** 2))),
            'within_5cm': float(np.sum(np.abs(valid_diff) < 0.05) / len(valid_diff) * 100),
            'within_10cm': float(np.sum(np.abs(valid_diff) < 0.10) / len(valid_diff) * 100),
            'comparison_points': int(len(valid_diff))
        }

    def find_deviations(self, tolerance_m: float = 0.1) -> List[Dict]:
        """Find areas with significant deviations"""
        deviations = []
        grid_size = 1.0

        # Grid comparison
        x_min = max(np.min(self.bim[:, 0]), np.min(self.drone.points[:, 0]))
        x_max = min(np.max(self.bim[:, 0]), np.max(self.drone.points[:, 0]))
        y_min = max(np.min(self.bim[:, 1]), np.min(self.drone.points[:, 1]))
        y_max = min(np.max(self.bim[:, 1]), np.max(self.drone.points[:, 1]))

        for x in np.arange(x_min, x_max, grid_size):
            for y in np.arange(y_min, y_max, grid_size):
                # Get average elevations
                bim_mask = (
                    (self.bim[:, 0] >= x) & (self.bim[:, 0] < x + grid_size) &
                    (self.bim[:, 1] >= y) & (self.bim[:, 1] < y + grid_size)
                )
                drone_mask = (
                    (self.drone.points[:, 0] >= x) & (self.drone.points[:, 0] < x + grid_size) &
                    (self.drone.points[:, 1] >= y) & (self.drone.points[:, 1] < y + grid_size)
                )

                if np.sum(bim_mask) > 0 and np.sum(drone_mask) > 0:
                    bim_z = np.mean(self.bim[bim_mask, 2])
                    drone_z = np.mean(self.drone.points[drone_mask, 2])
                    diff = drone_z - bim_z

                    if abs(diff) > tolerance_m:
                        deviations.append({
                            'x': x + grid_size / 2,
                            'y': y + grid_size / 2,
                            'bim_elevation': bim_z,
                            'actual_elevation': drone_z,
                            'deviation_m': diff,
                            'type': 'high' if diff > 0 else 'low'
                        })

        return deviations

Quick Reference

Measurement	Method	Accuracy
Stockpile Volume	Triangulated	±2-5%
Cut/Fill Volume	Grid comparison	±5%
Area Measurement	Orthomosaic	<1cm GSD
Elevation (DEM)	Photogrammetry	±2-5cm
Progress Tracking	Multi-temporal	Relative

Resources

• OpenDroneMap: https://www.opendronemap.org
• Pix4D: https://www.pix4d.com
• LAStools: https://rapidlasso.com/lastools/
• DDC Website: https://datadrivenconstruction.io

Next Steps

• See progress-monitoring-cv for image-based progress
• See bim-validation-pipeline for model comparison
• See data-visualization for 3D visualization