cblue_diagnostics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from pathlib import Path\n",
    "import geopandas as gpd\n",
    "from IPython.core.display import HTML, display\n",
    "\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src2016 = rasterio.open(r'z:\\St.JosephBay2016Composite.tif')\n",
    "src2017 = rasterio.open(r'z:\\St.JosephBay2017Composite.tif')\n",
    "\n",
    "array2016 = src2016.read(1)\n",
    "array2017 = src2017.read(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(src.transform)\n",
    "print(src.crs)\n",
    "print(src.crs.wkt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src.transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rasterio.plot import show_hist\n",
    "show_hist(src, bins=50, lw=0.0, stacked=False, alpha=0.3, title=\"Histogram\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax2016, ax2017) = plt.subplots(1, 2, figsize=(21,14))\n",
    "show((src2016), ax=ax2016)\n",
    "show((src2017), ax=ax2017)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyplot.imshow(array2017, cmap='viridis')\n",
    "pyplot.show() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(array, cmap='jet')\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 = Path(r'C:\\QAQC_contract\\marco_island\\tpu_output\\2016_429500e_2871000n.parquet')\n",
    "p2 = Path(r'C:\\QAQC_contract\\marco_island\\tpu_output\\2016_429500e_2871500n.parquet')\n",
    "\n",
    "df = pd.read_parquet(p1)\n",
    "df1 = df[['las_x', 'las_y', 'total_tvu']]\n",
    "\n",
    "df = pd.read_parquet(p2)\n",
    "df2 = df[['las_x', 'las_y', 'total_tvu']]\n",
    "\n",
    "df_tvu = pd.concat([df1, df2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr = df.corr()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(15, 15))\n",
    "im = ax.imshow(corr, cmap='seismic')\n",
    "\n",
    "fig.colorbar(im, ax=ax)\n",
    "\n",
    "# We want to show all ticks...\n",
    "ax.set_xticks(np.arange(len(corr.columns)))\n",
    "ax.set_yticks(np.arange(len(corr.columns)))\n",
    "\n",
    "# ... and label them with the respective list entries\n",
    "ax.set_xticklabels(corr.columns)\n",
    "ax.set_yticklabels(corr.columns)\n",
    "\n",
    "# Rotate the tick labels and set their alignment.\n",
    "pyplot.setp(ax.get_xticklabels(), rotation=45, ha=\"right\", rotation_mode=\"anchor\")\n",
    "\n",
    "ax.set_title(\"Harvest of local farmers (in tons/year)\")\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import holoviews as hv\n",
    "import dask.dataframe as dd\n",
    "import datashader as ds, datashader.geo\n",
    "\n",
    "from holoviews import opts\n",
    "from holoviews.operation.datashader import datashade, rasterize\n",
    "\n",
    "hv.extension('bokeh')\n",
    "\n",
    "points = hv.Points(df_tvu, ['las_x', 'las_y'])\n",
    "\n",
    "points1 = hv.Points(df1, ['las_x', 'las_y'])\n",
    "points2 = hv.Points(df2, ['las_x', 'las_y'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.cm import viridis, jet\n",
    "from bokeh.palettes import GnBu9, Turbo\n",
    "\n",
    "width = 600\n",
    "height = 900\n",
    "\n",
    "las_pts = datashade(points).opts(width=width, height=height, bgcolor=\"lightgray\")\n",
    "\n",
    "opts = dict(colorbar=True, width=width, height=height, cmap=Turbo[256], logz=True)\n",
    "tpu = rasterize(points, aggregator=ds.mean('total_tvu')).opts(**opts)\n",
    "\n",
    "las_pts + tpu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from holoviews.operation import histogram\n",
    "\n",
    "tvu_min = 0.05\n",
    "tvu_max = 0.13\n",
    "tvu_range = (tvu_min, tvu_max)\n",
    "\n",
    "hist_height = 200 \n",
    "hist_width = 800\n",
    "hist_alpha = 0.3\n",
    "num_bins = 200\n",
    "\n",
    "selected1 = points1.select(total_tvu=tvu_range)\n",
    "selected2 = points2.select(total_tvu=tvu_range)\n",
    "\n",
    "opts = dict(width=hist_width, height=hist_height, alpha=hist_alpha, line_color=None)\n",
    "\n",
    "hist1 = histogram(selected1, bin_range=tvu_range, bins=num_bins).opts(**opts)\n",
    "hist2 = histogram(selected2, bin_range=tvu_range, bins=num_bins).opts(**opts)\n",
    "\n",
    "hist1 * hist2\n",
    "#dummy"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}