{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "from scipy.interpolate import griddata\n",
    "import matplotlib.pyplot as plt\n",
    "import laspy\n",
    "\n",
    "triskele_path = Path('../triskele/python/')\n",
    "sys.path.append(str(triskele_path.resolve()))\n",
    "import triskele\n",
    "\n",
    "figsize = np.array((16, 3)) * 1.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load LiDAR data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_las(path):\n",
    "    path = Path(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": [
    "def auto_filter(data, treshold=.5):\n",
    "    tresholds = np.percentile(data, [treshold, 100 - treshold])\n",
    "\n",
    "    data[data < tresholds[0]] = tresholds[0]\n",
    "    data[data > tresholds[1]] = tresholds[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bulk_load(path, name=None):\n",
    "    data = {'file': path, 'name': name}\n",
    "\n",
    "    attributes = ['x', 'y', 'z', 'intensity', 'num_returns', 'scan_angle_rank', 'pt_src_id', 'gps_time']\n",
    "    for a in attributes:\n",
    "        data[a] = list()\n",
    "        \n",
    "    print('Load data...')\n",
    "    for f in find_las(path):\n",
    "        print('{}: '.format(f), end='')\n",
    "        infile = laspy.file.File(f)\n",
    "        for i, a in enumerate(attributes):\n",
    "            print('\\r {}: [{:3d}%]'.format(f, int(i/len(attributes) * 100)), end='')\n",
    "            data[a].extend(getattr(infile, a))\n",
    "        infile.close()\n",
    "        print('\\r {}: [Done]'.format(f))\n",
    "\n",
    "    \n",
    "    print('Create matrices...', end='')\n",
    "    for i, a in enumerate(attributes):\n",
    "        print('\\rCreate matrices: [{:3d}%]'.format(int(i/len(attributes) * 100)), end='')\n",
    "        data[a] = np.array(data[a])\n",
    "    print('\\rCreate matrices: [Done]')\n",
    "        \n",
    "    print('Filter data...', end='')\n",
    "    for i, a in enumerate(attributes):\n",
    "        print('\\rFilter data: [{:3d}%]'.format(int(i/len(attributes) * 100)), end='')\n",
    "        if a == 'x' or a == 'y':\n",
    "            continue\n",
    "        auto_filter(data[a])\n",
    "    print('\\rFilter data: [Done]')\n",
    "\n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C1 = bulk_load('../Data/lidar/C1/', 'C1')\n",
    "C2 = bulk_load('../Data/lidar/C2/', 'C2')\n",
    "C3 = bulk_load('../Data/lidar/C3/', 'C3')\n",
    "C123 = bulk_load('../Data/lidar/', 'C123') # Todo: Copy C1 C2 C3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define raster function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rasterize(coords, values, resolution=1, method='nearest'):\n",
    "    if isinstance(coords, np.ndarray):\n",
    "        xmin, xmax = coords[:,0].min(), coords[:,0].max()\n",
    "        ymin, ymax = coords[:,1].min(), coords[:,1].max()\n",
    "    else:\n",
    "        xmin, xmax = coords[0].min(), coords[0].max()\n",
    "        ymin, ymax = coords[1].min(), coords[1].max()\n",
    "\n",
    "    gridx, gridy = np.mgrid[xmin:xmax:resolution, ymin:ymax:resolution]\n",
    "    \n",
    "    return griddata(coords, values, (gridx, gridy), method=method).T # linear, nearest, cubic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def raster_show_and_save(coords, values, resolution=1, method='nearest', name=None, save_dir='../Res'):\n",
    "    \"\"\"Display and save rasters. If name is provided, will try to cache result to reload it later.\"\"\"\n",
    "    \n",
    "    if name is not None:\n",
    "        base_name = '{}/{}'.format(save_dir, name)\n",
    "        tiff_file = Path(base_name + '.tiff')\n",
    "        png_file = Path(base_name + '.png')\n",
    "        \n",
    "    if name is not None and tiff_file.exists() :\n",
    "        print ('WARNING: Loading cached result {}'.format(tiff_file))\n",
    "        raster = triskele.read(tiff_file)\n",
    "    else:         \n",
    "        raster = rasterize(coords, values, resolution, method)\n",
    "    \n",
    "    plt.figure(figsize=figsize)\n",
    "    plt.imshow(raster, origin='lower')\n",
    "    plt.colorbar()\n",
    "    plt.title(name)\n",
    "    plt.show()\n",
    "    \n",
    "    if name is not None:\n",
    "        plt.imsave(png_file, raster, origin='lower')\n",
    "        triskele.write(tiff_file, raster)\n",
    "        \n",
    "    return raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def raster_assistant(data, field='z', resolution='1', method='nearest', weight_field=None):\n",
    "    if weight_field is None:\n",
    "        name = '{}_{}_{}_{}'.format(data['name'], field, str(resolution).replace('.', '_'), method)\n",
    "        raster_show_and_save((data['x'], data['y']), data[field], resolution, method, name)\n",
    "    else:\n",
    "        name = '{}_{}_{}_{}_w{}'.format(data['name'], field, str(resolution).replace('.', '_'), method, weight_field)\n",
    "        raster_show_and_save((data['x'], data['y']), data[field] * data[weight_field], resolution, method, name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def raster_assistant_mp(data, field='z', resolution='1', method='nearest', weight_field=None):\n",
    "    if data == 'C1':\n",
    "        raster_assistant(C1, field, resolution, method)\n",
    "    if data == 'C2':\n",
    "        raster_assistant(C2, field, resolution, method)\n",
    "    if data == 'C3':\n",
    "        raster_assistant(C3, field, resolution, method)\n",
    "    if data == 'C123':\n",
    "        raster_assistant(C123, field, resolution, method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rl = raster_assistant(C1, resolution=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpolation Strategies for Global Griddata Rasterisation\n",
    "\n",
    "## Fine textured (visible) rasters\n",
    "\n",
    "### Nearest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "null = raster_assistant(C1, 'intensity', 0.5, 'nearest')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No interpolation, accurate values. On the downside, there is a lot of aliasing. Maybe that's because take all the points without selection on the z (altitude) axis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "null = raster_assistant(C1, 'intensity', 0.5, 'linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cubic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are pretty bad, same remark than before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "null = raster_assistant(C1, 'intensity', 0.5, 'cubic')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cubic create new values. Still, seems pretty useless with the aliasing. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flat area (DSM) rasters\n",
    "\n",
    "### Nearest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "null = raster_assistant(C1, 'z', 0.5, 'nearest')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No interpolation, accurate values. On the downside, there is a lot of aliasing. Maybe that's because take all the points without selection on the z (altitude) axis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "null = raster_assistant(C1, 'z', 0.5, 'linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cubic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are pretty bad, same remark than before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "null = raster_assistant(C1, 'z', 0.5, 'cubic')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Weight with altitude\n",
    "\n",
    "## Visible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "null = raster_assistant(C1, 'intensity', 0.5, 'nearest', 'z')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "null = raster_assistant(C1, 'intensity', 0.5, 'linear', 'z')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DSM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "null = raster_assistant(C1, 'z', 0.5, 'nearest', 'z')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "null = raster_assistant(C1, 'z', 0.5, 'linear', 'z')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Other"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from multiprocessing import Pool, Process, Queue\n",
    "import multiprocessing as mp\n",
    "#mp.set_start_method('spawn')\n",
    "\n",
    "pool = Pool(processes=20)\n",
    "\n",
    "job_args = list()\n",
    "for field in ('z', 'intensity', 'num_returns', 'gps_time'):\n",
    "    for data in ('C1', 'C2', 'C3', 'C123'):\n",
    "        job_args.append([data, field, 1, 'nearest'])\n",
    "        \n",
    "for i in pool.starmap(raster_assistant_mp, job_args):\n",
    "    pass        "
   ]
  }
 ],
 "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
}