{
  "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.6.0"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0,
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "bc5d8a40-770f-f003-a46c-9f2296b6249c",
        "_active": false
      },
      "source": "# NDVI for vegetation detection",
      "outputs": [],
      "execution_count": null,
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "3cd94403-ccb7-7249-9737-e53091735a11",
        "_active": false
      },
      "outputs": [],
      "source": "# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nfrom subprocess import check_output\nprint(check_output([\"ls\", \"../input\"]).decode(\"utf8\"))\n\n# Any results you write to the current directory are saved as output.",
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "41f183fb-ec11-6edc-24f4-b941d94bc215",
        "_active": false,
        "collapsed": false
      },
      "outputs": [],
      "source": "from collections import defaultdict\nimport cv2\nfrom shapely.geometry import MultiPolygon, Polygon\nimport shapely.wkt\nimport shapely.affinity\nimport numpy as np\nimport tifffile as tiff\n#import gdal\nfrom osgeo import gdal\n#import gdalconst\n\nimport matplotlib.pylab as plt\n%matplotlib inline",
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "2b4e9639-927b-91d8-54e5-8067794ac678",
        "_active": true
      },
      "outputs": [],
      "source": "gdal.AllRegister()\ngdal.VersionInfo()",
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "47ca4174-774a-5fda-f429-31bc7a534013",
        "_active": false
      },
      "outputs": [],
      "source": "DATA_3_BANDS='../input/three_band/'\nDATA_16_BANDS='../input/sixteen_band/'",
      "execution_state": "idle"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "36e21de1-acb5-b6f5-dfa8-cc0a86d764f9",
        "_active": false
      },
      "source": "Visualize all available bands",
      "outputs": [],
      "execution_count": null,
      "execution_state": "idle"
    },
    {
      "metadata": {
        "_cell_guid": "f12f4907-5ec3-4eab-e98e-67c8765da5a4",
        "_active": false,
        "collapsed": false
      },
      "source": null,
      "execution_count": null,
      "cell_type": "code",
      "outputs": [],
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "622a5763-fecd-1a10-d45f-3d3705538923",
        "_active": false
      },
      "outputs": [],
      "source": "image_id = \"6020_2_2\"\n\nfname_3b = DATA_3_BANDS + \"/{}.tif\".format(image_id)\nfname_pan = DATA_16_BANDS + \"/{}_P.tif\".format(image_id)\nfname_ms = DATA_16_BANDS + \"/{}_M.tif\".format(image_id)\nfname_swir = DATA_16_BANDS + \"/{}_A.tif\".format(image_id)\n\nimg_3b = gdal.Open(fname_3b, gdalconst.GA_ReadOnly)\nassert img_3b, \"WTF\"\nprint(\"Image 3 bands info :\", img_3b.RasterXSize, img_3b.RasterYSize,img_3b.RasterCount)\nimg_pan = gdal.Open(fname_pan, gdalconst.GA_ReadOnly)\nassert img_pan, \"WTF\"\nprint(\"Image Pan info :\", img_pan.RasterXSize, img_pan.RasterYSize,img_pan.RasterCount)\nimg_ms = gdal.Open(fname_ms, gdalconst.GA_ReadOnly)\nassert img_ms, \"WTF\"\nprint(\"Image MS info :\", img_ms.RasterXSize, img_ms.RasterYSize,img_ms.RasterCount)\nimg_swir = gdal.Open(fname_swir, gdalconst.GA_ReadOnly)\nassert img_swir, \"WTF\"\nprint(\"Image SWIR info :\", img_swir.RasterXSize, img_swir.RasterYSize,img_swir.RasterCount)",
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "124ecc1a-e533-c935-2497-e6b7268ad7f5",
        "_active": false
      },
      "outputs": [],
      "source": "def scale_percentile(matrix):\n    if len(matrix.shape) == 2:\n        matrix = matrix.reshape(matrix.shape + (1,))\n    w, h, d = matrix.shape\n    matrix = np.reshape(matrix, [w * h, d]).astype(np.float64)\n    # Get 2nd and 98th percentile\n    mins = np.percentile(matrix, 1, axis=0)\n    maxs = np.percentile(matrix, 99, axis=0) - mins\n    matrix = 255*(matrix - mins[None, :]) / maxs[None, :]\n    matrix = np.reshape(matrix, [w, h, d] if d > 1 else [w, h])\n    matrix = matrix.clip(0, 255).astype(np.uint8)\n    return matrix   \n    \n\ndef display_img_1b(img_1b, roi=None):   \n    img_1b_data = img_1b.ReadAsArray()\n    if roi is not None:\n        y,yh,x,xw = roi\n        img_1b_data = img_1b_data[y:yh,x:xw]\n    plt.figure(figsize=(8,4))\n    plt.imshow(scale_percentile(img_1b_data), cmap='gray')\n\n    \ndef display_img_3b(img_3b, roi=None):\n    img_3b_data = img_3b.ReadAsArray()\n    if roi is not None:\n        y,yh,x,xw = roi\n        img_3b_data = img_3b_data[:,y:yh,x:xw]\n    plt.figure(figsize=(8,4))\n    for i in [0,1,2]:\n        plt.subplot(1,3,i+1)\n        plt.imshow(scale_percentile(img_3b_data[i,:,:]), cmap='gray')\n        plt.title(\"Channel %i\" % i)\n    \ndef display_img_8b(img_ms, roi=None):\n    img_ms_data = img_ms.ReadAsArray()\n    if roi is not None:\n        y,yh,x,xw = roi\n        img_ms_data = img_ms_data[:,y:yh,x:xw]\n    plt.figure(figsize=(8,4))\n    for i in range(8):\n        plt.subplot(2,4,i+1)\n        plt.imshow(scale_percentile(img_ms_data[i,:,:]), cmap='gray')\n        plt.title(\"Channel %i\" % i)",
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "4b86b451-bc12-2937-19ea-de8159c83b3d",
        "_active": false
      },
      "outputs": [],
      "source": "from ipywidgets import interact, IntSlider\n\nfig = None\ndef interactive_visu(roi_x, roi_w, roi_y, roi_h):\n    roi = [roi_x,roi_x+roi_w,roi_y,roi_y+roi_h]\n    display_img_1b(img_pan, roi)\n    plt.suptitle(\"Pansharpened image ROI\")\n    display_img_3b(img_3b, roi)\n    plt.suptitle(\"3 bands image ROI\")\n    display_img_8b(img_ms, [int(r/4) for r in roi])\n    plt.suptitle(\"8 bands image ROI\")\n    display_img_8b(img_swir, [int(r/4/6.2) for r in roi])\n    _ = plt.suptitle(\"8 bands SWIR image ROI\")\n\n    \n_ = interact(interactive_visu,           \n         roi_x=IntSlider(value=100, min=0, max=3500, continuous_update=False), \n         roi_w=FloatSlider(value=200, min=150, max=350, continuous_update=False), \n         roi_y=FloatSlider(value=0, min=0, max=3500, continuous_update=False), \n         roi_h=FloatSlider(value=200, min=150, max=350, continuous_update=False))",
      "execution_state": "idle"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "1b2c3dae-e5e7-774f-351d-6d2643f8322a",
        "_active": false
      },
      "source": "Visualize 3-bands data ",
      "outputs": [],
      "execution_count": null,
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "ed5b8b9a-c37d-9e98-ddd9-a66764e13125",
        "_active": false
      },
      "outputs": [],
      "source": "def render_rgb(data_3_channels, bgr2rgb=False):    \n    out = np.zeros_like(data_3_channels).astype(np.uint8)\n    in_channels = [0,1,2]\n    out_channels = [2,1,0] if bgr2rgb else in_channels\n    for c, nc in zip(in_channels, out_channels):\n        band = data_3_channels[:,:,c]\n        min_value = np.percentile(band, 3)\n        max_value = np.percentile(band, 97)     \n        band[band < min_value] = min_value\n        band[band > max_value] = max_value\n        out[:,:,nc] = 255 * (band - min_value)/(max_value - min_value)    \n    return out",
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "82d15755-3435-717b-e842-afbd1b8be8c2",
        "_active": false
      },
      "outputs": [],
      "source": "!ls ../input/three_band/ | head",
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "0dafa663-ca1c-3262-2c29-6b2a5da53857",
        "_active": false
      },
      "outputs": [],
      "source": "rgb_image_ids = [\"6010_0_0\", \"6090_4_2\", \"6170_3_0\", \"6140_4_2\", \"6120_2_2\"]\n#rgb_image_ids = [\"6120_2_2\"]\n\nfor rgb_image_id in rgb_image_ids:\n    \n    test_rgb_image_filename = \"../input/three_band/{}.tif\".format(rgb_image_id)\n    test_rgb_image = gdal.Open(test_rgb_image_filename)\n    assert test_rgb_image is not None, \"WTF\"\n    print(\"File : \", test_rgb_image_filename)\n    print(\"Metadata list: \", test_rgb_image.GetMetadata_List())\n    print(\"Metadata domaines: \", test_rgb_image.GetMetadataDomainList())\n    print(\"Projection reference: \", test_rgb_image.GetProjectionRef())\n    print(\"Geotransform: \", test_rgb_image.GetGeoTransform())\n    band1 = test_rgb_image.GetRasterBand(1)\n    print(\"Pixel depth: \", gdal.GetDataTypeName(band1.DataType))\n    test_rgb_image_data = test_rgb_image.ReadAsArray().transpose([1,2,0])\n    test_rgb_image_data = render_rgb(test_rgb_image_data)\n    \n    plt.figure(figsize=(10,4))\n    plt.subplot(121)\n    plt.imshow(test_rgb_image_data)\n    plt.subplot(122)\n    plt.imshow(test_rgb_image_data[2900:3200,2000:2300,:])",
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "fdbd8583-e824-d818-1d9d-7f7581b20301",
        "_active": false
      },
      "outputs": [],
      "source": "!ls ../input/sixteen_band/ | head",
      "execution_state": "idle"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "87ddc0c0-b4b1-4a5f-1ef0-795cd23b3655",
        "_active": false
      },
      "source": "Visualize RGB channels from 16-bands images\nMultispectral image has `[6; 2; 3; 7; 4; 8; 5; 9; ]` WorldView-3 specified channels, which are \n\n- 1    Panchromatic – 450 to 800\n- 2    Coastal – 400 to 450\n- 3     Blue – 450 to 510\n- 4    Green – 510 to 580\n- 5    Yellow – 585 to 625\n- 6    Red – 630 to 690\n- 7    Red edge – 705 to 745\n- 8    NIR-1 – 770 to 895\n- 9    NIR-2 – 860 to 1040\n\nTo get RGB channels we should take `[6, 4, 3]` channels which correspond to `[0, 4, 2]` indices",
      "outputs": [],
      "execution_count": null,
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "c3c834e7-f1b9-5206-54da-cdb28d16edd9",
        "_active": false
      },
      "outputs": [],
      "source": "#mb_image_ids = [\"6010_0_0_M\", \"6090_4_2_M\", \"6170_3_0_M\", \"6140_4_2_M\", \"6120_2_2_M\"]\nmb_image_ids = [\"6120_2_2_M\"]\n\nrgb_channels = (0, 4, 2)\nfor mb_image_id in mb_image_ids:\n    \n    test_mb_image_filename = \"../input/sixteen_band/{}.tif\".format(mb_image_id)\n    test_mb_image = gdal.Open(test_mb_image_filename)\n    assert test_mb_image is not None, \"WTF\"\n    print(\"File : \", test_mb_image_filename)\n    print(\"Metadata list: \", test_mb_image.GetMetadata_List())\n    print(\"Metadata domaines: \", test_mb_image.GetMetadataDomainList())\n    print(\"Projection reference: \", test_mb_image.GetProjectionRef())\n    print(\"Geotransform: \", test_mb_image.GetGeoTransform())\n    band1 = test_mb_image.GetRasterBand(1)\n    print(\"Pixel depth: \", gdal.GetDataTypeName(band1.DataType))\n\n    test_mb_image_data = test_mb_image.ReadAsArray()  \n    test_rgb_image_data = test_mb_image_data[rgb_channels,:,:].transpose([1,2,0])\n    test_rgb_image_data = render_rgb(test_rgb_image_data)\n    \n    plt.figure(figsize=(10,4))\n    plt.subplot(121)\n    plt.imshow(test_rgb_image_data)\n    plt.subplot(122)\n    plt.imshow(test_rgb_image_data[200:300,200:300,:])",
      "execution_state": "idle"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "c382cff0-199f-bb33-6818-206438098c63",
        "_active": false
      },
      "source": "Compute NDVI ",
      "outputs": [],
      "execution_count": null,
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "b0988567-910b-37a7-bb48-17d499d3a287",
        "_active": false
      },
      "outputs": [],
      "source": "mb_image_ids = [\"6120_2_2_M\"]\n\nfor mb_image_id in mb_image_ids:\n    \n    test_mb_image_filename = \"../input/sixteen_band/{}.tif\".format(mb_image_id)\n    test_mb_image = gdal.Open(test_mb_image_filename)\n    assert test_mb_image is not None, \"WTF\"\n    print(\"File : \", test_mb_image_filename)\n    print(\"Metadata list: \", test_mb_image.GetMetadata_List())\n    print(\"Metadata domaines: \", test_mb_image.GetMetadataDomainList())\n    print(\"Projection reference: \", test_mb_image.GetProjectionRef())\n    print(\"Geotransform: \", test_mb_image.GetGeoTransform())\n    band1 = test_mb_image.GetRasterBand(1)\n    print(\"Pixel depth: \", gdal.GetDataTypeName(band1.DataType))\n\n    test_mb_image_data = test_mb_image.ReadAsArray()  \n    nir = test_mb_image_data[7,:,:]\n    vis = test_mb_image_data[6,:,:]\n    ndvi = (nir - vis) / (nir + vis)\n        \n    plt.figure(figsize=(10,4))\n    plt.subplot(121)\n    plt.imshow(ndvi)\n    plt.subplot(122)\n    plt.imshow(ndvi[200:300,200:300])",
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "9bd2dcfb-4dae-54e1-06d9-17ef412ccb91",
        "_active": false
      },
      "outputs": [],
      "source": "labels = [\n    None, \n    # 1\n    \"Buildings - large building, residential, non-residential, fuel storage facility, fortified building\",\n    # 2\n    \"Misc. Manmade structures\", \n    # 3\n    \"Road\", \n    # 4\n    \"Track - poor/dirt/cart track, footpath/trail\",\n    # 5\n    \"Trees - woodland, hedgerows, groups of trees, standalone trees\",\n    # 6\n    \"Crops - contour ploughing/cropland, grain (wheat) crops, row (potatoes, turnips) crops\",\n    # 7\n    \"Waterway\", \n    # 8\n    \"Standing water\",\n    # 9\n    \"Vehicle Large - large vehicle (e.g. lorry, truck,bus), logistics vehicle\",\n    # 10\n    \"Vehicle Small - small vehicle (car, van), motorbike\",    \n]",
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "6431a35a-4312-e942-951b-3b80e3e986ce",
        "_active": false
      },
      "outputs": [],
      "source": "train_wkt = pd.read_csv(\"../input/train_wkt_v4.csv\")\ntrain_wkt.head()",
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "8d929f7c-1d8a-37ab-8a5a-b5994a0041ce",
        "_active": false
      },
      "outputs": [],
      "source": "train_wkt[train_wkt[\"ImageId\"] == \"6120_2_2\"]",
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "982ecf2e-b051-1a50-cf22-e4fe365ce0ce",
        "_active": false
      },
      "outputs": [],
      "source": "GRID_SIZE = pd.read_csv('../input/grid_sizes.csv', names=['ImageId', 'Xmax', 'Ymin'], skiprows=1)\nGRID_SIZE.columns = ['ImageId','Xmax','Ymin']\nTRAIN_WKT = pd.read_csv('../input/train_wkt_v4.csv')\n\ndef get_grid_size(image_id):\n    x_max = GRID_SIZE[GRID_SIZE['ImageId']==image_id].Xmax.values[0]\n    y_min = GRID_SIZE[GRID_SIZE['ImageId']==image_id].Ymin.values[0]\n    return x_max, y_min\n\ndef get_scalers(image_shape, x_max, y_min):\n    h, w = image_shape  # they are flipped so that mask_for_polygons works correctly\n    w_ = w * (w / (w + 1))\n    h_ = h * (h / (h + 1))\n    return w_ / x_max, h_ / y_min\n\ndef generate_image_mask(image_id, class_type):\n    data_mask = (TRAIN_WKT[\"ImageId\"] == image_id) & (TRAIN_WKT[\"ClassType\"] == class_type)\n    poly = TRAIN_WKT[data_mask][\"MultipolygonWKT\"].values[0]\n    #print(\"poly=\", poly)\n    train_polygons = shapely.wkt.loads(poly)\n    \n    \n    return train_polygons\n\n\npoly = generate_image_mask(\"6120_2_2\", 5, train_wkt)\npoly",
      "execution_state": "idle"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "2c358f1d-2043-ffc5-0382-3a30c9f3f9f2",
        "_active": false
      },
      "outputs": [],
      "source": null,
      "execution_state": "idle"
    }
  ]
}