Ortho Image Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rasterio.plot import show\n",
    "import os\n",
    "import rasterio as rio\n",
    "import cv2\n",
    "import math\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "mpl.rcParams['figure.dpi']= 100\n",
    "plt.rcParams.update({'font.size': 6})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getSingleTifFromPath():\n",
    "    user_input = input('Enter the folder path to a sample GeoTIFF image:\\n')\n",
    "    test_path = f'{user_input}'\n",
    "    if os.path.exists(test_path): \n",
    "        return rio.open(os.path.join(test_path))\n",
    "    else:\n",
    "        print(\"The supplied folder path does not exist, please try again!\")\n",
    "        return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "test_image_full = getSingleTifFromPath()\n",
    "print(test_image_full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(test_image_full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_image_full.meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_image_full.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scale_min_max_contours(m, r_min, r_max, t_min, t_max):\n",
    "#     print(m[0], m[1])\n",
    "    x_scale = (m[0] - t_min)/(t_max - t_min) * (r_max - r_min) + r_min\n",
    "    y_scale = (m[1] - t_min)/(t_max - t_min) * (r_max - r_min) + r_min\n",
    "    return (int(x_scale), int(y_scale))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_contour_scale(dim, points):\n",
    "    t_min = round(np.array(points).min(),0)\n",
    "    t_max = round(np.array(points).max(),0)\n",
    "    return 0, dim-1, int(t_min), int(t_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_flat_contour_array(img):\n",
    "    width, length = img.shape[0], img.shape[1]\n",
    "    zero_matrix = np.zeros((width, length), int)\n",
    "    x, y = np.linspace(0, width, width), np.linspace(0, length, length)\n",
    "    X, Y = np.meshgrid(x, y)\n",
    "    cs = plt.contour(X, Y, img)\n",
    "    contour_points = []\n",
    "    for i, collection in enumerate(cs.collections[1:]):\n",
    "        for path in collection.get_paths():\n",
    "            if path.to_polygons():\n",
    "                for npoly, polypoints in enumerate(path.to_polygons()):\n",
    "                    for polypoint in polypoints:\n",
    "                        contour_points.append((polypoint[0], polypoint[1]))\n",
    "                        \n",
    "    r_min, r_max, t_min, t_max = get_contour_scale(width, contour_points)\n",
    "    print(r_min, r_max, t_min, t_max)\n",
    "    contour_points = np.apply_along_axis(scale_min_max_contours, 1, contour_points, r_min, r_max, t_min, t_max)\n",
    "    \n",
    "#     contour_points = np.apply_along_axis(scale_min_max_contours, 1, contour_points, 0, 99, 0, 100)\n",
    "    for point in contour_points:\n",
    "        zero_matrix[point[0]][point[1]] = 1\n",
    "    plt.imshow(zero_matrix.T)\n",
    "    plt.show()\n",
    "#     plt.hist(zero_matrix.T)\n",
    "#     plt.show()\n",
    "    return zero_matrix.T.flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotRGBNIR(img_path):\n",
    "    rio_img = rio.open(img_path)\n",
    "    plt.subplot(221),plt.imshow(rio_img.read(1), 'Reds'),plt.title(\"Red Band\")\n",
    "    plt.subplot(222),plt.imshow(rio_img.read(2), 'Blues'),plt.title(\"Green Band\")\n",
    "    plt.subplot(223),plt.imshow(rio_img.read(3), 'Greens'),plt.title(\"Blue Band\")\n",
    "    plt.subplot(224),plt.imshow(rio_img.read(4), 'Spectral_r'),plt.title(\"NIR Band\")\n",
    "    plt.subplots_adjust(hspace=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getFolderPath():\n",
    "    user_input = input('Enter the folder path to where sample TIF images are located:\\n')\n",
    "    test_path = f'{user_input}'\n",
    "    if os.path.exists(test_path): \n",
    "        return test_path\n",
    "    else:\n",
    "        print(\"The supplied folder path does not exist, please try again!\")\n",
    "        return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this folder path is expecting the \"annotation subset\" folder\n",
    "folder_path = getAllImagesFromPath()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is a sample of spliced HRO images representing nature and man made structures\n",
    "sample_nature = '10seg730865_700-1400.tif'\n",
    "sample_nature2 = '10seg985865_2100-4800.tif'\n",
    "sample_nature3 = '10seg985850_2400-3100.tif'\n",
    "sample_pool = '10seg985850_300-3700.tif'\n",
    "sample_house = '10seg745865_4700-1000.tif'\n",
    "sample_road = '10seg745880_400-1200.tif'\n",
    "sample_parkinglot = '10seg745880_4600-4900.tif'\n",
    "sample_solar = '10seg895850_2400-200.tif'\n",
    "\n",
    "sample_set = [['nature1', sample_nature, '1'],\n",
    "              ['nature2', sample_nature2, '1'],\n",
    "              ['nature3', sample_nature3, '1'],\n",
    "              ['house', sample_house, '0'],\n",
    "              ['road', sample_road, '0'],\n",
    "              ['parking', sample_parkinglot, '0'],\n",
    "              ['solar', sample_solar, '0'],\n",
    "              ['pool', sample_pool, '0']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotRGBNIR(os.path.join(folder_path, 'cluster_0', sample_solar))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotRGBNIR(os.path.join(folder_path, 'cluster_1', sample_nature))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_img = rio.open(os.path.join(folder_path, 'cluster_0', sample_house))\n",
    "for ii in range(1, sample_img.meta['count']+1):\n",
    "    print('-'*80)\n",
    "    print(np.percentile(sample_img.read(ii), 5))\n",
    "    print(np.median(sample_img.read(ii)))\n",
    "    print(np.percentile(sample_img.read(ii), 95))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(sample_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(sample_img.read(4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_img = rio.open(os.path.join(folder_path, 'cluster_0', sample_solar))\n",
    "for ii in range(1, sample_img.meta['count']+1):\n",
    "    print('-'*80)\n",
    "    print(np.percentile(sample_img.read(ii), 5))\n",
    "    print(np.median(sample_img.read(ii)))\n",
    "    print(np.percentile(sample_img.read(ii), 95))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(sample_img.read(4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(sample_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_img = rio.open(os.path.join(folder_path, 'cluster_1', sample_nature2))\n",
    "for ii in range(1, sample_img.meta['count']+1):\n",
    "    print('-'*80)\n",
    "    print(np.percentile(sample_img.read(ii), 5))\n",
    "    print(np.median(sample_img.read(ii)))\n",
    "    print(np.percentile(sample_img.read(ii), 95))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(sample_img.read(4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(sample_img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate pixelwise histogram of average pixel value by specified color band"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generatePixelRowDistribution(path, cluster, tif, band):\n",
    "    image = rio.open(os.path.join(path, cluster, tif))\n",
    "    n = image.meta['width']\n",
    "    pixel_dist = [0]*n\n",
    "    for row in image.read(band):\n",
    "        for cell in range(len(row)):\n",
    "            pixel_dist[cell] += row[cell]\n",
    "    pixel_dist = [round(val/n, 0) for val in pixel_dist]\n",
    "    plt.title(f'Row-wise pixel value distribution for band {band} - {tif} - {cluster}')\n",
    "    plt.plot(pixel_dist)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for sample in sample_set:\n",
    "    generatePixelRowDistribution(folder_path, f'cluster_{sample[2]}', sample[1], 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using PCA to determine if components can be built from image contours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flat_contours = []\n",
    "for sample in sample_set:\n",
    "    print(f'Opening image and gathering contours for {sample[0]}')\n",
    "    rio_img = rio.open(os.path.join(folder_path, f'cluster_{sample[2]}', sample[1]))\n",
    "    show(rio_img)\n",
    "    cv_img = cv2.imread(os.path.join(folder_path, f'cluster_{sample[2]}', sample[1]), 2)\n",
    "    flat_contours.append(build_flat_contour_array(cv_img))\n",
    "    rio_img.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for li in flat_contours:\n",
    "    plt.hist(li, bins=10)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.array(flat_contours).T.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "pca = PCA(n_components=8)\n",
    "pca.fit(np.array(flat_contours).T)\n",
    "print(plt.plot(pca.explained_variance_ratio_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pca.explained_variance_ratio_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_contours = []\n",
    "tif_dir = os.path.join(folder_path, 'cluster_0')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Warning, this next cell takes a long time to run because it generates contours with CV2 library on potentially thousands of TIF images, depending on the folder you specify"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for image in [tif for tif in os.listdir(tif_dir) if tif.__contains__('.tif')]:\n",
    "    cv_img = cv2.imread(os.path.join(tif_dir, image), 2)\n",
    "    all_contours.append(build_flat_contour_array(cv_img))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca2 = PCA(.8)\n",
    "pca_mat = pca2.fit_transform(np.array(all_contours))\n",
    "print(pca2.explained_variance_ratio_)\n",
    "plt.plot(pca2.explained_variance_ratio_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(pca2.explained_variance_ratio_[0:1400])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(pca2.explained_variance_ratio_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca2.components_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Everything below this point is experimental and code gaps exist, run at your own risk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "df = pd.DataFrame()\n",
    "df['contours'] = flat_contours\n",
    "df = pd.concat([df[col].apply(pd.Series) for col in df.columns], axis=1, ignore_index=True)\n",
    "col_count = df.shape[1]\n",
    "df['target'] = [0,0,0,1,1,1,1,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head(8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "km = KMeans(n_clusters=2)\n",
    "contour_cols = np.linspace(0,9999, 10000)\n",
    "X = df[contour_cols].values\n",
    "km.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['cluster'] = km.labels_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['diff'] = df.apply(lambda row: 1 if row['target'] == row['cluster'] else 0, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_set[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_list = [sample[1] for sample in sample_set]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['file'] = file_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.loc[df['file']=='10seg730865_700-1400.tif']['file']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rio_solar = rio.open(os.path.join(folder_path, sample_solar))\n",
    "rio_forest = rio.open(os.path.join(folder_path, sample_nature))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(5,5))\n",
    "show((rio_solar, 1), cmap='Greys_r', interpolation='none', ax=ax)\n",
    "show((rio_solar, 1), contour=True, ax=ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(5,5))\n",
    "show((rio_forest, 1), cmap='Greys_r', interpolation='none', ax=ax)\n",
    "show((rio_forest, 1), contour=True, ax=ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = np.linspace(0,100, 100), np.linspace(0,100, 100)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "plt.contourf(X, Y, img, 20, cmap='RdGy')\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "contours = plt.contour(X, Y, img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_cv = cv2.imread(os.path.join(folder_path, sample_solar))\n",
    "imgray = cv2.cvtColor(img_cv, cv2.COLOR_BGR2GRAY)\n",
    "ret, thresh = cv2.threshold(imgray, 127, 255, 0)\n",
    "contours, threshold = cv2.findContours(thresh, cv2.THRESH_TOZERO, cv2.CHAIN_APPROX_SIMPLE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv2.drawContours(img_cv, im2, -1, (0,255,0), 3)\n",
    "cv2.imshow('Contours', img_cv)\n",
    "cv2.waitKey(0)\n",
    "# cv2.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cs = plt.contourf(rio_forest, levels=[10, 30, 50],\n",
    "#     colors=['#808080', '#A0A0A0', '#C0C0C0'], extend='both')\n",
    "# cs.cmap.set_over('red')\n",
    "# cs.cmap.set_under('blue')\n",
    "# cs.changed()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# im = cv2.imread('test.jpg')\n",
    "# imgray = cv2.cvtColor(im,cv2.COLOR_BGR2GRAY)\n",
    "# ret,thresh = cv2.threshold(imgray,127,255,0)\n",
    "# image, contours, hierarchy = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forest_img = cv2.imread(os.path.join(folder_path, sample_nature), 2)\n",
    "edges = cv2.Canny(forest_img,25,255)\n",
    "\n",
    "plt.subplot(121),plt.imshow(forest_img,cmap = 'gray')\n",
    "plt.title('Original Image'), plt.xticks([]), plt.yticks([])\n",
    "plt.subplot(122),plt.imshow(edges,cmap = 'gray')\n",
    "plt.title('Edge Image'), plt.xticks([]), plt.yticks([])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "build_flat_contour_array(forest_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# >>> cs = plt.contour(x,y,m, [9.5])\n",
    "# >>> cs.collections[0].get_paths()\n",
    "# p = cs.collections[0].get_paths()[0]\n",
    "# v = p.vertices\n",
    "# x = v[:,0]\n",
    "# y = v[:,1]\n",
    "\n",
    "# cs = plt.contour(X, Y, forest_img)\n",
    "# p = cs.collections.get_paths()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parking_img = cv2.imread(os.path.join(folder_path, sample_parkinglot), 2)\n",
    "edges = cv2.Canny(parking_img,25,255)\n",
    "\n",
    "plt.subplot(121),plt.imshow(parking_img,cmap = 'gray')\n",
    "plt.title('Original Image'), plt.xticks([]), plt.yticks([])\n",
    "plt.subplot(122),plt.imshow(edges,cmap = 'gray')\n",
    "plt.title('Edge Image'), plt.xticks([]), plt.yticks([])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cs = plt.contour(X, Y, parking_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "build_flat_contour_array(parking_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "house_img = cv2.imread(os.path.join(folder_path, sample_house), 2)\n",
    "edges = cv2.Canny(house_img,25,255)\n",
    "\n",
    "plt.subplot(121),plt.imshow(house_img,cmap = 'gray')\n",
    "plt.title('Original Image'), plt.xticks([]), plt.yticks([])\n",
    "plt.subplot(122),plt.imshow(edges,cmap = 'gray')\n",
    "plt.title('Edge Image'), plt.xticks([]), plt.yticks([])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cs = plt.contour(X, Y, house_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "build_flat_contour_array(house_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "road_img = cv2.imread(os.path.join(folder_path, sample_road), 2)\n",
    "edges = cv2.Canny(road_img,25,200)\n",
    "\n",
    "plt.subplot(121),plt.imshow(road_img,cmap = 'gray')\n",
    "plt.title('Original Image'), plt.xticks([]), plt.yticks([])\n",
    "plt.subplot(122),plt.imshow(edges,cmap = 'gray')\n",
    "plt.title('Edge Image'), plt.xticks([]), plt.yticks([])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cs = plt.contour(X, Y, road_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "build_flat_contour_array(road_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "solar_img = cv2.imread(os.path.join(folder_path, sample_solar), 2)\n",
    "edges = cv2.Canny(solar_img,100,255)\n",
    "\n",
    "plt.subplot(121),plt.imshow(solar_img,cmap = 'gray')\n",
    "plt.title('Original Image'), plt.xticks([]), plt.yticks([])\n",
    "plt.subplot(122),plt.imshow(edges,cmap = 'gray')\n",
    "plt.title('Edge Image'), plt.xticks([]), plt.yticks([])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cs = plt.contour(X, Y, solar_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "build_flat_contour_array(solar_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ret,thresh1 = cv2.threshold(img,50,255,cv2.THRESH_BINARY)\n",
    "ret,thresh2 = cv2.threshold(img,50,255,cv2.THRESH_BINARY_INV)\n",
    "ret,thresh3 = cv2.threshold(img,50,255,cv2.THRESH_TRUNC)\n",
    "ret,thresh4 = cv2.threshold(img,50,255,cv2.THRESH_TOZERO)\n",
    "ret,thresh5 = cv2.threshold(img,50,255,cv2.THRESH_TOZERO_INV)\n",
    "\n",
    "titles = ['Original Image','BINARY','BINARY_INV','TRUNC','TOZERO','TOZERO_INV']\n",
    "images = [img, thresh1, thresh2, thresh3, thresh4, thresh5]\n",
    "\n",
    "for i in range(6):\n",
    "    plt.subplot(2,3,i+1),plt.imshow(images[i],'gray')\n",
    "    plt.title(titles[i])\n",
    "    plt.xticks([]),plt.yticks([])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thresh3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ret,th1 = cv2.threshold(img,50,255,cv2.THRESH_BINARY)\n",
    "th2 = cv2.adaptiveThreshold(img,255,cv2.ADAPTIVE_THRESH_MEAN_C,\\\n",
    "            cv2.THRESH_BINARY,11,2)\n",
    "th3 = cv2.adaptiveThreshold(img,255,cv2.ADAPTIVE_THRESH_GAUSSIAN_C,\\\n",
    "            cv2.THRESH_BINARY,11,2)\n",
    "\n",
    "titles = ['Original Image', 'Global Thresholding (v = 127)',\n",
    "            'Adaptive Mean Thresholding', 'Adaptive Gaussian Thresholding']\n",
    "images = [img, th1, th2, th3]\n",
    "\n",
    "for i in range(4):\n",
    "    plt.subplot(2,2,i+1),plt.imshow(images[i],'gray')\n",
    "    plt.title(titles[i])\n",
    "    plt.xticks([]),plt.yticks([])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# global thresholding\n",
    "ret1,th1 = cv2.threshold(img,127,255,cv2.THRESH_BINARY)\n",
    "\n",
    "# Otsu's thresholding\n",
    "ret2,th2 = cv2.threshold(img,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)\n",
    "\n",
    "# Otsu's thresholding after Gaussian filtering\n",
    "blur = cv2.GaussianBlur(img,(5,5),0)\n",
    "ret3,th3 = cv2.threshold(blur,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)\n",
    "\n",
    "# plot all the images and their histograms\n",
    "images = [img, 0, th1,\n",
    "          img, 0, th2,\n",
    "          blur, 0, th3]\n",
    "titles = ['Original Noisy Image','Histogram','Global Thresholding (v=127)',\n",
    "          'Original Noisy Image','Histogram',\"Otsu's Thresholding\",\n",
    "          'Gaussian filtered Image','Histogram',\"Otsu's Thresholding\"]\n",
    "\n",
    "for i in range(3):\n",
    "    plt.subplot(3,3,i*3+1),plt.imshow(images[i*3],'gray')\n",
    "    plt.title(titles[i*3]), plt.xticks([]), plt.yticks([])\n",
    "    plt.subplot(3,3,i*3+2),plt.hist(images[i*3].ravel(),256)\n",
    "    plt.title(titles[i*3+1]), plt.xticks([]), plt.yticks([])\n",
    "    plt.subplot(3,3,i*3+3),plt.imshow(images[i*3+2],'gray')\n",
    "    plt.title(titles[i*3+2]), plt.xticks([]), plt.yticks([])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for sample in sample_set:\n",
    "#     im = rio.open(os.path.join(folder_path, sample[1]))\n",
    "#     show_hist(im, bins=50, lw=0.0, stacked=False, alpha=0.3,\n",
    "#          histtype='stepfilled', title=f'{sample[0]} Histogram')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_path = 'D://Raster//HRO//parsed_ortho_files//bay_area_ca_0.3m_tiles//'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_files = os.listdir(image_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_forest = '10seg745835_0-1100.tif'\n",
    "test_house = '10seg925850_1100-100.tif'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_forest_image = rio.open(os.path.join(image_path, test_forest))\n",
    "test_house_image = rio.open(os.path.join(image_path, test_house))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.imshow(test_image.read(2))\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(test_forest_image)\n",
    "# show((src, 2), cmap='viridis')\n",
    "# show(src.read(2), transform=src.transform, cmap='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show((test_forest_image,1), cmap='viridis')\n",
    "# show(test_forest_image.read(2), transform=test_forest_image.transform, cmap='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(test_house_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show((test_house_image,1), cmap='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show((test_house_image,2), cmap='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show((test_house_image,3), cmap='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show((test_house_image,4), cmap='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(5,5))\n",
    "show((test_house_image, 1), cmap='Greys_r', interpolation='none', ax=ax)\n",
    "show((test_house_image, 1), contour=True, ax=ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(5,5))\n",
    "show((test_forest_image, 1), cmap='Greys_r', interpolation='none', ax=ax)\n",
    "show((test_forest_image, 1), contour=True, ax=ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# >>> fig, ax = pyplot.subplots(1, figsize=(12, 12))\n",
    "# >>> show((src, 1), cmap='Greys_r', interpolation='none', ax=ax)\n",
    "# <matplotlib.axes._subplots.AxesSubplot object at 0x...>\n",
    "# >>> show((src, 1), contour=True, ax=ax)\n",
    "# <matplotlib.axes._subplots.AxesSubplot object at 0x...>\n",
    "# >>> pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rasterio.plot import show_hist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_hist(test_house_image, bins=50, lw=0.0, stacked=False, alpha=0.3,\n",
    "         histtype='stepfilled', title='House Histogram')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_hist(test_forest_image, bins=50, lw=0.0, stacked=False, alpha=0.3,\n",
    "         histtype='stepfilled', title='Forest Histogram')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# >>> from rasterio.plot import show_hist\n",
    "# >>> show_hist(\n",
    "# ...     src, bins=50, lw=0.0, stacked=False, alpha=0.3,\n",
    "# ...     histtype='stepfilled', title=\"Histogram\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_image.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import measure\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_house_image.read(1)/255"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = np.ogrid[-np.pi:np.pi:100j, -np.pi:np.pi:100j]\n",
    "r = np.sin(np.exp((np.sin(x)**3 + np.cos(y)**2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "contours = measure.find_contours(r, .8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "house_contours = measure.find_contours(test_house_image.read(1)/255, .2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# house_contours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.imshow(test_house_image.read(1), cmap=plt.cm.gray)\n",
    "\n",
    "for contour in house_contours:\n",
    "    ax.plot(contour[:, 1], contour[:, 0], linewidth=2)\n",
    "\n",
    "ax.axis('image')\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forest_contours = measure.find_contours(test_forest_image.read(1)/255, .2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.imshow(test_forest_image.read(1), cmap=plt.cm.gray)\n",
    "\n",
    "for contour in forest_contours:\n",
    "    ax.plot(contour[:, 1], contour[:, 0], linewidth=2)\n",
    "\n",
    "ax.axis('image')\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nature_list = ['10seg925850_1000-4600.tif', '10seg745835_0-1100.tif', '10seg925850_1100-1300.tif', '10seg925850_1100-2000.tif', '10seg925850_1100-3900.tif', '10seg955835_4000-100.tif']\n",
    "house_list = ['10seg925850_1100-300.tif', '10seg925850_1100-100.tif', '10seg925850_1200-300.tif', '10seg925850_1200-400.tif', '10seg955835_3900-3900.tif','10seg955835_3900-3800.tif']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_contours(in_arr, const=.2):\n",
    "    return measure.find_contours(in_arr/255, const)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_img_contours(img, band, contours):\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.imshow(img.read(band), cmap=plt.cm.gray)\n",
    "\n",
    "    for contour in contours:\n",
    "        ax.plot(contour[:, 1], contour[:, 0], linewidth=2)\n",
    "\n",
    "    ax.axis('image')\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    plt.show()"
   ]
  },
  {
   "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",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}