lidar-coloration.Rmd

---
title: "Coloring lidar point clouds with RGB imagery in R"
author: "Max Joseph"
date: "June 30, 2016"
output:
  html_document:
    keep_md: true
---

```{r setup, echo = FALSE}
library(knitr)
library(rgl)
knit_hooks$set(webgl = hook_webgl)
```

# Acquire lidar and RGB data from NEON

We will use some sample data provided by the National Ecological Observatory Network Airborne Observation Platform (NEON AOP), available [here](http://www.neonscience.org/data-resources/get-data/airborne-data). 

```{r test-bash, results='hide'}
library(rLiDAR)
library(rgl)
library(raster)

# download data if necessary
neon_dir <- list.files(getwd(), pattern = 'NEON_AOP_sample_data')
if (length(neon_dir) == 0) {
  system("sample_d=AOP_sample_data.zip")
  system("sample_prefix=http://neonhighered.org/AOPSampleData/")
  system("wget $sample_prefix$sample_d")
  system("unzip $sample_d && rm $sample_d")
}
```

# Load lidar point cloud and visualize

```{r load-lidar}
gen_lidar_path <- function() {
  sample_dir <- list.files(getwd(), pattern = 'NEON_AOP_sample_data')
  lidar_dir <- 'LiDAR'
  type <- 'Discrete_LiDAR'
  which_lidar <- 'Point_Cloud'
  filename <- '2013_SJER_AOP_point_cloud_unclassified.las'
  file.path(sample_dir, lidar_dir, type, which_lidar, filename)
}

path_to_lidar <- gen_lidar_path()
lidar_d <- readLAS(path_to_lidar)
str(lidar_d)
```

# Load RGB imagery and visualize

```{r load-rgb}
gen_rgb_path <- function() {
  sample_dir <- list.files(getwd(), pattern = 'NEON_AOP_sample_data')
  rgb_dir <- 'RGB_Camera'
  filename <- '2013_SJER_AOP_Camera_sample.tif'
  file.path(sample_dir, rgb_dir, filename)
}

path_to_rgb <- gen_rgb_path()
rgb_image <- stack(path_to_rgb)
plotRGB(rgb_image)
```

# Combine lidar with RGB imagery

To link the lidar data to the RGB image, the coordinate systems must match (this requirement is satisfied in the NEON data).
Let's visualize the locations of our lidar points.

```{r testgl, webgl=TRUE}
rgb_vals <- extract(rgb_image, lidar_d[, c('X', 'Y')])
color_vector <- rgb(red = rgb_vals[, 1], 
                    green = rgb_vals[, 2], 
                    blue = rgb_vals[, 3], 
                    maxColorValue = 255)
plot3d(lidar_d[, 1:3], col = color_vector,
       xlab = 'Easting', ylab = 'Northing', zlab = 'Elevation')
```