Work on rasters

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
+}