AMAPVox-Scripts

R scripts for automated voxelization and Plant Area Density (PAD) calculation from LiDAR data using AMAPVox.

Quick Start

Extract Plant Area Density (PAD) from your RiSCAN project in three simple steps.

⚡ Fast Preview Mode (Recommended for Testing)

Best for: Quick testing, parameter exploration, initial runs

# Fast preview run (10-15 minutes)
# Uses 10% of shots, 1m voxels
CONFIG_XML=configs/LSS_2021_07_core_preview.xml Rscript scripts/run_voxelization_batch.R

# Then compute PAD from preview voxels
VOX_PATH=data/voxel/LSS_2021_07_core_preview.vox Rscript scripts/get_voxel_pad.R

Speed improvements:

• 10% decimation (vs 50%) → 5x faster
• 1.0m voxels (vs 0.5m) → 8x fewer voxels
• Expected time: 10-15 minutes

📊 Full Resolution Mode

Best for: Final analysis, publication-quality results

# Step 1: Voxelize the RiSCAN project (full resolution, ~30-45 minutes)
Rscript scripts/run_voxelization_batch.R

# Step 2: Compute PAD from voxels
VOX_PATH=data/voxel/LSS_2021_07_core.vox Rscript scripts/get_voxel_pad.R

# Step 3 (optional): Analyze PAD data
VOX_FILE=data/voxel/voxel_with_pad.vox Rscript scripts/example_pad_analysis.R

Processing time: 30-45 minutes

🖱️ GUI Mode

Best for: First-time users, exploring parameters, visual feedback

# Launches GUI, then computes PAD after you save the voxel file
CREATE_VOX=true Rscript scripts/get_voxel_pad.R

📁 Existing Voxel File

Best for: When you already have a .vox file

VOX_PATH=path/to/your.vox Rscript scripts/get_voxel_pad.R

What's Included

✅ R scripts for voxelization and PAD computation ✅ XML configuration for RiSCAN projects ✅ DTM generation from LAS/LAZ point clouds ✅ Multi-threaded processing support ✅ Comprehensive documentation

What Each Script Does

run_voxelization_batch.R

• Reads RiSCAN project (up to 66 scan positions)
• Processes all scans with their transformations
• Computes transmittance and attenuation per voxel
• Outputs: .vox file (~100MB)

Processing time: 10-60 minutes (depends on point density and settings)

get_voxel_pad.R

• Reads voxel file
• Computes PAD from transmittance/attenuation
• Uses spherical leaf angle distribution
• Outputs: .vox file with PAD variables

Processing time: < 1 minute

create_dtm.R

• Generates Digital Terrain Model from LAS/LAZ point clouds
• Multi-threaded ground classification (PMF/CSF algorithms)
• Exports DTM in ESRI ASCII grid format
• See scripts/DTM_README.md for details

Usage:

THREADS=16 ALGORITHM=csf Rscript scripts/create_dtm.R input.las output.asc

example_pad_analysis.R

• Reads voxel file with PAD
• Computes Plant Area Index (PAI)
• Generates statistics and summaries
• Exports sample data to CSV

Processing time: < 1 minute

Configuration

Configuration Files

Two configurations are provided:

configs/LSS_2021_07_core_preview.xml - Fast preview

• 1.0m voxel resolution
• 10% decimation (10% of shots)
• ~5-10 minute processing time
• Good for testing and exploration

configs/LSS_2021_07_core.xml - Full resolution

• 0.5m voxel resolution
• 50% decimation (50% of shots)
• ~20-40 minute processing time (with 3 parallel tasks)
• Good for final analysis

Parallel Processing

Control CPU usage with environment variables:

# Recommended: Balanced performance (16 cores)
NUM_TASKS=2 THREADS_PER_TASK=8 Rscript scripts/run_voxelization_batch.R

# Conservative: Most stable (8 cores, default)
NUM_TASKS=1 THREADS_PER_TASK=8 Rscript scripts/run_voxelization_batch.R

# Aggressive: May be unstable with some datasets (24 cores)
NUM_TASKS=3 THREADS_PER_TASK=8 Rscript scripts/run_voxelization_batch.R

⚠️ Threading Limitations:

• The native RXP library has thread safety issues at very high parallelism
• Using >24 threads (e.g., NUM_TASKS=2 THREADS_PER_TASK=14) may cause JVM crashes
• Recommended max: NUM_TASKS=3 with THREADS_PER_TASK=8 (24 cores total)
• Safe default: NUM_TASKS=1 with THREADS_PER_TASK=8 (8 cores total)

Voxel Resolution

Edit config XML file to adjust resolution:

<!-- High resolution (0.5m voxels, slower): -->
<voxelspace ... resolution="(0.5, 0.5, 0.5)" split="(240, 240, 400)" ... />

<!-- Medium resolution (1.0m voxels, faster): -->
<voxelspace ... resolution="(1.0, 1.0, 1.0)" split="(120, 120, 200)" ... />

<!-- Low resolution (2.0m voxels, very fast): -->
<voxelspace ... resolution="(2.0, 2.0, 2.0)" split="(60, 60, 100)" ... />

Note: split must match spatial extent divided by resolution.

Decimation (Speed vs Accuracy)

<!-- Full data (slowest, most accurate): -->
<parameters decimation-factor="1.0" />

<!-- Half data (balanced): -->
<parameters decimation-factor="0.5" />

<!-- 10% data (fastest, good for preview): -->
<parameters decimation-factor="0.1" />

Output Variables

Only the minimum variables needed for PAD are enabled by default:

• number_of_shots - pulse count per voxel
• mean_angle - mean beam angle
• estimated_transmittance - transmittance estimate
• attenuation_fpl_unbiased_mle - Free Path Length attenuation
• attenuation_ppl_mle - Path Length attenuation
• ground_distance - height above ground (requires DTM)

Enable additional variables in configs/LSS_2021_07_core.xml as needed.

Output Files

Voxel file structure

data/voxel/
├── LSS_2021_07_core.vox       # Voxelization output
└── voxel_with_pad.vox         # With PAD variables added

PAD variables

After running get_voxel_pad.R, your voxel file will contain:

• pad_transmittance - PAD from transmittance
• pad_attenuation_FPL_unbiasedMLE - PAD from Free Path Length
• pad_attenuation_PPL_MLE - PAD from Path Length

All in units of m²/m³ (plant area per volume)

Troubleshooting

JVM Crash (SIGSEGV)

Symptom: Process crashes with "fatal error", "SIGSEGV", or creates hs_err_pid*.log

Cause: Native RXP library has thread safety issues at high parallelism (>24 threads)

Solution: Reduce parallel tasks:

# Use default safe settings (8 cores)
Rscript scripts/run_voxelization_batch.R

# Or specify explicitly
NUM_TASKS=1 THREADS_PER_TASK=8 Rscript scripts/run_voxelization_batch.R

Do NOT use: NUM_TASKS × THREADS_PER_TASK > 24

"Out of Memory" (OutOfMemoryError)

Solution: Increase Java heap size (default is 8GB):

# Edit scripts/run_voxelization_batch.R, line 96:
jvm.options = "-Xms16384m"  # 16GB instead of 8GB

Or reduce voxel resolution:

<!-- In XML config, change from 0.5m to 1.0m -->
<voxelspace ... resolution="(1.0, 1.0, 1.0)" split="(120, 120, 200)" ... />

"Java not found"

# Check Java installation
java -version

# Should show Java 8+ (64-bit)
# If not installed, install Java 8 or higher

"Voxel file not found"

Make sure run_voxelization_batch.R completed successfully and check the output path.

"No transmittance or attenuation variables"

Your voxelization settings may not have enabled these variables. Check the XML configuration.

DTM Issues

If you have problems with the DTM file:

• Check that it's in proper ESRI ASCII grid format
• Verify the file with head data/dtm/your_dtm.asc
• See scripts/DTM_README.md for creating valid DTM files

Documentation

• scripts/README.md - Detailed script documentation
• scripts/DTM_README.md - DTM generation guide
• AMAPVox documentation: https://amapvox.org
• R package help: ?AMAPVox::plantAreaDensity
• PAD estimators paper: https://doi.org/10.23708/1AJNMP

Next Steps

1. Visualize PAD: Use ggplot2 to plot vertical profiles
2. Compute LAI: Use AMAPVox::plantAreaIndex()
3. Export data: Use data.table::fwrite() for CSV export
4. 3D visualization: Use the AMAPVox 3D viewer

See scripts/example_pad_analysis.R for examples.

Requirements

• R (3.6+)
• R packages: AMAPVox, data.table
• For DTM generation: lidR, terra, future
• Java (8+, 64-bit) for AMAPVox

License

Scripts are provided as-is for use with AMAPVox. See AMAPVox license at https://amapvox.org