From d7f1163c06615e56a24f054bcbd2368b6bca7195 Mon Sep 17 00:00:00 2001 From: Karamaz0V1 Date: Fri, 31 Aug 2018 15:55:52 +0200 Subject: [PATCH] Try some raster combination --- Notebooks/Raster Factory.ipynb | 6 +- Notebooks/Raster combination.ipynb | 246 +++++++++++++++++++++++++++++ 2 files changed, 250 insertions(+), 2 deletions(-) create mode 100644 Notebooks/Raster combination.ipynb diff --git a/Notebooks/Raster Factory.ipynb b/Notebooks/Raster Factory.ipynb index 9587aec..46b8bca 100644 --- a/Notebooks/Raster Factory.ipynb +++ b/Notebooks/Raster Factory.ipynb @@ -43,7 +43,7 @@ "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', 'C1', filter_treshold=.01, 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)" ] }, @@ -230,6 +230,8 @@ " 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", @@ -241,7 +243,7 @@ " 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', 'C123'):\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", diff --git a/Notebooks/Raster combination.ipynb b/Notebooks/Raster combination.ipynb new file mode 100644 index 0000000..fdc081b --- /dev/null +++ b/Notebooks/Raster combination.ipynb @@ -0,0 +1,246 @@ +{ + "cells": [ + { + "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": [ + "# Raster Combination\n", + "\n", + "## DSM - DTM " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "C12 = ra.bulk_load(['../Data/lidar/C1', '../Data/lidar/C2'], 'C12', filter_treshold=.01, clip_treshold=.1, dtype=np.float32)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dsm = ra.rasterize_cache('z', C12, .5, 'nearest', False, '../Res/enrichment_rasters/')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dtm = ra.rasterize_cache('z', C12, .5, 'nearest', True, '../Res/enrichment_rasters/')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=figsize)\n", + "plt.imshow(dsm)\n", + "plt.colorbar()\n", + "plt.title('DSM')\n", + "plt.show()\n", + "\n", + "plt.figure(figsize=figsize)\n", + "plt.imshow(dtm)\n", + "plt.colorbar()\n", + "plt.title('DTM')\n", + "plt.show()\n", + "\n", + "plt.figure(figsize=figsize)\n", + "plt.imshow(dsm - dtm)\n", + "plt.colorbar()\n", + "plt.title('DSM - DTM')\n", + "plt.show()\n", + "\n", + "plt.imsave('../Res/dsm-dtm.png', dsm - dtm)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## NDVI\n", + "\n", + "With \n", + "\n", + "$NDVI = \\frac{NIR - Red}{NIR + Red}$\n", + "\n", + "and the wavelenth of the Titan\n", + "\n", + "| Lazer | Wavelenght | Color |\n", + "| --- | ---: | --- |\n", + "| C1 | 1550 nm | IR? |\n", + "| C2 | 1064 nm | NIR |\n", + "| C3 | 532 nm | Green |\n", + "\n", + "we can compute a NDVI like intensity raster with\n", + "\n", + "$NDVI_{like} = \\frac{C1 - C2}{C1 + C2}$" + ] + }, + { + "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)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "C1_raster = ra.rasterize_cache('intensity', C1, .5, 'linear', False, '../Res/enrichment_rasters/')\n", + "C2_raster = ra.rasterize_cache('intensity', C2, .5, 'linear', False, '../Res/enrichment_rasters/')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ndvi = (C1_raster - C2_raster) / (C1_raster + C2_raster)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=figsize)\n", + "plt.imshow(C1_raster)\n", + "plt.title('C1')\n", + "plt.colorbar()\n", + "plt.show()\n", + "\n", + "plt.figure(figsize=figsize)\n", + "plt.imshow(C2_raster)\n", + "plt.title('C2')\n", + "plt.colorbar()\n", + "plt.show()\n", + "\n", + "plt.figure(figsize=figsize)\n", + "plt.imshow(ndvi)\n", + "plt.title('NDVI')\n", + "plt.colorbar()\n", + "plt.show()\n", + "\n", + "plt.imsave('../Res/ndvi_linear.png', ndvi)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Normalized NDVI" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "C1_raster /= C1_raster.max()\n", + "C2_raster /= C2_raster.max()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ndvi = (C1_raster - C2_raster) / (C1_raster + C2_raster)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=figsize)\n", + "plt.imshow(C1_raster)\n", + "plt.title('C1')\n", + "plt.colorbar()\n", + "plt.show()\n", + "\n", + "plt.figure(figsize=figsize)\n", + "plt.imshow(C2_raster)\n", + "plt.title('C2')\n", + "plt.colorbar()\n", + "plt.show()\n", + "\n", + "plt.figure(figsize=figsize)\n", + "plt.imshow(ndvi)\n", + "plt.title('NDVI')\n", + "plt.colorbar()\n", + "plt.show()\n", + "\n", + "plt.imsave('../Res/ndvi_normalized_linear.png', ndvi)" + ] + }, + { + "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 +}