#!/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 visualize it ([original notebook here](https://github.com/tsloan1377/montreal_open_data/blob/master/lidar_datashader_blog.ipynb)).
# 
# This notebook depends on the LAStools for working with LiDAR - we only really need the LASzip tool, but we've installed the entire suite. You can see what other tools are available from the filemanager in the [LAStools/bin](Lastools/bin) directory.

# ## Imports

# In[12]:


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 use `curl` to get the `.laz`, the compressed lidar data from the city of Montreal which we save as mtl.lz

# In[1]:


get_ipython().system('curl http://depot.ville.montreal.qc.ca/geomatique/lidar_aerien/2015/300-5046_2015.laz > mtl.laz')
    


# Now we need to make the LASzip file executable; you'll need to do this to the other tools individually if you wish to use any of the other tools in that bin, by the way.

# In[6]:


get_ipython().system('chmod 755 LAStools/bin/laszip')


# We unzip the lidar data, specifying the input and outputs:

# In[9]:


get_ipython().system('LAStools/bin/laszip -i mtl.laz -o mtl.las')


# ## Load LIDAR file with LasPy

# In[10]:


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


# In[13]:


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

# In[14]:


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")


# ## Separate the different classes of pixel

# In[15]:


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

# 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[16]:


# 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

# 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[ ]: