diff --git a/Notebooks/Custom rasters from LiDAR data.ipynb b/Notebooks/Custom rasters from LiDAR data.ipynb new file mode 100644 index 0000000..80104e7 --- /dev/null +++ b/Notebooks/Custom rasters from LiDAR data.ipynb @@ -0,0 +1,212 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import laspy\n", + "\n", + "figsize = (16, 9)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Load LAS files data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def find_las(path):\n", + " print(path)\n", + " las = list()\n", + " \n", + " if path.is_dir():\n", + " for child in path.iterdir():\n", + " las.extend(find_las(child))\n", + " \n", + " if path.is_file() and path.suffix == '.las':\n", + " las.append(path)\n", + " \n", + " return las" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "files = ['../Data/lidar/C1/272056_3289689.las',\n", + " '../Data/lidar/C1/272652_3289689.las',\n", + " '../Data/lidar/C1/273248_3289689.las',\n", + " '../Data/lidar/C1/273844_3289689.las']\n", + "\n", + "folder = Path('../Data/lidar/')\n", + "files = find_las(folder)\n", + "\n", + "xdata = list()\n", + "ydata = list()\n", + "idata = list()\n", + "\n", + "for file in files:\n", + " infile = laspy.file.File(file)\n", + " print('file: {}, point count: {}'.format(file, infile.header.get_count()))\n", + " xdata.append(infile.x)\n", + " ydata.append(infile.y)\n", + " idata.append(infile.intensity)\n", + " \n", + "data = np.array((np.concatenate(xdata), np.concatenate(ydata), np.concatenate(idata))).T\n", + "data.shape" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Display a 2D hist" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "resolution = 2\n", + "ground_size = np.array((data[:,0].max() - data[:,0].min(), data[:,1].max() - data[:,1].min()))\n", + "display(ground_size)\n", + "\n", + "bins = ground_size / resolution" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=figsize)\n", + "hist, xedges, yedges = np.histogram2d(data[:,0], data[:,1], bins=bins)\n", + "plt.imshow(hist.T, origin='lower', extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.imsave('../Res/hist.png', hist.T, origin='lower')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Display intensity map" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=figsize)\n", + "hist_i, xedges, yedges = np.histogram2d(data[:,0], data[:,1], bins=bins, weights=data[:,2])\n", + "plt.imshow(hist_i.T / hist.T, origin='lower', extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It seems that we have no data at this resolution. Should interpolate ?\n", + "\n", + "### Better intensity display" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "i = hist_i.T / hist.T\n", + "f = np.logical_not(np.isnan(i))\n", + "plt.hist(i[f], bins=1000)\n", + "plt.show()\n", + "\n", + "i[f].min(), i[f].max()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "imax_disp = 10000\n", + "\n", + "idisp = i\n", + "idisp[i > imax_disp] = imax_disp\n", + "\n", + "plt.figure(figsize=figsize)\n", + "plt.imshow(idisp, origin='lower', extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.imsave('../Res/hist_int.png', idisp, origin='lower')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.scatter(infile.x, infile.y)\n", + "plt.show()" + ] + } + ], + "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" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}