1_Road_classification_VGG.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a full explanation of the code, visit http://ataspinar.com/2017/12/04/using-convolutional-neural-networks-to-detect-features-in-sattelite-images/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.insert(0,'..')\n",
    "\n",
    "import os\n",
    "import json\n",
    "import pickle\n",
    "import requests\n",
    "import shutil\n",
    "import random\n",
    "import numpy as np\n",
    "\n",
    "from collections import defaultdict\n",
    "from scipy import ndimage\n",
    "from datetime import date\n",
    "from datetime import datetime\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#We are using owslib to download images from a WMS Service\n",
    "#install with 'pip install owslib'\n",
    "\n",
    "import owslib\n",
    "\n",
    "from owslib.wms import WebMapService\n",
    "\n",
    "#pyshp is necessary for loading and saving shapefiles\n",
    "#install with 'pip install pyshp'\n",
    "import shapefile\n",
    "\n",
    "# Install opencv with 'pip install opencv-python'\n",
    "import cv2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_min = 90000\n",
    "y_min = 427000\n",
    "dx, dy = 200, 200\n",
    "no_tiles_x = 100\n",
    "no_tiles_y = 100\n",
    "total_no_tiles = no_tiles_x * no_tiles_y\n",
    "\n",
    "x_max = x_min + no_tiles_x * dx\n",
    "y_max = y_min + no_tiles_y * dy\n",
    "bounding_box = [x_min, y_min, x_max, y_max]\n",
    "\n",
    "TILE_FOLDER = \"./datasets/image_tiles_200/\"\n",
    "URL_TILES = \"https://geodata.nationaalgeoregister.nl/luchtfoto/rgb/wms?request=GetCapabilities\"\n",
    "\n",
    "URL_SHP = 'https://www.rijkswaterstaat.nl/apps/geoservices/geodata/dmc/nwb-wegen/geogegevens/shapefile/Nederland_totaal/01-08-2017/Wegvakken/Wegvakken.shp'\n",
    "URL_PRF = 'https://www.rijkswaterstaat.nl/apps/geoservices/geodata/dmc/nwb-wegen/geogegevens/shapefile/Nederland_totaal/01-08-2017/Wegvakken/Wegvakken.prj'\n",
    "URL_DBF = 'https://www.rijkswaterstaat.nl/apps/geoservices/geodata/dmc/nwb-wegen/geogegevens/shapefile/Nederland_totaal/01-08-2017/Wegvakken/Wegvakken.dbf'\n",
    "URL_SHX = 'https://www.rijkswaterstaat.nl/apps/geoservices/geodata/dmc/nwb-wegen/geogegevens/shapefile/Nederland_totaal/01-08-2017/Wegvakken/Wegvakken.shx'\n",
    "\n",
    "URLS_SHAPEFILES = [URL_SHP, URL_PRF, URL_DBF, URL_SHX]\n",
    "\n",
    "DATA_FOLDER = \"./data/nwb_wegvakken/\"\n",
    "\n",
    "json_filename = DATA_FOLDER + '2017_09_wegvakken.json'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Downloading the image tiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wms = WebMapService(URL_TILES, version='1.1.1')\n",
    "\n",
    "if not os.path.exists(TILE_FOLDER):\n",
    "    os.makedirs(TILE_FOLDER)\n",
    "\n",
    "for ii in range(25,no_tiles_x):\n",
    "    print(ii)\n",
    "    for jj in range(0,no_tiles_y):\n",
    "        ll_x_ = x_min + ii*dx\n",
    "        ll_y_ = y_min + jj*dy\n",
    "        bbox = (ll_x_, ll_y_, ll_x_ + dx, ll_y_ + dy) \n",
    "        img = wms.getmap(layers=['Actueel_ortho25'], srs='EPSG:28992', bbox=bbox, size=(256, 256), format='image/jpeg', transparent=True)\n",
    "        filename = \"{}{}_{}_{}_{}.jpg\".format(TILE_FOLDER, bbox[0], bbox[1], bbox[2], bbox[3])\n",
    "        out = open(filename, 'wb')\n",
    "        out.write(img.read())\n",
    "        out.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1b. Downloading the shapefiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(DATA_FOLDER):\n",
    "    os.makedirs(DATA_FOLDER)\n",
    "\n",
    "for url in URLS_SHAPEFILES:\n",
    "    filename = url.split('/')[-1]\n",
    "    print(\"Downloading file {}\".format(filename))\n",
    "    r = requests.get(url, stream=True)\n",
    "    if r.status_code == 200:\n",
    "        with open(DATA_FOLDER + filename, 'wb') as f:\n",
    "            r.raw.decode_content = True\n",
    "            shutil.copyfileobj(r.raw, f)     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Loading shapefile and converting to (GEO)Json"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def json_serial(obj):\n",
    "    \"\"\"JSON serializer for objects not serializable by default json code\"\"\"\n",
    "\n",
    "    if isinstance(obj, (datetime, date)):\n",
    "        serial = obj.isoformat()\n",
    "        return serial\n",
    "    if isinstance(obj, bytes):\n",
    "        return {'__class__': 'bytes',\n",
    "                '__value__': list(obj)}\n",
    "    raise TypeError (\"Type %s not serializable\" % type(obj))\n",
    "\n",
    "reader = shapefile.Reader(DATA_FOLDER + 'Wegvakken.shp')\n",
    "fields = reader.fields[1:]\n",
    "field_names = [field[0] for field in fields]\n",
    "\n",
    "buffer = []\n",
    "for sr in reader.shapeRecords()[:500000]:\n",
    "    atr = dict(zip(field_names, sr.record))\n",
    "    geom = sr.shape.__geo_interface__\n",
    "    buffer.append(dict(type=\"Feature\", geometry=geom, properties=atr)) \n",
    "\n",
    "\n",
    "json_file = open(json_filename, \"w\")\n",
    "json_file.write(json.dumps({\"type\": \"FeatureCollection\", \"features\": buffer}, indent=2, default=json_serial) + \"\\n\")\n",
    "json_file.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Declaring some variables and methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dict_roadtype = {\n",
    "    \"G\": 'Gemeente',\n",
    "    \"R\": 'Rijk',\n",
    "    \"P\": 'Provincie',\n",
    "    \"W\": 'Waterschap',\n",
    "    'T': 'Andere wegbeheerder',\n",
    "    '' : 'leeg'\n",
    "}\n",
    "\n",
    "dict_roadtype_to_color = {\n",
    "    \"G\": 'red',\n",
    "    \"R\": 'blue',\n",
    "    \"P\": 'green',\n",
    "    \"W\": 'magenta',\n",
    "    'T': 'yellow',\n",
    "    '' : 'leeg'\n",
    "}\n",
    "\n",
    "FEATURES_KEY = 'features'\n",
    "PROPERTIES_KEY = 'properties'\n",
    "GEOMETRY_KEY = 'geometry'\n",
    "COORDINATES_KEY = 'coordinates'\n",
    "WEGSOORT_KEY = 'WEGBEHSRT'\n",
    "\n",
    "MINIMUM_NO_POINTS_PER_TILE = 4\n",
    "POINTS_PER_METER = 0.1\n",
    "\n",
    "INPUT_FOLDER_TILES = './datasets/image_tiles_200/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_to_dict(d1, d2, coordinates, rtype):\n",
    "    coordinate_ll_x = int((coordinates[0] // dx)*dx)\n",
    "    coordinate_ll_y = int((coordinates[1] // dy)*dy)\n",
    "    coordinate_ur_x = int((coordinates[0] // dx)*dx + dx)\n",
    "    coordinate_ur_y = int((coordinates[1] // dy)*dy + dy)\n",
    "    tile = \"{}_{}_{}_{}.jpg\".format(coordinate_ll_x, coordinate_ll_y, coordinate_ur_x, coordinate_ur_y)\n",
    "    \n",
    "    rel_coord_x = (coordinates[0] - coordinate_ll_x) / dx\n",
    "    rel_coord_y = (coordinates[1] - coordinate_ll_y) / dy\n",
    "    value = (rtype, rel_coord_x, rel_coord_y)\n",
    "    d1[tile].append(value)\n",
    "    d2[rtype].add(tile)\n",
    "\n",
    "def coord_is_in_bb(coord, bb):\n",
    "    x_min = bb[0]\n",
    "    y_min = bb[1]\n",
    "    x_max = bb[2]\n",
    "    y_max = bb[3]\n",
    "    return coord[0] > x_min and coord[0] < x_max and coord[1] > y_min and coord[1] < y_max\n",
    "\n",
    "def retrieve_roadtype(elem):\n",
    "    return elem[PROPERTIES_KEY][WEGSOORT_KEY]\n",
    "   \n",
    "def retrieve_coordinates(elem):\n",
    "    return elem[GEOMETRY_KEY][COORDINATES_KEY]\n",
    "\n",
    "def eucledian_distance(p1, p2):\n",
    "    diff = np.array(p2)-np.array(p1)\n",
    "    return np.linalg.norm(diff)\n",
    "\n",
    "def calculate_intermediate_points(p1, p2, no_points):\n",
    "    dx = (p2[0] - p1[0]) / (no_points + 1)\n",
    "    dy = (p2[1] - p1[1]) / (no_points + 1)\n",
    "    return [[p1[0] + i * dx, p1[1] +  i * dy] for i in range(1, no_points+1)]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Map contents of shapefile to the tiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename_wegvakken = json_filename\n",
    "dict_wegvakken = json.load(open(filename_wegvakken))[FEATURES_KEY]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_tile_contents = defaultdict(list)\n",
    "d_roadtype_tiles = defaultdict(set)\n",
    "    \n",
    "for elem in dict_wegvakken:\n",
    "    coordinates = retrieve_coordinates(elem)\n",
    "    rtype = retrieve_roadtype(elem)\n",
    "    coordinates_in_bb = [coord for coord in coordinates if coord_is_in_bb(coord, bounding_box)]\n",
    "    if len(coordinates_in_bb)==1:\n",
    "        coord = coordinates_in_bb[0]\n",
    "        add_to_dict(d_tile_contents, d_roadtype_tiles, coord, rtype)\n",
    "    if len(coordinates_in_bb)>1:\n",
    "        add_to_dict(d_tile_contents, d_roadtype_tiles, coordinates_in_bb[0], rtype)\n",
    "        for ii in range(1,len(coordinates_in_bb)):\n",
    "            previous_coord = coordinates_in_bb[ii-1]\n",
    "            coord = coordinates_in_bb[ii]\n",
    "            add_to_dict(d_tile_contents, d_roadtype_tiles, coord, rtype)\n",
    "            \n",
    "            dist = eucledian_distance(previous_coord, coord)\n",
    "            no_intermediate_points = int(dist*POINTS_PER_METER)           \n",
    "            intermediate_coordinates = calculate_intermediate_points(previous_coord, coord, no_intermediate_points)\n",
    "            for intermediate_coord in intermediate_coordinates:\n",
    "                add_to_dict(d_tile_contents, d_roadtype_tiles, intermediate_coord, rtype)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4b. Visualize results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = 95000\n",
    "y0 = 430000\n",
    "\n",
    "fig, axarr = plt.subplots(nrows=11,ncols=11, figsize=(16,16))\n",
    "\n",
    "for ii in range(0,11):\n",
    "    for jj in range(0,11):\n",
    "        ll_x = x0 + ii*dx\n",
    "        ll_y = y0 + jj*dy\n",
    "        ur_x = ll_x + dx\n",
    "        ur_y = ll_y + dy\n",
    "        tile = \"{}_{}_{}_{}.jpg\".format(ll_x, ll_y, ur_x, ur_y)\n",
    "        filename = INPUT_FOLDER_TILES + tile        \n",
    "        tile_contents = d_tile_contents[tile]\n",
    "        \n",
    "        ax = axarr[10-jj, ii]\n",
    "        image = ndimage.imread(filename)\n",
    "        rgb_image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)\n",
    "        ax.imshow(rgb_image)\n",
    "        ax.get_xaxis().set_visible(False)\n",
    "        ax.get_yaxis().set_visible(False)\n",
    "\n",
    "        for elem in tile_contents:\n",
    "            color = dict_roadtype_to_color[elem[0]]\n",
    "            x = elem[1]*256\n",
    "            y = (1-elem[2])*256\n",
    "            ax.scatter(x,y,c=color,s=10)\n",
    "\n",
    "plt.subplots_adjust(wspace=0, hspace=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = 95400\n",
    "y0 = 432000\n",
    "POINTS_PER_METER = 0\n",
    "\n",
    "fig, axarr = plt.subplots(nrows=11,ncols=11, figsize=(16,16))\n",
    "\n",
    "for ii in range(0,11):\n",
    "    for jj in range(0,11):\n",
    "        ll_x = x0 + ii*dx\n",
    "        ll_y = y0 + jj*dy\n",
    "        ur_x = ll_x + dx\n",
    "        ur_y = ll_y + dy\n",
    "        tile = \"{}_{}_{}_{}.jpg\".format(ll_x, ll_y, ur_x, ur_y)\n",
    "        filename = INPUT_FOLDER_TILES + tile\n",
    "        tile_contents = d_tile_contents[tile]\n",
    "        \n",
    "        ax = axarr[10-jj, ii]\n",
    "        image = ndimage.imread(filename)\n",
    "        rgb_image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)\n",
    "        ax.imshow(rgb_image)\n",
    "        ax.get_xaxis().set_visible(False)\n",
    "        ax.get_yaxis().set_visible(False)\n",
    "\n",
    "        for elem in tile_contents:\n",
    "            color = dict_roadtype_to_color[elem[0]]\n",
    "            x = elem[1]*256\n",
    "            y = (1-elem[2])*256\n",
    "            ax.scatter(x,y,c=color,s=10)\n",
    "\n",
    "plt.subplots_adjust(wspace=0, hspace=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4c. Some statistics about the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"There are {} tiles containing roads.\".format(len(d_tile_contents.keys())))\n",
    "\n",
    "for rtype in d_roadtype_tiles.keys():\n",
    "    roadtype = dict_roadtype[rtype]\n",
    "    no_tiles = len(d_roadtype_tiles[rtype])\n",
    "    print(\"Of roadtype {} ({}) there are {} tiles.\".format(rtype, roadtype, no_tiles))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Prepare dataset for CNN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def accuracy(predictions, labels):\n",
    "    return (100.0 * np.sum(np.argmax(predictions, 1) == np.argmax(labels, 1)) / predictions.shape[0])\n",
    "\n",
    "def onehot_encode_labels(labels):\n",
    "    list_possible_labels = list(np.unique(labels))\n",
    "    encoded_labels = map(lambda x: list_possible_labels.index(x), labels)\n",
    "    return encoded_labels\n",
    "\n",
    "def randomize(dataset, labels1, labels2, labels3):\n",
    "    permutation = np.random.permutation(dataset.shape[0])\n",
    "    print(permutation.shape)\n",
    "    print(dataset.shape)\n",
    "    print(labels1.shape)\n",
    "    print(labels2.shape)\n",
    "    print(labels3.shape)\n",
    "    randomized_dataset = dataset[permutation, :, :, :]\n",
    "    randomized_labels1 = labels1[permutation]\n",
    "    randomized_labels2 = labels2[permutation]\n",
    "    randomized_labels3 = labels3[permutation]\n",
    "    return randomized_dataset, randomized_labels1, randomized_labels2, randomized_labels3\n",
    "\n",
    "def one_hot_encode(np_array, num_unique_labels):\n",
    "    return (np.arange(num_unique_labels) == np_array[:,None]).astype(np.float32)\n",
    "\n",
    "def reformat_data(dataset, labels1, labels2, labels3):\n",
    "    dataset, labels1, labels2, labels3 = randomize(dataset, labels1, labels2, labels3)\n",
    "    num_unique_labels1 = len(np.unique(labels1))\n",
    "    num_unique_labels2 = len(np.unique(labels2))\n",
    "    labels1 = one_hot_encode(labels1, num_unique_labels1)\n",
    "    labels2 = one_hot_encode(labels2, num_unique_labels2)\n",
    "    print(dataset.shape)\n",
    "    print(labels1.shape)\n",
    "    print(labels2.shape)\n",
    "    print(labels3.shape)\n",
    "    return dataset, labels1, labels2, labels3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_width = 256\n",
    "image_height = 256\n",
    "image_depth = 3\n",
    "total_no_images = 10000\n",
    "\n",
    "image_files = os.listdir(INPUT_FOLDER_TILES)\n",
    "\n",
    "dataset = np.ndarray(shape=(total_no_images, image_width, image_height, image_depth), dtype=np.float32)\n",
    "labels_roadtype = []\n",
    "labels_roadpresence = np.ndarray(total_no_images, dtype=np.float32)\n",
    "labels_filename = []\n",
    "\n",
    "# Create labels\n",
    "for counter, image in enumerate(image_files):\n",
    "    filename = INPUT_FOLDER_TILES + image\n",
    "    labels_filename.append(image)\n",
    "    if image in list(d_tile_contents.keys()):\n",
    "        tile_contents = d_tile_contents[image]\n",
    "        roadtypes = sorted(list(set([elem[0] for elem in tile_contents])))\n",
    "        roadtype = \"_\".join(roadtypes)\n",
    "        labels_roadpresence[counter] = 1\n",
    "    else:\n",
    "        roadtype = ''\n",
    "        labels_roadpresence[counter] = 0\n",
    "    labels_roadtype.append(roadtype)\n",
    "\n",
    "    image_data = ndimage.imread(filename).astype(np.float32)\n",
    "    dataset[counter, :, :] = image_data\n",
    "    if counter % 1000 == 0:\n",
    "        print(\"{} images have been loaded.\".format(counter))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "labels_filename = np.array(labels_filename)\n",
    "labels_roadtype_ohe = np.array(list(onehot_encode_labels(labels_roadtype)))\n",
    "print(\"Randomizing dataset...\")\n",
    "dataset, labels_roadpresence, labels_roadtype_ohe, labels_filename = reformat_data(dataset, labels_roadpresence, labels_roadtype_ohe, labels_filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save dataset to pickle file\n",
    "start_train_dataset = 0\n",
    "start_valid_dataset = 120\n",
    "start_test_dataset = 160\n",
    "total_no_images = 200\n",
    "\n",
    "output_pickle_file = './data/sattelite_dataset.pickle'\n",
    "\n",
    "f = open(output_pickle_file, 'wb')\n",
    "save = {\n",
    "'train_dataset': dataset[start_train_dataset:start_valid_dataset,:,:,:],\n",
    "'train_labels_roadtype': labels_roadtype[start_train_dataset:start_valid_dataset],\n",
    "'train_labels_roadpresence': labels_roadpresence[start_train_dataset:start_valid_dataset],\n",
    "'train_labels_filename': labels_filename[start_train_dataset:start_valid_dataset],\n",
    "'valid_dataset': dataset[start_valid_dataset:start_test_dataset,:,:,:],\n",
    "'valid_labels_roadtype': labels_roadtype[start_valid_dataset:start_test_dataset],\n",
    "'valid_labels_roadpresence': labels_roadpresence[start_valid_dataset:start_test_dataset],\n",
    "'valid_labels_filename': labels_filename[start_valid_dataset:start_test_dataset],\n",
    "'test_dataset': dataset[start_test_dataset:total_no_images,:,:,:],\n",
    "'test_labels_roadtype': labels_roadtype[start_test_dataset:total_no_images],\n",
    "'test_labels_roadpresence': labels_roadpresence[start_test_dataset:total_no_images],\n",
    "'test_labels_filename': labels_filename[start_test_dataset:total_no_images]\n",
    "}\n",
    "pickle.dump(save, f, pickle.HIGHEST_PROTOCOL)\n",
    "f.close()\n",
    "\n",
    "print(\"\\nsaved dataset to {}\".format(output_pickle_file))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. The Convolutional neural network part"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ubuntu/anaconda3/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n",
      "/home/ubuntu/anaconda3/lib/python3.6/site-packages/matplotlib/__init__.py:1067: UserWarning: Duplicate key in file \"/home/ubuntu/.config/matplotlib/matplotlibrc\", line #2\n",
      "  (fname, cnt))\n",
      "/home/ubuntu/anaconda3/lib/python3.6/site-packages/matplotlib/__init__.py:1067: UserWarning: Duplicate key in file \"/home/ubuntu/.config/matplotlib/matplotlibrc\", line #3\n",
      "  (fname, cnt))\n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "from cnn_models.vggnet16 import *\n",
    "from sklearn.metrics import roc_curve\n",
    "from sklearn.metrics import auc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "pickle_file = './data/sattelite_dataset.pickle'\n",
    "f = open(pickle_file, 'rb')\n",
    "save = pickle.load(f)\n",
    "\n",
    "train_dataset = save['train_dataset'].astype(dtype = np.float32)\n",
    "train_labels = save['train_labels_roadpresence'].astype(dtype = np.float32)\n",
    "valid_dataset = save['valid_dataset'].astype(dtype = np.float32)\n",
    "valid_labels = save['valid_labels_roadpresence'].astype(dtype = np.float32)\n",
    "test_dataset = save['test_dataset'].astype(dtype = np.float32)\n",
    "test_labels = save['test_labels_roadpresence'].astype(dtype = np.float32)\n",
    "\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 1.]\n"
     ]
    }
   ],
   "source": [
    "num_labels = len(np.unique(train_labels))\n",
    "print(np.unique(train_labels))\n",
    "image_width = 256\n",
    "image_height = 256\n",
    "image_depth = 3\n",
    "num_steps = 10\n",
    "display_step = 1\n",
    "learning_rate = 0.0001\n",
    "batch_size = 16\n",
    "lambda_loss_amount = 0.0015"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "STARTING WITH SATTELITE\n",
      "<tf.Variable 'Variable_31:0' shape=(2,) dtype=float32_ref>\n",
      "WARNING:tensorflow:From <ipython-input-4-4f92221d41cc>:25: softmax_cross_entropy_with_logits (from tensorflow.python.ops.nn_ops) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "\n",
      "Future major versions of TensorFlow will allow gradients to flow\n",
      "into the labels input on backprop by default.\n",
      "\n",
      "See @{tf.nn.softmax_cross_entropy_with_logits_v2}.\n",
      "\n",
      "Initialized with learning_rate 0.0001  model \n",
      "step 0000 : loss is 2147222912.00, accuracy on training set 50.00 %, accuracy on test set 37.50 accuracy on valid set 18.75 %\n",
      "step 0001 : loss is 2657171712.00, accuracy on training set 56.25 %, accuracy on test set 68.75 accuracy on valid set 62.50 %\n",
      "step 0002 : loss is 1961304192.00, accuracy on training set 50.00 %, accuracy on test set 50.00 accuracy on valid set 68.75 %\n",
      "step 0003 : loss is 3072617216.00, accuracy on training set 50.00 %, accuracy on test set 50.00 accuracy on valid set 68.75 %\n",
      "step 0004 : loss is 3190644736.00, accuracy on training set 31.25 %, accuracy on test set 50.00 accuracy on valid set 56.25 %\n",
      "step 0005 : loss is 1370354432.00, accuracy on training set 56.25 %, accuracy on test set 62.50 accuracy on valid set 12.50 %\n",
      "step 0006 : loss is 2114763648.00, accuracy on training set 50.00 %, accuracy on test set 37.50 accuracy on valid set 50.00 %\n",
      "step 0007 : loss is 1925198336.00, accuracy on training set 56.25 %, accuracy on test set 50.00 accuracy on valid set 43.75 %\n",
      "step 0008 : loss is 694073024.00, accuracy on training set 68.75 %, accuracy on test set 50.00 accuracy on valid set 50.00 %\n",
      "step 0009 : loss is 471897152.00, accuracy on training set 56.25 %, accuracy on test set 56.25 accuracy on valid set 56.25 %\n"
     ]
    }
   ],
   "source": [
    "train_accuracies, test_accuracies, valid_accuracies = [], [], []\n",
    "train_preds, test_preds, valid_preds = [], [], []\n",
    " \n",
    "print(\"STARTING WITH SATTELITE\")\n",
    "graph = tf.Graph()\n",
    "with graph.as_default():\n",
    "    #1) First we put the input data in a tensorflow friendly form. \n",
    "    tf_train_dataset = tf.placeholder(tf.float32, shape=(batch_size, image_width, image_height, image_depth))\n",
    "    tf_train_labels = tf.placeholder(tf.float32, shape = (batch_size, num_labels))\n",
    "    tf_test_dataset = tf.placeholder(tf.float32, shape=(batch_size, image_width, image_height, image_depth))\n",
    "    tf_test_labels = tf.placeholder(tf.float32, shape = (batch_size, num_labels))\n",
    "    tf_valid_dataset = tf.placeholder(tf.float32, shape=(batch_size, image_width, image_height, image_depth))\n",
    "    tf_valid_labels = tf.placeholder(tf.float32, shape = (batch_size, num_labels))\n",
    " \n",
    "    #2) Then, the weight matrices and bias vectors are initialized\n",
    "    variables = variables_vggnet16()\n",
    "    print(variables[\"b16\"])\n",
    "    #3. The model used to calculate the logits (predicted labels)\n",
    "    model = model_vggnet16\n",
    "    \n",
    "    logits = model(tf_train_dataset, variables)\n",
    " \n",
    "    #4. then we compute the softmax cross entropy between the logits and the (actual) labels\n",
    "    l2 = lambda_loss_amount * sum(tf.nn.l2_loss(tf_var) for tf_var in tf.trainable_variables())\n",
    "    loss = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(logits=logits, labels=tf_train_labels)) + l2\n",
    " \n",
    "    #learning_rate = tf.train.exponential_decay(0.05, global_step, 1000, 0.85, staircase=True)\n",
    "    #5. The optimizer is used to calculate the gradients of the loss function \n",
    "    optimizer = tf.train.AdamOptimizer(learning_rate=learning_rate).minimize(loss)\n",
    " \n",
    "    # Predictions for the training, validation, and test data.\n",
    "    train_prediction = tf.nn.softmax(logits)\n",
    "    test_prediction = tf.nn.softmax(model(tf_test_dataset, variables))\n",
    "    valid_prediction = tf.nn.softmax(model(tf_valid_dataset, variables))\n",
    " \n",
    " \n",
    "with tf.Session(graph=graph) as session:\n",
    "    test_counter = 0\n",
    "    tf.global_variables_initializer().run()\n",
    "    print('Initialized with learning_rate', learning_rate, \" model \")\n",
    "    for step in range(num_steps):\n",
    "        #Since we are using stochastic gradient descent, we are selecting  small batches from the training dataset,\n",
    "        #and training the convolutional neural network each time with a batch. \n",
    "        offset = (step * batch_size) % (train_labels.shape[0] - batch_size)\n",
    "        batch_data = train_dataset[offset:(offset + batch_size), :,  :]\n",
    "        batch_labels = train_labels[offset:(offset + batch_size), :]\n",
    "        feed_dict = {tf_train_dataset : batch_data, tf_train_labels : batch_labels}\n",
    "        \n",
    "        _, l, predictions = session.run([optimizer, loss, train_prediction], feed_dict=feed_dict)\n",
    "        train_accuracy = accuracy(predictions, batch_labels)\n",
    "        train_accuracies.append(train_accuracy)\n",
    "        train_preds.append(predictions)\n",
    " \n",
    "        if step % display_step == 0:\n",
    "            offset2 = (test_counter * batch_size) % (test_labels.shape[0] - batch_size)\n",
    "            test_dataset_batch = test_dataset[offset2:(offset2 + batch_size), :, :]\n",
    "            test_labels_batch = test_labels[offset2:(offset2 + batch_size), :]\n",
    "            feed_dict2 = {tf_test_dataset : test_dataset_batch, tf_test_labels : test_labels_batch}\n",
    "            \n",
    "            test_prediction_ = session.run(test_prediction, feed_dict=feed_dict2)\n",
    "            test_accuracy = accuracy(test_prediction_, test_labels_batch)\n",
    "            test_accuracies.append(test_accuracy)\n",
    "            test_preds.append(test_prediction_)\n",
    "            \n",
    "            valid_dataset_batch = valid_dataset[offset2:(offset2 + batch_size), :, :]\n",
    "            valid_labels_batch = valid_labels[offset2:(offset2 + batch_size), :]\n",
    "            feed_dict3 = {tf_valid_dataset : valid_dataset_batch, tf_valid_labels : valid_labels_batch}\n",
    "            \n",
    "            valid_prediction_ = session.run(valid_prediction, feed_dict=feed_dict3)\n",
    "            valid_accuracy = accuracy(valid_prediction_, valid_labels_batch)\n",
    "            valid_accuracies.append(valid_accuracy)\n",
    "            valid_preds.append(valid_prediction_)\n",
    " \n",
    "            message = \"step {:04d} : loss is {:06.2f}, accuracy on training set {:02.2f} %, accuracy on test set {:02.2f} accuracy on valid set {:02.2f} %\".format(step, l, train_accuracy, test_accuracy, valid_accuracy)\n",
    "            print(message)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test metrics\n",
    "fpr_test, tpr_test, thresholds_test = roc_curve((np.argmax(test_preds,1)), (np.argmax(test_labels,1)))\n",
    "auc_test = auc(fpr_test, tpr_test)\n",
    "\n",
    "plt.figure(1)\n",
    "plt.plot([0, 1], [0, 1], 'k--')\n",
    "plt.plot(fpr_test, tpr_test, label='Alexnet (area = {:.3f})'.format(auc_test))\n",
    "plt.xlabel('False positive rate')\n",
    "plt.ylabel('True positive rate')\n",
    "plt.title('ROC curve')\n",
    "plt.legend(loc='best')\n",
    "plt.show()\n",
    "# Zoom in view of the upper left corner.\n",
    "plt.figure(2)\n",
    "plt.xlim(0, 0.2)\n",
    "plt.ylim(0.8, 1)\n",
    "plt.plot([0, 1], [0, 1], 'k--')\n",
    "plt.plot(fpr_test, tpr_test, label='Alexnet (area = {:.3f})'.format(auc_test))\n",
    "plt.xlabel('False positive rate')\n",
    "plt.ylabel('True positive rate')\n",
    "plt.title('ROC curve (zoomed in at top left)')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizing Accuracies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def average_points(points, stepsize = 10):\n",
    "    averaged_points = []\n",
    "    for ii in range(stepsize,len(points),stepsize):\n",
    "        subsection  = points[ii-stepsize:ii]\n",
    "        average = np.nanmean(subsection)\n",
    "        averaged_points.append(average)\n",
    "    return averaged_points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_steps = 10\n",
    "ylimit = [0,100]\n",
    "labels = ['Train accuracy', 'Test accuracy', 'Validation accuracy']\n",
    "ylabel = \"Accuracy [%]\"\n",
    "xlabel = \"Number of Iterations\"\n",
    "title = \"Accuracy of road detection in Aerial Images\"\n",
    "colors = ['r', 'g', 'b']\n",
    "\n",
    "list_accuracies = [train_accuracies, test_accuracies, valid_accuracies]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "ax.set_ylim(ylimit)\n",
    "ax.set_ylabel(ylabel, fontsize=16)\n",
    "ax.set_xlabel(xlabel, fontsize=16)\n",
    "ax.set_title(title, fontsize=20)\n",
    "\n",
    "\n",
    "for ii, accuracies in enumerate(list_accuracies):\n",
    "    color = colors[ii]\n",
    "    label = labels[ii]\n",
    "    if ii > 0:\n",
    "        y_values = accuracies\n",
    "        x_values = range(0,num_steps, 10)\n",
    "        ax.plot(x_values, y_values, '.-{}'.format(color), label = label)\n",
    "    else:\n",
    "        y_values_ = accuracies\n",
    "        y_values = average_points(y_values_, 5)\n",
    "        x_values = range(1,len(y_values_),5)\n",
    "        ax.plot(x_values, y_values, '.{}'.format(color), label = label)\n",
    "ax.legend(loc='lower right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}