{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Custom Raster Gen for LD2DAPs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sys.path.append(\"..\")\n",
    "import rasterizer\n",
    "import raster_assistant as ra\n",
    "\n",
    "sys.path.append('../triskele/python/')\n",
    "import triskele\n",
    "\n",
    "figsize = np.array((16, 3)) * 1.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load LAS data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C1 = ra.bulk_load('../Data/lidar/C1', 'C1', filter_treshold=.01, clip_treshold=.1, dtype=np.float32)\n",
    "C2 = ra.bulk_load('../Data/lidar/C2', 'C2', filter_treshold=.01, clip_treshold=.1, dtype=np.float32)\n",
    "C3 = ra.bulk_load('../Data/lidar/C3', 'C3', filter_treshold=.2, clip_treshold=.1, dtype=np.float32)\n",
    "C12 = ra.bulk_load(['../Data/lidar/C1', '../Data/lidar/C2'], 'C12', filter_treshold=.01, clip_treshold=.1, dtype=np.float32)\n",
    "C123 = ra.bulk_load('../Data/lidar', 'C123', filter_treshold=.08, clip_treshold=.1, dtype=np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster = ra.rasterize_cache('num_returns', C1, .5, 'nearest', False, cache_dir='../Res/enrichment_rasters')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster, origin='upper')\n",
    "plt.show()\n",
    "plt.imsave('../Res/raster_validation.png', raster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Raster Validation\n",
    "\n",
    "### Rasterize some data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster = rasterizer.rasterize(C1.spatial, C1.intensity, 0.5, dtype=np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster, origin='upper')\n",
    "plt.show()\n",
    "plt.imsave('../Res/raster_validation.png', raster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(raster.flatten(), bins=1000)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write TIFF file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "triskele.write('../Res/validation.tiff', raster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare with DFC dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfc = triskele.read('../Data/phase1_rasters/Intensity_C1/UH17_GI1F051_TR.tif')\n",
    "our = triskele.read('../Res/validation.tiff')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Filter DFC with same parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ra.extremum_clip(dfc, treshold=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Display Stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfc.shape, our.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfc.dtype, our.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(dfc.flatten(), bins=1000, label='DFC')\n",
    "plt.hist(our.flatten(), bins=1000, label='Our', alpha=.8)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Display Rasters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axs = plt.subplots(2, figsize=figsize * 2)\n",
    "\n",
    "axs[0].imshow(dfc)\n",
    "axs[0].set_title('DFC')\n",
    "axs[1].imshow(our)\n",
    "axs[1].set_title('Our')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Raster Pool Process"
   ]
  },
  {
   "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",
    "def rasterize_cache_mp(data_var, field, res, method, reverse, cache):\n",
    "    print('FROM WORKER: Start processing {} {} with {} resolution and {} interpolation {}'.format(data_var, field, res, method, 'resversed' if reverse else ''))\n",
    "    if data_var == 'C1':\n",
    "        ra.rasterize_cache(field, C1, res, method, reverse, cache)    \n",
    "    if data_var == 'C2':\n",
    "        ra.rasterize_cache(field, C2, res, method, reverse, cache)\n",
    "    if data_var == 'C3':\n",
    "        ra.rasterize_cache(field, C3, res, method, reverse, cache)\n",
    "    if data_var == 'C12':\n",
    "        ra.rasterize_cache(field, C12, res, method, reverse, cache)\n",
    "    if data_var == 'C123':\n",
    "        ra.rasterize_cache(field, C123, res, method, reverse, cache)\n",
    "    \n",
    "pool = Pool(processes=9)\n",
    "\n",
    "job_args = list()\n",
    "\n",
    "for res in (0.5, 1., 2., 3., 5., 10., .1):\n",
    "    for reverse in (False, True):\n",
    "        for inter in ('linear', 'nearest', 'cubic-clip'):\n",
    "            for field in ('z', 'intensity', 'num_returns'):\n",
    "                for data in ('C1', 'C2', 'C3', 'C12', 'C123'):\n",
    "                    job_args.append([data, field, res, inter, reverse, '../Res/enrichment_rasters'])\n",
    "                    \n",
    "print(\"Job list length: {}\".format(len(job_args)))\n",
    "                    \n",
    "for i in pool.starmap(rasterize_cache_mp, job_args):\n",
    "    pass"
   ]
  },
  {
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}