Drone Orthomosaic Analysis and Agricultural Field Phenotyping in Python
Dronelytics is a comprehensive Python package for processing multispectral drone imagery and LAS/LAZ point clouds. It provides production-ready tools for vegetation analysis, automated plot detection, spectral data extraction, digital elevation modeling, and 3D visualization.
Built for researchers, agronomists, and GIS professionals who need fast, reliable analysis of drone data.
- Precision Agriculture - Crop health monitoring and variable-rate application mapping
- Crop Phenotyping - Automated trait extraction from experimental plots
- Plant Breeding - Large-scale genotype evaluation and selection
- Environmental Monitoring - Wetland and vegetation change detection
- Crop Insurance - Automated damage assessment and loss evaluation
- Remote Sensing Research - Multispectral algorithm development and validation
- β RGB, RGB+NIR, and RGB+NIR+Red-Edge imagery support
- β GeoTIFF-based workflows with metadata preservation
- β Flexible multi-band configuration
- β Efficient memory management for large files
10+ vegetation indices including:
- NDVI - Normalized Difference Vegetation Index (general vegetation health)
- NDRE - Normalized Difference Red Edge Index (crop stress detection)
- GNDVI - Green NDVI (chlorophyll content)
- VARI - Visible Atmospherically Resistant Index (early season analysis)
- EVI - Enhanced Vegetation Index (canopy structure)
- Plus 5 more: SAVI, MSAVI, ARVI, CVI, OSAVI
- Boundary detection from orthomosaics
- Three segmentation methods (watershed, quickshift, felzenszwalb)
- Plot-level statistics computation
- GeoJSON export for GIS integration
- Pixel-level vegetation index values
- Per-plot summary statistics (mean, std, min, max, median)
- Integration of multiple indices
- CSV and Excel export
- Digital Terrain Model (DTM) generation
- Digital Surface Model (DSM) generation
- Canopy Height Model (CHM) calculation
- 3D mesh generation from LAS/LAZ files
- GeoTIFF and STL export formats
- CSV (pixel-level data)
- Excel with multiple sheets (pixel data + statistics + metadata)
- GeoJSON (vector boundaries)
- GeoTIFF (raster models)
- STL (3D mesh models)
pip install dronelyticspip install dronelytics[pointcloud]git clone https://github.com/Lalitgis/dronelytics.git
cd dronelytics
pip install -e ".[pointcloud]"python -c "from dronelytics import Orthomosaic, VegetationIndices, PlotSegmentation; print('β Dronelytics installed successfully')"from dronelytics import Orthomosaic, VegetationIndices
# Load orthomosaic with band configuration
band_config = {
'red': 1,
'green': 2,
'blue': 3,
'nir': 4,
'red_edge': 5 # Optional: for 5-band imagery
}
ortho = Orthomosaic('field_image.tif', band_config=band_config)
# Calculate vegetation indices
vi = VegetationIndices(ortho)
ndvi = vi.calculate_ndvi()
print(f"NDVI Mean: {ndvi.mean:.4f} Β± {ndvi.std:.4f}")
print(f"Valid Pixels: {ndvi.pixel_count}")
# Other indices
ndre = vi.calculate_ndre()
gndvi = vi.calculate_gndvi()
vari = vi.calculate_vari()from dronelytics import PlotSegmentation
segmentation = PlotSegmentation(ortho)
# Detect plots using watershed method (recommended for NDVI)
result = segmentation.segment_by_ndvi(
method='watershed',
min_plot_size=100
)
print(f"Detected {result.num_plots} plots")
print(f"Boundaries: {len(result.boundaries)}")
# Export as GeoJSON
gdf = segmentation.get_boundaries_geodataframe()
gdf.to_file('plot_boundaries.geojson', driver='GeoJSON')from dronelytics import PixelExtraction
extraction = PixelExtraction(ortho)
# Extract vegetation indices for each plot
result = extraction.extract_from_segmentation(
segmentation_result=result,
vegetation_indices={
'ndvi': ndvi.values,
'ndre': ndre.values,
'gndvi': gndvi.values
},
raw_bands=['red', 'nir']
)
# Get summary statistics
summary = extraction.summarize_by_plot()
summary.to_csv('plot_statistics.csv', index=False)
# Export to Excel with multiple sheets
from dronelytics import ExcelExporter
exporter = ExcelExporter()
exporter.export(
output_path='analysis_results.xlsx',
pixel_data=extraction.to_dataframe(),
plot_stats=extraction.get_plot_statistics(),
vegetation_indices={'ndvi': ndvi, 'ndre': ndre}
)from dronelytics import PointCloudProcessor
pc = PointCloudProcessor('lidar_data.las')
# Generate elevation models
dtm, dtm_meta = pc.generate_dem(resolution=0.1, method='ground')
dsm, dsm_meta = pc.generate_dem(resolution=0.1, method='all')
# Generate canopy height model
chm, chm_meta = pc.generate_chm(resolution=0.1)
# Export as GeoTIFF
pc.export_dem_as_tiff(dtm, 'dtm.tif', resolution=0.1)
pc.export_dem_as_tiff(chm, 'chm.tif', resolution=0.1)
# Create 3D mesh
mesh = pc.create_3d_mesh(point_type='vegetation')
pc.export_mesh_as_stl(mesh, 'canopy_3d.stl')The band_config parameter maps band names to their position in your image.
Multispectral Drone (4-band: RGB+NIR)
band_config = {
'red': 1,
'green': 2,
'blue': 3,
'nir': 4
}Multispectral Drone (5-band: RGB+NIR+Red-Edge)
band_config = {
'red': 1,
'green': 2,
'blue': 3,
'nir': 4,
'red_edge': 5 # Use 'red_edge' (with underscore)
}Satellite Imagery (Landsat 8)
band_config = {
'blue': 2,
'green': 3,
'red': 4,
'nir': 5,
'swir_1': 6,
'swir_2': 7
}# Using GDAL (command line):
gdalinfo your_file.tif | grep "Band"
# Using Python:
import rasterio
with rasterio.open('your_file.tif') as src:
print(f"Number of bands: {src.count}")
for i in range(1, src.count + 1):
print(f"Band {i}: {src.descriptions[i-1]}")| Index | Formula | Best For | Range |
|---|---|---|---|
| NDVI | (NIR - RED) / (NIR + RED) | General vegetation health | -1 to +1 |
| NDRE | (NIR - RedEdge) / (NIR + RedEdge) | Crop stress detection | -1 to +1 |
| GNDVI | (NIR - GREEN) / (NIR + GREEN) | Chlorophyll content | -1 to +1 |
| VARI | (GREEN - RED) / (GREEN + RED - BLUE) | Early season analysis | -1 to +1 |
| EVI | 2.5 * (NIR - RED) / (NIR + 6RED - 7.5BLUE + 1) | Canopy structure | -1 to +1 |
| SAVI | 1.5 * (NIR - RED) / (NIR + RED + 0.5) | Soil-adjusted | -1 to +1 |
| Method | Input | Best For | Speed | Quality |
|---|---|---|---|---|
| watershed | NDVI/index | Regular, distinct plots | Fast | Excellent |
| felzenszwalb | NDVI/index | Complex, irregular boundaries | Fast | Very Good |
| quickshift | RGB image (3-band) | Fine-scale details | Slow | Excellent |
Use when you have RGB data and want fine-scale segmentation:
# First, prepare RGB data (not NDVI)
import numpy as np
rgb = np.stack([ortho.get_band('red'),
ortho.get_band('green'),
ortho.get_band('blue')], axis=2)
# Normalize to 0-1 range
rgb = (rgb - rgb.min()) / (rgb.max() - rgb.min())
# Use custom implementation
from skimage import segmentation
labeled = segmentation.quickshift(rgb, kernel_size=15, max_dist=40)Use when you have NDVI or other vegetation indices:
result = segmentation.segment_by_ndvi(
method='watershed',
min_plot_size=100 # Minimum pixels per plot
)Use for varying plot sizes or irregular boundaries:
result = segmentation.segment_by_ndvi(
method='felzenszwalb',
scale=500, # Larger = fewer, larger regions
sigma=0.5 # Smoothing parameter
)import logging
from dronelytics import (
Orthomosaic, VegetationIndices, PlotSegmentation,
PixelExtraction, ExcelExporter, setup_logger
)
# Setup logging
logging.basicConfig(level=logging.INFO)
logger = setup_logger('field_analysis')
try:
# 1. Load orthomosaic
logger.info("Loading orthomosaic...")
ortho = Orthomosaic('field_image.tif', band_config={
'red': 1,
'green': 2,
'blue': 3,
'nir': 4,
'red_edge': 5
})
# 2. Calculate vegetation indices
logger.info("Calculating vegetation indices...")
vi = VegetationIndices(ortho)
ndvi = vi.calculate_ndvi()
ndre = vi.calculate_ndre()
gndvi = vi.calculate_gndvi()
logger.info(f"NDVI: Mean={ndvi.mean:.3f}, Std={ndvi.std:.3f}")
# 3. Detect plot boundaries
logger.info("Detecting plot boundaries...")
segmentation = PlotSegmentation(ortho)
seg_result = segmentation.segment_by_ndvi(method='watershed', min_plot_size=100)
logger.info(f"Detected {seg_result.num_plots} plots")
# 4. Extract spectral data
logger.info("Extracting pixel values...")
extraction = PixelExtraction(ortho)
ext_result = extraction.extract_from_segmentation(
segmentation_result=seg_result,
vegetation_indices={'ndvi': ndvi.values, 'ndre': ndre.values, 'gndvi': gndvi.values},
raw_bands=['red', 'nir']
)
logger.info(f"Extracted {len(ext_result.pixels)} pixels")
# 5. Generate reports
logger.info("Generating reports...")
extraction.to_csv('pixel_data.csv')
summary = extraction.summarize_by_plot()
summary.to_csv('plot_summary.csv', index=False)
exporter = ExcelExporter()
exporter.export(
output_path='field_analysis.xlsx',
pixel_data=extraction.to_dataframe(),
plot_stats=extraction.get_plot_statistics(),
vegetation_indices={'ndvi': ndvi, 'ndre': ndre, 'gndvi': gndvi},
orthomosaic_metadata=ortho.metadata
)
# 6. Export boundaries
logger.info("Exporting boundaries...")
gdf = segmentation.get_boundaries_geodataframe()
gdf.to_file('plot_boundaries.geojson', driver='GeoJSON')
logger.info("β Analysis complete!")
except Exception as e:
logger.error(f"Analysis failed: {e}", exc_info=True)
finally:
ortho.close()For files larger than 2GB, consider tiling:
# Process in tiles to manage memory
ortho = Orthomosaic('large_file.tif', band_config=band_config)
try:
# Your processing code here
...
finally:
ortho.clear_cache() # Free memory
ortho.close()import glob
from pathlib import Path
tif_files = glob.glob('data/*.tif')
for tif_file in tif_files:
logger.info(f"Processing {Path(tif_file).stem}...")
ortho = Orthomosaic(tif_file, band_config=band_config)
vi = VegetationIndices(ortho)
ndvi = vi.calculate_ndvi()
# Export results
output_file = f"results/{Path(tif_file).stem}_ndvi.tif"
# ... save results
ortho.close()# Define your own formula
def custom_index(ortho):
red = ortho.get_band('red').astype(float)
nir = ortho.get_band('nir').astype(float)
# Your custom formula
index = (nir - red) / (nir + red + 1e-8)
return index
custom = custom_index(ortho)Solution: Install with point cloud support:
pip install dronelytics[pointcloud]
# or manually:
pip install scikit-image geopandas shapely openpyxlProblem: Band named 'rededge' but code expects 'red_edge'
Solution: Use correct band name with underscore:
band_config = {
# ...
'red_edge': 5 # Use 'red_edge' (with underscore)
}Problem: Trying to use quickshift method with NDVI data (1-band)
Solution: Use watershed instead or prepare RGB data:
# Option 1 (Recommended):
result = segmentation.segment_by_ndvi(method='watershed')
# Option 2: Use RGB data with quickshift
rgb = np.stack([ortho.get_band('red'),
ortho.get_band('green'),
ortho.get_band('blue')], axis=2)Solution 1: Use smaller tile size
ortho = Orthomosaic('file.tif', band_config=band_config)
ortho.clear_cache() # Free memory between operationsSolution 2: Process in tiles or use a machine with more RAM
Check:
from pathlib import Path
file = Path('your_file.tif')
print(f"File exists: {file.exists()}")
print(f"Absolute path: {file.absolute()}")Orthomosaic
ortho = Orthomosaic(filepath, band_config)
ortho.get_band(name) # Get single band as numpy array
ortho.close() # Free resources
ortho.clear_cache() # Clear cached dataVegetationIndices
vi = VegetationIndices(ortho)
ndvi = vi.calculate_ndvi() # Returns IndexResult
ndre = vi.calculate_ndre()
gndvi = vi.calculate_gndvi()
# ... and morePlotSegmentation
seg = PlotSegmentation(ortho)
result = seg.segment_by_ndvi(method='watershed')
gdf = seg.get_boundaries_geodataframe() # Export to GeoJSONPixelExtraction
ext = PixelExtraction(ortho)
result = ext.extract_from_segmentation(segmentation_result, vegetation_indices)
df = ext.to_dataframe()
ext.to_csv('output.csv')PointCloudProcessor
pc = PointCloudProcessor('data.las')
dtm, dtm_meta = pc.generate_dem(resolution=0.1, method='ground')
mesh = pc.create_3d_mesh()
pc.export_mesh_as_stl(mesh, 'output.stl')| Format | Use Case | Function |
|---|---|---|
| CSV | Pixel-level data, easy to analyze | extraction.to_csv() |
| Excel | Summary reports, multi-sheet | ExcelExporter.export() |
| GeoJSON | GIS integration, vector boundaries | gdf.to_file(..., 'GeoJSON') |
| GeoTIFF | Raster models, elevation data | pc.export_dem_as_tiff() |
| STL | 3D visualization, 3D printing | pc.export_mesh_as_stl() |
- numpy (numerical computing)
- pandas (data manipulation)
- rasterio (geospatial I/O)
- scipy (scientific computing)
- matplotlib (visualization)
- scikit-image (image processing algorithms)
- geopandas (vector data handling)
- shapely (geometry operations)
- openpyxl (Excel file handling)
- laspy (LAS/LAZ file reading)
- pyvista (3D visualization)
- Always define band_config explicitly - Don't rely on defaults
- Use watershed for NDVI - Quickshift requires RGB preprocessing
- Test on small files first - Before processing large datasets
- Handle exceptions - Use try/finally to ensure ortho.close() is called
- Enable logging - Helps debug issues and monitor progress
- Export results - Keep original data and generated indices separate
Contributions are welcome! Please:
- Fork the repository
- Create a feature branch
- Add tests for new functionality
- Submit a pull request
MIT License - see LICENSE file for details
- Documentation: https://github.com/Lalitgis/dronelytics/wiki
- Issues: https://github.com/Lalitgis/dronelytics/issues
- Email: lalitiaas@gmail.com
If you use Dronelytics in your research, please cite:
@software{dronelytics2026,
author = {BC, Lalit},
title = {Dronelytics: Comprehensive Drone Orthomosaic Analysis and Agricultural Field Phenotyping Toolkit},
year = {2026},
version = {1.0.3},
url = {https://github.com/Lalitgis/dronelytics}
}- GPU acceleration for large-scale processing
- Web interface for remote processing
- Machine learning model integration
- Hyperspectral image support
- Real-time drone feed processing
- Cloud storage integration (AWS S3, Azure)
Made with β€οΈ for agricultural researchers and GIS professionals