drone-site-survey

name	drone-site-survey
description	Process drone survey data for construction sites. Generate orthomosaics, DEMs, point clouds, calculate volumes, track progress, and integrate with BIM models for comparison.
homepage	https://datadrivenconstruction.io
metadata	{"openclaw":{"emoji":"🚀","os":["darwin","linux","win32"],"homepage":"https://datadrivenconstruction.io","requires":{"bins":"[Truncated]"}}}

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:  =   

     ():
         .gsd ==   .camera:
            
            sensor_width = .camera.sensor_width_mm
            focal_length = .camera.focal_length_mm
            image_width = .camera.image_width_px
            altitude = .flight_altitude

            .gsd = (altitude * sensor_width) / (focal_length * image_width) *   


 :
    orthomosaic_path: [] = 
    dem_path: [] = 
    dsm_path: [] = 
    point_cloud_path: [] = 
    report_path: [] = 
    statistics:  = field(default_factory=)

 :
    

     ():
        .output_dir = Path(output_dir)
        .output_dir.mkdir(parents=, exist_ok=)

     () -> ProcessingResult:
        
        result = ProcessingResult()
        result.statistics[] = flight.flight_id
        result.statistics[] = (flight.images)
        result.statistics[] = flight.gsd
        result.statistics[] = flight.date.isoformat()

        
        

         generate_ortho:
            result.orthomosaic_path = (.output_dir / )
            result.statistics[] = flight.gsd

         generate_dem:
            result.dem_path = (.output_dir / )
            result.dsm_path = (.output_dir / )

         generate_pointcloud:
            result.point_cloud_path = (.output_dir / )

        
        result.report_path = (.output_dir / )
         (result.report_path, )  f:
            json.dump(result.statistics, f, indent=)

         result

     () -> PointCloud:
        
        
        
        n_points = 
        points = np.random.rand(n_points, ) * [, , ]
        colors = np.random.randint(, , (n_points, ), dtype=np.uint8)

         PointCloud(points=points, colors=colors)

     () -> :
        
        
        pc1 = .extract_point_cloud(survey1.point_cloud_path)
        pc2 = .extract_point_cloud(survey2.point_cloud_path)

        
        

        comparison = {
            : survey1.statistics.get(),
            : survey2.statistics.get(),
            : (pc2.points) - (pc1.points),
            : []
        }

         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[:, :],
            design_surface[:, ],
            (xx, yy),
            method=
        )

        
        diff = actual_z - design_z
        cell_area = grid_size ** 

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

         {
            : (cut_volume),
            : (fill_volume),
            : (net_volume),
            :   net_volume >   ,
            : grid_size,
            : ((x_max - x_min) * (y_max - y_min))
        }

     () -> VolumeResult:
        
        
         scipy.spatial  ConvexHull
        hull = ConvexHull(.points[:, :])
        boundary_points = .points[hull.vertices]

         base_method == :
            reference_z = np.(boundary_points[:, ])
         base_method == :
            reference_z = np.mean(boundary_points[:, ])
         base_method == :
            
            reference_z = np.(.points[:, ])
        :
            reference_z = np.(.points[:, ])

        
        volume = ._triangulated_volume(reference_z)

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

         VolumeResult(
            volume_m3=volume,
            area_m2=hull_area,
            method= + base_method,
            reference_plane=,
            confidence=
        )

     () -> :
        
        
        points_2d = .points[:, :]
        tri = Delaunay(points_2d)

        volume = 
         simplex  tri.simplices:
            
            p1 = .points[simplex[]]
            p2 = .points[simplex[]]
            p3 = .points[simplex[]]

            
            avg_height = (p1[] + p2[] + p3[]) /  - reference_z
             avg_height > :
                
                area =  * (
                    (p2[] - p1[]) * (p3[] - p1[]) -
                    (p3[] - p1[]) * (p2[] - p1[])
                )
                volume += area * avg_height

         volume

     () -> :
        
        hull = ConvexHull(points_2d)
         hull.volume  

     () -> []:
        
        z_min = np.(.points[:, ])
        z_max = np.(.points[:, ])

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

        contours = []
         level  contour_levels:
            
            contours.append({
                : (level),
                :   level %  ==   
            })

         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({
                : survey.date.isoformat(),
                : i + ,
                : cut_fill[],
                : cut_fill[],
                : cut_fill[]
            })

        
         (progress) > :
            initial = progress[]
             p  progress[:]:
                 initial[] != :
                    p[] = (
                        (initial[] - p[]) /
                        (initial[]) * 
                    )
                :
                    p[] = 

         progress

     () -> :
        
        pc1 = .surveys[survey1_idx].point_cloud
        pc2 = .surveys[survey2_idx].point_cloud

        
        grid_size =   

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

        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 = {
            : [],  
            : [],  
            : 
        }

         i  ((x_bins) - ):
             j  ((y_bins) - ):
                
                cell_x = (x_bins[i] + x_bins[i + ]) / 
                cell_y = (y_bins[j] + y_bins[j + ]) / 

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

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

                     diff > threshold_m:
                        changes[].append({
                            : cell_x, : cell_y, : diff
                        })
                     diff < -threshold_m:
                        changes[].append({
                            : cell_x, : cell_y, : diff
                        })
                    :
                        changes[] += 

        changes[] = {
            : (changes[]),
            : (changes[]),
            : changes[],
            : .surveys[survey1_idx].date.isoformat(),
            : .surveys[survey2_idx].date.isoformat()
        }

         changes

     () -> :
        
         pandas  pd

        report_data = {
            : .project_name,
            : (.surveys),
            : {
                : .surveys[].date.isoformat()  .surveys  ,
                : .surveys[-].date.isoformat()  .surveys  
            }
        }

         design_surface   :
            report_data[] = .calculate_earthwork_progress(design_surface)

         pd.ExcelWriter(output_path, engine=)  writer:
            
            pd.DataFrame([report_data]).to_excel(writer, sheet_name=, index=)

            
            survey_data = [{
                : s.date,
                : s.orthomosaic_path
            }  s  .surveys]
            pd.DataFrame(survey_data).to_excel(writer, sheet_name=, index=)

            
               report_data:
                pd.DataFrame(report_data[]).to_excel(
                    writer, sheet_name=, index=
                )

         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(.bim[:, :], .bim[:, ], (xx, yy), method=)
        drone_z = griddata(
            .drone.points[:, :],
            .drone.points[:, ],
            (xx, yy),
            method=
        )

        
        diff = drone_z - bim_z

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

         {
            : (np.mean(valid_diff)),
            : (np.std(valid_diff)),
            : (np.(valid_diff)),
            : (np.(valid_diff)),
            : (np.sqrt(np.mean(valid_diff ** ))),
            : (np.(np.(valid_diff) < ) / (valid_diff) * ),
            : (np.(np.(valid_diff) < ) / (valid_diff) * ),
            : ((valid_diff))
        }

     () -> []:
        
        deviations = []
        grid_size = 

        
        x_min = (np.(.bim[:, ]), np.(.drone.points[:, ]))
        x_max = (np.(.bim[:, ]), np.(.drone.points[:, ]))
        y_min = (np.(.bim[:, ]), np.(.drone.points[:, ]))
        y_max = (np.(.bim[:, ]), np.(.drone.points[:, ]))

         x  np.arange(x_min, x_max, grid_size):
             y  np.arange(y_min, y_max, grid_size):
                
                bim_mask = (
                    (.bim[:, ] >= x) & (.bim[:, ] < x + grid_size) &
                    (.bim[:, ] >= y) & (.bim[:, ] < y + grid_size)
                )
                drone_mask = (
                    (.drone.points[:, ] >= x) & (.drone.points[:, ] < x + grid_size) &
                    (.drone.points[:, ] >= y) & (.drone.points[:, ] < y + grid_size)
                )

                 np.(bim_mask) >   np.(drone_mask) > :
                    bim_z = np.mean(.bim[bim_mask, ])
                    drone_z = np.mean(.drone.points[drone_mask, ])
                    diff = drone_z - bim_z

                     (diff) > tolerance_m:
                        deviations.append({
                            : x + grid_size / ,
                            : y + grid_size / ,
                            : bim_z,
                            : drone_z,
                            : diff,
                            :   diff >   
                        })

         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