{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "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)).\n", "\n", "Make sure you've worked through the [Montreal data demo](Demo%20using%20Montreal%20LiDAR%20data.ipynb) first." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Imports" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from mpl_toolkits.mplot3d import Axes3D\n", "\n", "import os\n", "import imageio\n", "from laspy.file import File\n", "import datashader as ds\n", "import datashader.transfer_functions as tf\n", "from matplotlib import cm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "!LAStools/bin/laszip -i aveburtytiles/SU1477_P_8065_20120113_20120113.laz -o avebury1.las" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load LIDAR file with LasPy" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "sample_data = 'avebury1.las'\n", "export_path = 'export//'" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
XYZclass
04146878117785233175112
14146871617785165175173
24146873417785159175223
34146879317785220175121
44146904417785222175052
54146898017785156175153
64146891217785083175072
74146884717785016175141
84146885317784998175112
94146891417785062175112
104146897717785128175143
114146903817785192175143
124146931217785213175121
134146924417785142175082
144146917817785073175082
154146911317785004175112
164146904617784935175128
174146898017784866175132
184146898417784849175152
194146904517784912175122
204146910717784976175112
214146916717785039175082
224146923017785104175111
234146929117785169175098
244146935417785235175143
254146958617785216175042
264146951917785144174998
274146945617785080175133
284146938817785008175072
294146932317784941175141
...............
8500324149985717755556181382
8500334149991417755615181382
8500344149997217755674181352
8500354150000017755624181372
8500364149994217755565181372
8500374149988417755507181412
8500384149982717755447181402
8500394149976917755388181422
8500404149971017755328181493
8500414149965317755270181462
8500424149959517755211181501
8500434149953817755152181482
8500444149948017755093181502
8500454149942317755035181522
8500464149963217755061181502
8500474149968917755119181538
8500484149974817755178181482
8500494149980517755237181452
8500504149986217755295181471
8500514149992117755355181412
8500524149998017755413181392
8500534149999217755346181402
8500544149993417755287181442
8500554149987517755228181483
8500564149981917755170181428
8500574149976117755110181492
8500584149970217755051181548
8500594149990617755054181518
8500604149996517755114181448
8500614149998717755065181472
\n", "

850062 rows × 4 columns

\n", "
" ], "text/plain": [ " X Y Z class\n", "0 41468781 17785233 17511 2\n", "1 41468716 17785165 17517 3\n", "2 41468734 17785159 17522 3\n", "3 41468793 17785220 17512 1\n", "4 41469044 17785222 17505 2\n", "5 41468980 17785156 17515 3\n", "6 41468912 17785083 17507 2\n", "7 41468847 17785016 17514 1\n", "8 41468853 17784998 17511 2\n", "9 41468914 17785062 17511 2\n", "10 41468977 17785128 17514 3\n", "11 41469038 17785192 17514 3\n", "12 41469312 17785213 17512 1\n", "13 41469244 17785142 17508 2\n", "14 41469178 17785073 17508 2\n", "15 41469113 17785004 17511 2\n", "16 41469046 17784935 17512 8\n", "17 41468980 17784866 17513 2\n", "18 41468984 17784849 17515 2\n", "19 41469045 17784912 17512 2\n", "20 41469107 17784976 17511 2\n", "21 41469167 17785039 17508 2\n", "22 41469230 17785104 17511 1\n", "23 41469291 17785169 17509 8\n", "24 41469354 17785235 17514 3\n", "25 41469586 17785216 17504 2\n", "26 41469519 17785144 17499 8\n", "27 41469456 17785080 17513 3\n", "28 41469388 17785008 17507 2\n", "29 41469323 17784941 17514 1\n", "... ... ... ... ...\n", "850032 41499857 17755556 18138 2\n", "850033 41499914 17755615 18138 2\n", "850034 41499972 17755674 18135 2\n", "850035 41500000 17755624 18137 2\n", "850036 41499942 17755565 18137 2\n", "850037 41499884 17755507 18141 2\n", "850038 41499827 17755447 18140 2\n", "850039 41499769 17755388 18142 2\n", "850040 41499710 17755328 18149 3\n", "850041 41499653 17755270 18146 2\n", "850042 41499595 17755211 18150 1\n", "850043 41499538 17755152 18148 2\n", "850044 41499480 17755093 18150 2\n", "850045 41499423 17755035 18152 2\n", "850046 41499632 17755061 18150 2\n", "850047 41499689 17755119 18153 8\n", "850048 41499748 17755178 18148 2\n", "850049 41499805 17755237 18145 2\n", "850050 41499862 17755295 18147 1\n", "850051 41499921 17755355 18141 2\n", "850052 41499980 17755413 18139 2\n", "850053 41499992 17755346 18140 2\n", "850054 41499934 17755287 18144 2\n", "850055 41499875 17755228 18148 3\n", "850056 41499819 17755170 18142 8\n", "850057 41499761 17755110 18149 2\n", "850058 41499702 17755051 18154 8\n", "850059 41499906 17755054 18151 8\n", "850060 41499965 17755114 18144 8\n", "850061 41499987 17755065 18147 2\n", "\n", "[850062 rows x 4 columns]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "inFile = File(sample_data, mode='r')\n", "df = pd.DataFrame() \n", "df['X'] = inFile.X \n", "df['Y'] = inFile.Y \n", "df['Z'] = inFile.Z\n", "df['class'] = inFile.classification\n", "display(df)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot the tile using datashader" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Separate the different classes of pixel" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cvs = ds.Canvas(plot_width=1000, plot_height=1000)\n", "agg = cvs.points(df, 'X', 'Y', ds.mean('Z'))\n", "img = tf.shade(agg)#, cmap=['lightblue', 'darkblue'], how='log')\n", "tf.set_background(tf.shade(agg, cmap=cm.inferno),\"black\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a dataframe containing only the lidar voxels for buildings.\n", "class_df = df.loc[df['class'] == 2]\n", "\n", "# Visualize with datashader\n", "cvs = ds.Canvas(plot_width=1000, plot_height=1000)\n", "agg = cvs.points(class_df, 'X', 'Y', ds.mean('Z'))\n", "img = tf.shade(agg)#, how='log')\n", "tf.set_background(tf.shade(agg, cmap=cm.inferno),\"black\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Combine multiple classes of voxels that contain levels of vegetation\n", "veg_df = df.loc[(df['class'] > 2) & (df['class'] < 6)]\n", "\n", "# Visualize with datashader\n", "cvs = ds.Canvas(plot_width=1000, plot_height=1000)\n", "agg = cvs.points(veg_df, 'X', 'Y', ds.mean('Z'))\n", "img = tf.shade(agg)#, how='log')\n", "tf.set_background(tf.shade(agg, cmap=cm.inferno),\"black\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create a 3D surface visualization\n", "\n", "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.\n", "\n", "You will know the code is finished when the `[*]` at the left of the code block changes to a number." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/srv/conda/lib/python3.6/site-packages/ipykernel_launcher.py:34: DeprecationWarning: The binary mode of fromstring is deprecated, as it behaves surprisingly on unicode inputs. Use frombuffer instead\n" ] } ], "source": [ "# Use entire image containing olympic stadium, etc.\n", "X = df['X']\n", "Y = df['Y']\n", "Z = df['Z']\n", "\n", "# Downsample x and y\n", "ds_factor = 500\n", "ds_x = X[::ds_factor] \n", "ds_y = Y[::ds_factor] \n", "ds_z = Z[::ds_factor] \n", "\n", "##### Export the gif\n", "frames = []\n", "identifier = 'bigO_tile_downsample_' + str(ds_factor) + 'x_surface_lidar' \n", "\n", "if not os.path.exists(export_path):\n", " os.makedirs(export_path)\n", "\n", "fig = plt.figure(figsize = (10,10)) \n", "ax = fig.add_subplot(111, projection='3d')\n", "\n", "# Plot the surface.\n", "surf = ax.plot_trisurf(ds_x, ds_y, ds_z, cmap=cm.inferno,\n", " linewidth=0, antialiased=False)\n", "\n", "for angle in range(0, 360):\n", " ax.view_init(60, angle) # Higher angle than usually used of 30.\n", " ax.set_axis_off()\n", "\n", " # Draw the figure\n", " fig.canvas.draw()\n", "\n", " # Convert to numpy array, and append to list\n", " np_fig = np.fromstring(fig.canvas.tostring_rgb(), dtype=np.uint8, sep='')\n", " np_fig = np_fig.reshape(fig.canvas.get_width_height()[::-1] + (3,))\n", " frames.append(np_fig)\n", "\n", "imageio.mimsave(export_path + identifier + '.gif', frames)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }