#!/usr/bin/env python
# coding: utf-8

# In this demo, we follow [Tyler Sloan](https://quorumetrix.blogspot.com/2018/06/visualizing-lidar-data-with-datashader.html)'s walkthrough where he downloads LiDAR data from the City of Montreal (featuring the 1976 Olympic Stadium) and visualizes it ([original notebook here](https://github.com/tsloan1377/montreal_open_data/blob/master/lidar_datashader_blog.ipynb)).
# 
# Make sure you've worked through the [Montreal data demo](Demo%20using%20Montreal%20LiDAR%20data.ipynb) first.

# ## Imports

# In[4]:


import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

import os
import imageio
from laspy.file import File
import datashader as ds
import datashader.transfer_functions as tf
from matplotlib import cm


# We already downloaded the tiles around the [site of Avebury](https://en.wikipedia.org/wiki/Avebury) from the [data.gov.uk](https://environment.data.gov.uk/ds/survey/index.jsp#/survey?grid=SU17) service. These tiles are in the [aveburytiles](/aveburytiles) folder. Unzip them, and explore! Try changing the various `class` parameters in the visualizations. Reuse code blocks

# In[23]:


get_ipython().system('LAStools/bin/laszip -i aveburtytiles/SU1477_P_8065_20120113_20120113.laz -o avebury1.las')


# ## Load LIDAR file with LasPy

# In[24]:


sample_data = 'avebury1.las'
export_path = 'export//'


# In[25]:


inFile = File(sample_data, mode='r')
df = pd.DataFrame() 
df['X'] = inFile.X 
df['Y'] = inFile.Y 
df['Z'] = inFile.Z
df['class'] = inFile.classification
display(df)


# ## Plot the tile using datashader

# ## Separate the different classes of pixel

# In[ ]:


cvs = ds.Canvas(plot_width=1000, plot_height=1000)
agg = cvs.points(df, 'X', 'Y', ds.mean('Z'))
img = tf.shade(agg)#, cmap=['lightblue', 'darkblue'], how='log')
tf.set_background(tf.shade(agg, cmap=cm.inferno),"black")


# In[ ]:


# Create a dataframe containing only the lidar voxels for buildings.
class_df = df.loc[df['class'] == 2]

# Visualize with datashader
cvs = ds.Canvas(plot_width=1000, plot_height=1000)
agg = cvs.points(class_df, 'X', 'Y', ds.mean('Z'))
img = tf.shade(agg)#, how='log')
tf.set_background(tf.shade(agg, cmap=cm.inferno),"black")


# In[ ]:


# Combine multiple classes of voxels that contain levels of vegetation
veg_df = df.loc[(df['class'] > 2) & (df['class'] < 6)]

# Visualize with datashader
cvs = ds.Canvas(plot_width=1000, plot_height=1000)
agg = cvs.points(veg_df, 'X', 'Y', ds.mean('Z'))
img = tf.shade(agg)#, how='log')
tf.set_background(tf.shade(agg, cmap=cm.inferno),"black")


# ## Create a 3D surface visualization
# 
# The code below is computationally intensive; it will take some time. The results will be written to a new folder called 'export'; you can open that folder by clicking on the jupyter logo at top and then clicking on the 'export' folder.
# 
# You will know the code is finished when the `[*]` at the left of the code block changes to a number.

# In[ ]:


# Use entire image containing olympic stadium, etc.
X = df['X']
Y = df['Y']
Z = df['Z']

# Downsample x and y
ds_factor = 500
ds_x = X[::ds_factor] 
ds_y = Y[::ds_factor] 
ds_z = Z[::ds_factor] 

##### Export the gif
frames = []
identifier = 'bigO_tile_downsample_' + str(ds_factor) + 'x_surface_lidar'  

if not os.path.exists(export_path):
    os.makedirs(export_path)

fig = plt.figure(figsize = (10,10)) 
ax = fig.add_subplot(111, projection='3d')

# Plot the surface.
surf = ax.plot_trisurf(ds_x, ds_y, ds_z, cmap=cm.inferno,
                       linewidth=0, antialiased=False)

for angle in range(0, 360):
    ax.view_init(60, angle) # Higher angle than usually used of 30.
    ax.set_axis_off()

    # Draw the figure
    fig.canvas.draw()

    # Convert to numpy array, and append to list
    np_fig = np.fromstring(fig.canvas.tostring_rgb(), dtype=np.uint8, sep='')
    np_fig = np_fig.reshape(fig.canvas.get_width_height()[::-1] + (3,))
    frames.append(np_fig)

imageio.mimsave(export_path + identifier + '.gif', frames)


# In[ ]: