ld2daps/Notebooks/GDAL vs Matplotlib.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GDAL vs Matplotlib (vs OpenCV vs LibTIFF)\n",
    "\n",
    "We have in data fusion contest dataset TIFF file with 32 bit per sample. For now TRISKELE only work with less than 17 bps: `BOOST_ASSERT (bits < 17);`. I want to ensure that we can oppen such data with Python.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import gdal\n",
    "from pathlib import Path\n",
    "import subprocess\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import cv2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = Path('../Data/phase1_rasters/DSM_C12/UH17c_GEF051_TR.tif')\n",
    "ofile = Path('../Res/python_out.tif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "info = subprocess.Popen(['tiffinfo', file], stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "print(info.communicate()[0].decode())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat_data = plt.imread(file)\n",
    "mat_data.shape, mat_data.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdl_data = gdal.Open(str(file))\n",
    "gdl_data.GetMetadata(), gdl_data.RasterCount"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdl_data.GetRasterBand(1).ReadAsArray().shape, gdl_data.GetRasterBand(1).ReadAsArray().dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## There.\n",
    "\n",
    "`matplotlib` is derping around with bps. \n",
    "\n",
    "Maybe each byte is split `[a, b, c, d]` as $V = a 2^{24} + b 2^{16} + c 2^8 + d$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(np.power(2, (np.arange(4)[::-1] * 8)) * mat_data).sum(axis=2).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster = gdl_data.GetRasterBand(1).ReadAsArray()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(32,9))\n",
    "plt.imshow(raster)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(raster.reshape(-1), 100, log=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets\n",
    "\n",
    "def disp(x):\n",
    "    display(x)\n",
    "    plt.figure(figsize=(32,9))\n",
    "    plt.imshow(raster * (raster < x))\n",
    "    plt.colorbar()\n",
    "    plt.show()\n",
    "    plt.hist(raster.reshape(-1) * (raster.reshape(-1) < x), 100, log=True)\n",
    "    plt.show()\n",
    "    \n",
    "interact(disp, x=widgets.IntSlider(min=raster.min(),max=80,value=80));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_df = pd.DataFrame(raster.reshape(-1))\n",
    "raster_df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rasterf = raster.copy()\n",
    "rasterf[rasterf == rasterf.max()] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rasterf_df = pd.DataFrame(rasterf.reshape(-1))\n",
    "rasterf_df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(rasterf)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write TIFF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### OpenCV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_cv2 = cv2.imread(str(file), -1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(raster_cv2 != raster).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "str(ofile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls ../Res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rasterf_im = np.expand_dims(rasterf, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv2.imwrite(str(ofile), rasterf_im)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GDAL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdl_data.RasterXSize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = gdl_data.RasterXSize\n",
    "rows = gdl_data.RasterYSize\n",
    "\n",
    "format = 'GTiff'\n",
    "driver = gdal.GetDriverByName(format)\n",
    "dst_ds = driver.Create(str(ofile), cols, rows, 1, gdal.GDT_Byte)\n",
    "dst_ds.GetRasterBand(1).WriteArray(rasterf, 0, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dst_ds = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Test "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tiffinfo ../Res/python_out.tif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(rasterf)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_ofile = cv2.imread(str(ofile))\n",
    "display(test_ofile.shape)\n",
    "todf = pd.DataFrame(test_ofile.reshape(-1,3))\n",
    "display(todf.describe())\n",
    "plt.imshow(test_ofile[:,:,0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LibTiff\n",
    "\n",
    "### Write"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import libtiff\n",
    "\n",
    "tiff = libtiff.TIFFimage(rasterf.astype(np.float16), description='test')\n",
    "tiff.write_file(ofile)#, compression='lzw')\n",
    "del tiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tiffinfo ../Res/python_out.tif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ogdl_data = gdal.Open(str(ofile))\n",
    "ogdl_data.GetMetadata(), ogdl_data.RasterCount"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "otest = ogdl_data.GetRasterBand(1).ReadAsArray()\n",
    "otest.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(otest)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(otest != rasterf).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Précision au millimètre :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.isclose(otest, rasterf, rtol=1e-3).sum() / rasterf.size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tiff_r = libtiff.TIFF.open(file, mode='r')\n",
    "tiff_r.read_image().shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(tiff_r.read_image())\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(tiff_r.read_image() != raster).sum()"
   ]
  }
 ],
 "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
}