{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "14190218-4ba1-e5af-2acd-7d1e5932cd1b"
      },
      "source": [
        "\n",
        "## End-to-end object-based solution for DSTL\n",
        "\n",
        "There are typically two approaches for geo-image-segmentation: pixel-based and object-based and I'm just surprised that the latter was rarely mentioned in the forum and there seems no kernels available for this competition so I decided to share my object-based solution and hopefully it will be of help.\n",
        "\n",
        "Again, this competition comes with tons of challenges mostly in programmming/engineering which have been a pain for me. I'm grateful for those who shared their scripts/solutions. Without them, I weren't be able to make a single valid submission.\n",
        "\n",
        "\n",
        "\n",
        "This solution was inspried by the following two articles:\n",
        "\n",
        "* Python for Object Based Image Analysis (OBIA)\n",
        "https://www.machinalis.com/blog/obia/\n",
        "\n",
        "* A Python-Based Open Source System for Geographic Object-Based Image Analysis (GEOBIA) Utilizing Raster Attribute Tables\n",
        "http://www.mdpi.com/2072-4292/6/7/6111/htm\n",
        "\n",
        "Many ideas/functions/tools were borrowed from Konstantin Lopuhin's great kernel:\n",
        "https://www.kaggle.com/lopuhin/dstl-satellite-imagery-feature-detection/full-pipeline-demo-poly-pixels-ml-poly\n",
        " "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "39d981f2-eb3b-84e3-3b74-887646521bdb"
      },
      "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",
        "\n",
        "import numpy as np # linear algebra\n",
        "import 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",
        "\n",
        "from subprocess import check_output\n",
        "print(check_output([\"ls\", \"../input\"]).decode(\"utf8\"))\n",
        "\n",
        "# Any results you write to the current directory are saved as output."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "12e5041e-da88-077e-4c6b-23a53cb729b7"
      },
      "source": [
        "Prepping"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "a507651c-c342-57fa-2627-1ee6e9d3ae54"
      },
      "outputs": [],
      "source": [
        "import pandas as pd\n",
        "\n",
        "GRID_SIZE = pd.read_csv('../input/grid_sizes.csv', names=['ImageId', 'Xmax', 'Ymin'], skiprows=1)\n",
        "GRID_SIZE.columns = ['ImageId','Xmax','Ymin']\n",
        "TRAIN_WKT = pd.read_csv('../input/train_wkt_v4.csv')\n",
        "\n",
        "## Exclude empty polygons\n",
        "# TRAIN_WKT[TRAIN_WKT['MultipolygonWKT']!='MULTIPOLYGON EMPTY']\n",
        "\n",
        "CLASSES = {\n",
        "        1 : 'Buildings',\n",
        "        2 : 'Misc',\n",
        "        3 : 'Road',\n",
        "        4 : 'Track',\n",
        "        5 : 'Trees',\n",
        "        6 : 'Crops',\n",
        "        7 : 'Waterway',\n",
        "        8 : 'Standing water',\n",
        "        9 : 'Vehicle Large',\n",
        "        10 : 'Vehicle Small',\n",
        "        }\n",
        "\n",
        "COLORS = {\n",
        "        1 : '0.7',\n",
        "        2 : '0.4',\n",
        "        3 : '#b35806',\n",
        "        4 : '#dfc27d',\n",
        "        5 : '#1b7837',\n",
        "        6 : '#a6dba0',\n",
        "        7 : '#74add1',\n",
        "        8 : '#4575b4',\n",
        "        9 : '#f46d43',\n",
        "        10: '#d73027',\n",
        "        }\n",
        "ZORDER = {\n",
        "        1 : 5,\n",
        "        2 : 5,\n",
        "        3 : 4,\n",
        "        4 : 1,\n",
        "        5 : 3,\n",
        "        6 : 2,\n",
        "        7 : 7,\n",
        "        8 : 8,\n",
        "        9 : 9,\n",
        "        10: 10,\n",
        "        }    \n",
        "\n",
        "sample_submission = pd.read_csv('../input/sample_submission.csv')\n",
        "test_image_ids = sample_submission.ImageId.unique()\n",
        "train_image_ids = TRAIN_WKT.ImageId.unique()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "a1a5d285-1fc7-8866-6deb-273ab12347db"
      },
      "source": [
        "### Run following command to create working folders if needed"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "69b44c19-0448-c03f-5ff3-4d938e01d180"
      },
      "outputs": [],
      "source": [
        "mkdir ../input/feature;mkdir ../input/feature_train;mkdir ../input/feature/segment"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "188d5713-4bd5-4792-d6db-1afe83358abe"
      },
      "source": [
        "## Image segmentation\n",
        "\n",
        "* We will use RSGISLib for image segmentation\n",
        "\n",
        "    * Webiste: http://www.rsgislib.org/index.html\n",
        "    * Installation:\n",
        "    \n",
        "        * http://www.rsgislib.org/download.html#binary-downloads \n",
        "    \n",
        "        * https://groups.google.com/forum/#!searchin/rsgislib-support/MAC$20INSTALLATION%7Csort:relevance/rsgislib-support/OqnN9y--ff0/R-Hevkw2BAAJ\n",
        "\n",
        "\n",
        "**IMPORTANT NOTE:** RSGISlib installation will create a new Python environment osgeoenv and you'll need to reinstall necessary packages in this environment if you want to use it to continue subsquent works. Other wise, once image segmentation is done, you can switch back to your normal environment to continue following works."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "b6701dc5-7b5c-f768-fd39-160e81b3b6c4"
      },
      "outputs": [],
      "source": [
        "from rsgislib.segmentation import segutils\n",
        "import time\n",
        "def run_image_segmentation(input_image, output_image, num_clusters=60, min_pixels=100, dist_thres=100):\n",
        "    '''\n",
        "    This funcation will segment input image and create a mask file where the pixel value is segment id\n",
        "    '''\n",
        "    print ('Running segmentation for image \"%s\" ...' % (input_image))    \n",
        "    start = time.time()\n",
        "    segutils.runShepherdSegmentation(input_image, output_image, \n",
        "                                     gdalformat='GTiff',\n",
        "                                     numClusters=num_clusters, \n",
        "                                     minPxls=min_pixels, \n",
        "                                     noStats=True, ## have to setup this param otherwise there will be an error\n",
        "                                     distThres=dist_thres,\n",
        "                                     processInMem=True)\n",
        "    print ('Segmentation for image \"%s\" finished in %d seconds.' % (input_image, time.time()-start))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "021c8c1f-4b59-afcd-a2f0-8ee7e99f5360"
      },
      "outputs": [],
      "source": [
        "# imageIds in a DataFrame\n",
        "all_image_ids = GRID_SIZE.ImageId.unique()\n",
        "\n",
        "for image_id in all_image_ids:\n",
        "    input_image = '../input/sixteen_band/'+image_id+'_M.tif'\n",
        "    output_image = '../input/segment/'+image_id+'_M_SEG.tif'\n",
        "    run_image_segmentation(input_image, output_image)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "82afec72-b351-8cf7-c0fd-eea6d22c0e6e"
      },
      "source": [
        "## Feature extraction"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "d24cca2d-a16e-ced7-5ef2-d97e470c2926"
      },
      "outputs": [],
      "source": [
        "import tifffile as tiff\n",
        "import numpy as np\n",
        "import warnings\n",
        "from scipy import stats as sc_stats\n",
        "\n",
        "def get_poly_features(poly_pixels):\n",
        "    \"\"\"For each band, compute: min, max, mean, variance, skewness, kurtosis\"\"\"\n",
        "    features = []\n",
        "    if len(poly_pixels.shape)<2:\n",
        "        return None\n",
        "    n_pixels, n_bands = poly_pixels.shape\n",
        "    for b in range(n_bands):\n",
        "        stats = sc_stats.describe(poly_pixels[:,b])\n",
        "        band_stats = list(stats.minmax) + list(stats)[2:]\n",
        "        if n_pixels == 1:\n",
        "            # scipy.stats.describe raises a Warning and sets variance to nan\n",
        "            band_stats[3] = 0.0  # Replace nan with something (zero)\n",
        "        features += band_stats\n",
        "    return features\n",
        "\n",
        "def segment_to_features(segment_mask, image):\n",
        "    '''\n",
        "    Extract features from polygons generated by segmentation.\n",
        "    '''\n",
        "    image_shape = image.shape[:2]\n",
        "    poly_ids = np.unique(segment_mask)\n",
        "    with warnings.catch_warnings():\n",
        "        warnings.simplefilter(\"ignore\")\n",
        "\n",
        "        features = []\n",
        "        for poly_id in poly_ids:\n",
        "            poly_pixels = image[segment_mask==poly_id]\n",
        "            poly_features = get_poly_features(poly_pixels)\n",
        "            features.append(poly_features)     \n",
        "    return features, poly_ids \n",
        "\n",
        "def image_to_array(image_id,image_type='M'):\n",
        "    '''\n",
        "    image to array\n",
        "    '''\n",
        "    image_type ='M'\n",
        "    if image_type =='3':\n",
        "        image = tiff.imread('../input/three_band/{}.tif'.format(image_id)).transpose([1, 2, 0])\n",
        "    elif image_type =='M':\n",
        "        image = tiff.imread('../input/sixteen_band/{}_M.tif'.format(image_id)).transpose([1, 2, 0])\n",
        "    elif image_type =='A':\n",
        "        image = tiff.imread('../input/sixteen_band/{}_A.tif'.format(image_id)).transpose([1, 2, 0])\n",
        "    elif image_type =='P':\n",
        "        image = tiff.imread('../input/sixteen_band/{}_P.tif'.format(image_id))\n",
        "#     image = exposure.rescale_intensity(image)\n",
        "    return image\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "61059b88-820d-49b8-a373-bbb610acb2b9"
      },
      "source": [
        "## Extract features for testing images"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "919b2cde-a8f0-0724-98c4-7281a43eacbd"
      },
      "outputs": [],
      "source": [
        "## Feature extraction for testing images\n",
        "for image_id in test_image_ids:\n",
        "    image = image_to_array(image_id,image_type='M')\n",
        "    segment_mask = tiff.imread('../input/segment/'+image_id+'_M_SEG.tif')\n",
        "    start = time.time()\n",
        "\n",
        "    features, segment_ids = segment_to_features(segment_mask, image)\n",
        "    feature_df = pd.DataFrame(features\n",
        "                ,columns=['b1_min','b1_max','b1_mean','b1_variance','b1_skewness','b1_kurtosis',\n",
        "                          'b2_min','b2_max','b2_mean','b2_variance','b2_skewness','b2_kurtosis',\n",
        "                          'b3_min','b3_max','b3_mean','b3_variance','b3_skewness','b3_kurtosis',\n",
        "                          'b4_min','b4_max','b4_mean','b4_variance','b4_skewness','b4_kurtosis',\n",
        "                          'b5_min','b5_max','b5_mean','b5_variance','b5_skewness','b5_kurtosis',\n",
        "                          'b6_min','b6_max','b6_mean','b6_variance','b6_skewness','b6_kurtosis',\n",
        "                          'b7_min','b7_max','b7_mean','b7_variance','b7_skewness','b7_kurtosis',\n",
        "                          'b8_min','b8_max','b8_mean','b8_variance','b8_skewness','b8_kurtosis'\n",
        "                         ])\n",
        "    feature_df['segment_id'] = segment_ids\n",
        "    feature_df['image_id'] = image_id\n",
        "    feature_df.to_csv('../input/feature/'+image_id+'.csv',index = False)\n",
        "    print ('Feature extraction for image %s finished in %d seconds' % (image_id,time.time() - start))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "63c1a6b7-5998-9604-f5de-940b7ba06dda"
      },
      "outputs": [],
      "source": [
        "## Extract features for training images\n",
        "### Run following commands if required pacakges have not been installed in osgoenv environment\n",
        "\n",
        "!source activate osgeoenv;pip install shapely;pip install opencv-python;pip install rasterio"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "9764ded8-4982-7d73-56ad-23799f461549"
      },
      "outputs": [],
      "source": [
        "## !source activate osgeoenv;pip install shapely\n",
        "import shapely.wkt\n",
        "import shapely.affinity\n",
        "from shapely.geometry import MultiPolygon, Polygon\n",
        "import cv2\n",
        "from collections import Counter\n",
        "# make sure rasterio was imported after shapely otherwise it may cause kernel error!!!\n",
        "import rasterio.features\n",
        "\n",
        "def get_grid_size(image_id):\n",
        "    '''\n",
        "    '''\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",
        "\n",
        "def get_scalers(image_shape, x_max, y_min):\n",
        "    '''\n",
        "    To provide scalers that will be used to scale predicted polygons\n",
        "    '''\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",
        "\n",
        "def image_to_features(image_id,image_type='M', class_type = 0):\n",
        "    image = image_to_array(image_id)\n",
        "    image_shape = image.shape[:2]\n",
        "    x_max, y_min = get_grid_size(image_id)\n",
        "    x_scaler, y_scaler = get_scalers(image_shape, x_max, y_min)\n",
        "\n",
        "    start = time.time()\n",
        "    features = np.zeros((0,48)) ## 48 is the number of features from scipy.stats (6) * number of bands (8)\n",
        "    labels = np.array([])\n",
        "    poly_ids = np.array([])\n",
        "    mask = np.zeros(image_shape)\n",
        "    if class_type==0: ##all classes\n",
        "        for cls in CLASSES:\n",
        "            train_polygons = TRAIN_WKT[(TRAIN_WKT['ImageId']==image_id) & \n",
        "                                       (TRAIN_WKT['ClassType']==cls)].MultipolygonWKT.values[0]\n",
        "            if train_polygons != 'MULTIPOLYGON EMPTY':\n",
        "                ## Check if polygons is empty\n",
        "                train_polygons = shapely.wkt.loads(train_polygons)\n",
        "                ## Scale polygons based on image grid size\n",
        "                train_polygons = shapely.affinity.scale(train_polygons,\n",
        "                                                        xfact=x_scaler,\n",
        "                                                        yfact=y_scaler,\n",
        "                                                        origin=(0, 0, 0))\n",
        "                poly_features, pids, poly_mask = poly_to_features(train_polygons, image, cls)\n",
        "                features = np.vstack((features, poly_features))\n",
        "                labels = np.hstack((labels, np.full((len(poly_features)),cls, dtype=np.int)))\n",
        "                poly_ids = np.hstack((poly_ids, pids))\n",
        "                mask = np.max(np.stack((mask,poly_mask)),axis=0)\n",
        "    print (\"Feature extracted in %d seconds\" % (time.time()-start) )        \n",
        "    return features,labels,poly_ids,mask\n",
        "\n",
        "def mask_for_polygons(polygons, image_size):\n",
        "    image_mask = np.zeros(image_size, np.uint8)\n",
        "    if not polygons:\n",
        "        return image_mask\n",
        "    int_coords = lambda x: np.array(x).round().astype(np.int32)\n",
        "    exteriors = [int_coords(poly.exterior.coords) for poly in polygons]\n",
        "    interiors = [int_coords(pi.coords) for poly in polygons\n",
        "                 for pi in poly.interiors]\n",
        "    cv2.fillPoly(image_mask, exteriors, 1)\n",
        "    cv2.fillPoly(image_mask, interiors, 0)\n",
        "    return image_mask\n",
        "\n",
        "def poly_to_mask(polygons, image_shape, class_type):\n",
        "    mask = np.zeros(image_shape, np.uint8)    \n",
        "    if not polygons:\n",
        "        return image\n",
        "    int_coords = lambda x: np.array(x).round().astype(np.int32)\n",
        "\n",
        "    exteriors = []\n",
        "    interiors = []\n",
        "    for pid, poly in enumerate(polygons):\n",
        "        poly_mask = np.zeros(image_shape, np.uint8)\n",
        "        exteriors=[int_coords(poly.exterior.coords)]\n",
        "        interiors = []\n",
        "        for pi in poly.interiors:\n",
        "            interiors.append(int_coords(pi.coords) )\n",
        "\n",
        "        cv2.fillPoly(poly_mask, exteriors, 1)\n",
        "\n",
        "        cv2.fillPoly(poly_mask, interiors, 0)\n",
        "        poly_id = (pid + 1) * 100 + class_type\n",
        "        poly_mask = poly_mask * poly_id\n",
        "        mask = np.max(np.stack((mask, poly_mask)),axis=0)\n",
        "    return mask \n",
        "\n",
        "def poly_to_features(polygons, image, class_type):\n",
        "    '''\n",
        "    Extract features for training image. \n",
        "    Polygons are extracted from training WKT then converted to masks.\n",
        "    Featuers are statistical metrics of each polygon\n",
        "    '''\n",
        "    image_shape = image.shape[:2]\n",
        "    poly_mask = poly_to_mask(polygons, image_shape, class_type)\n",
        "    poly_ids = np.unique(poly_mask)\n",
        "    poly_ids = poly_ids[poly_ids != 0]\n",
        "    with warnings.catch_warnings():\n",
        "        warnings.simplefilter(\"ignore\")\n",
        "\n",
        "        features = []\n",
        "        for poly_id in poly_ids:\n",
        "            poly_pixels = image[poly_mask==poly_id]\n",
        "            poly_features = get_poly_features(poly_pixels)\n",
        "            features.append(poly_features)     \n",
        "    return features, poly_ids, poly_mask \n",
        "\n",
        "def image_to_train(image_id, class_type, image_type='3'):\n",
        "    # Get grid size: x_max and y_min\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",
        "\n",
        "    # Load train poly with shapely\n",
        "    train_polygons = shapely.wkt.loads(TRAIN_WKT[(TRAIN_WKT['ImageId']==image_id) & \n",
        "                                                (TRAIN_WKT['ClassType']==class_type)].MultipolygonWKT.values[0])\n",
        "\n",
        "    # Read image with tiff\n",
        "    if image_type =='3':\n",
        "        image = tiff.imread('../input/three_band/{}.tif'.format(image_id)).transpose([1, 2, 0])\n",
        "    if image_type =='M':\n",
        "        image = tiff.imread('../input/sixteen_band/{}_M.tif'.format(image_id)).transpose([1, 2, 0])\n",
        "    if image_type =='A':\n",
        "        image = tiff.imread('../input/sixteen_band/{}_A.tif'.format(image_id)).transpose([1, 2, 0])\n",
        "    if image_type =='P':\n",
        "        image = tiff.imread('../input/sixteen_band/{}_P.tif'.format(image_id))\n",
        "    image_size = image.shape[:2]\n",
        "    x_scaler, y_scaler = get_scalers(image_size, x_max, y_min)\n",
        "\n",
        "    # Scale polygons\n",
        "    train_polygons_scaled = shapely.affinity.scale(train_polygons,\n",
        "                                                   xfact=x_scaler,\n",
        "                                                   yfact=y_scaler,\n",
        "                                                   origin=(0, 0, 0))\n",
        "\n",
        "    train_mask = mask_for_polygons(train_polygons_scaled, image_size)\n",
        "    if image_type =='3':\n",
        "        X = image.reshape(-1, 3).astype(np.float32)\n",
        "    if image_type =='M':\n",
        "        X = image.reshape(-1, 8).astype(np.float32)\n",
        "    if image_type =='A':\n",
        "        X = image.reshape(-1, 8).astype(np.float32)\n",
        "    if image_type =='P':\n",
        "        X = image\n",
        "    y = train_mask.reshape(-1)\n",
        "    return train_mask*class_type\n",
        "\n",
        "def most_common(lst):\n",
        "    data = Counter(lst)\n",
        "    return data.most_common(1)[0][0]\n",
        "\n",
        "def train_image_to_feature(image_id):\n",
        "    # original image tif\n",
        "    image = image_to_array(image_id)\n",
        "    # Segmented training image mask\n",
        "    image_segment_mask = tiff.imread('../input/segment/'+image_id+'_M_SEG.tif')\n",
        "    # Training image mask by classes\n",
        "    image_class_mask = np.max([ image_to_train(image_id,c,'M') for c in CLASSES],axis=0)\n",
        "\n",
        "    start = time.time()\n",
        "    # for each segment, set its class as the one which has most pixels\n",
        "    segment_class_mask=np.zeros(image_segment_mask.shape)\n",
        "    segment_ids = np.unique(image_segment_mask)\n",
        "    labels = []\n",
        "    features = []\n",
        "    for segment_id in segment_ids:\n",
        "    #         segment_class_mask[image_segment_mask==segment_id] = \n",
        "        # Labels\n",
        "        labels.append(most_common(image_class_mask[image_segment_mask==segment_id]))\n",
        "        # Features\n",
        "        segment_pixels = image[image_segment_mask==segment_id]\n",
        "        features.append(get_poly_features(segment_pixels))\n",
        "\n",
        "    return features, labels, segment_ids  "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "6696f85a-5853-498e-e359-abf679a117d1"
      },
      "outputs": [],
      "source": [
        "## Feature extraction for training images\n",
        "train_image_ids = TRAIN_WKT.ImageId.unique()\n",
        "for image_id in train_image_ids:\n",
        "    start = time.time()\n",
        "    features, class_types, segment_ids= train_image_to_feature(image_id)\n",
        "    feature_df = pd.DataFrame(features\n",
        "                ,columns=['b1_min','b1_max','b1_mean','b1_variance','b1_skewness','b1_kurtosis',\n",
        "                          'b2_min','b2_max','b2_mean','b2_variance','b2_skewness','b2_kurtosis',\n",
        "                          'b3_min','b3_max','b3_mean','b3_variance','b3_skewness','b3_kurtosis',\n",
        "                          'b4_min','b4_max','b4_mean','b4_variance','b4_skewness','b4_kurtosis',\n",
        "                          'b5_min','b5_max','b5_mean','b5_variance','b5_skewness','b5_kurtosis',\n",
        "                          'b6_min','b6_max','b6_mean','b6_variance','b6_skewness','b6_kurtosis',\n",
        "                          'b7_min','b7_max','b7_mean','b7_variance','b7_skewness','b7_kurtosis',\n",
        "                          'b8_min','b8_max','b8_mean','b8_variance','b8_skewness','b8_kurtosis'\n",
        "                         ])\n",
        "    feature_df['segment_id'] = segment_ids\n",
        "    feature_df['image_id'] = image_id\n",
        "    feature_df['class_type'] = class_types\n",
        "    feature_df.to_csv('../input/feature_train/'+image_id+'.csv',index = False)\n",
        "    print ('Feature extraction for image %s finished in %d seconds' % (image_id,time.time() - start))\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "b0d22e61-254e-dd06-9edc-aa877e366a1d"
      },
      "source": [
        "### Load training data"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "ca2c9bac-b19e-bd14-6592-71cfa35ef007"
      },
      "outputs": [],
      "source": [
        "train_df = pd.DataFrame()\n",
        "for image_id in train_image_ids:\n",
        "#     print (image_id)\n",
        "    train_df = pd.concat([train_df, pd.read_csv('../input/feature_train/'+image_id+'.csv')],axis=0)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "d8e8897a-d2d1-18ea-819e-b199d0673af5"
      },
      "outputs": [],
      "source": [
        "full_cols = train_df.columns.tolist()\n",
        "full_cols.remove('segment_id')\n",
        "full_cols.remove('image_id')\n",
        "full_cols.remove('class_type')\n",
        "\n",
        "target = 'class_type'"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "10ab7617-7b7a-9c32-fc4b-3d804de4ea0c"
      },
      "source": [
        "## Training"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "f865fb00-3fd2-71d8-a6a3-db20b27b9718"
      },
      "outputs": [],
      "source": [
        "from sklearn.model_selection import train_test_split\n",
        "import xgboost as xgb\n",
        "\n",
        "start = time.time()\n",
        "clf = xgb.XGBClassifier(n_estimators=720, learning_rate = 0.1, max_depth=5)\n",
        "clf.fit(train_df[full_cols].values,train_df[target].values)\n",
        "print (time.time()-start)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "456d8358-c2d2-48e6-65be-5b2a3042ee8a"
      },
      "source": [
        "## Prediction"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "ba647e84-dfbc-81f7-a89e-8d1df7e5f305"
      },
      "outputs": [],
      "source": [
        "from matplotlib import pyplot as plt\n",
        "from matplotlib import colors\n",
        "\n",
        "\n",
        "def pixels_to_poly(image,mask):\n",
        "    poly=[]\n",
        "    for vec in rasterio.features.shapes(image,mask):\n",
        "        poly.append(shapely.geometry.geo.shape(vec[0]))\n",
        "    poly = MultiPolygon(poly)\n",
        "    return poly\n",
        "\n",
        "def predict_image(image_id, clf, full_cols,plot_image = False):\n",
        "    test_df = pd.read_csv('../input/feature/'+image_id+'.csv')\n",
        "\n",
        "\n",
        "    test_x = test_df[full_cols].values\n",
        "    test_segment_id = test_df['segment_id'].values\n",
        "\n",
        "    pred_test_segment_y = clf.predict(test_x)\n",
        "    start = time.time()\n",
        "    test_segment_mask = tiff.imread('../input/segment/{}_M_SEG.tif'.format(image_id))\n",
        "    test_pred_mask = np.zeros(test_segment_mask.shape,dtype=np.int32)\n",
        "\n",
        "    for pid, cls in zip(test_segment_id,pred_test_segment_y):\n",
        "        test_pred_mask[test_segment_mask==pid] = np.int(cls)\n",
        "    \n",
        "    if plot_image:\n",
        "        cmap = colors.ListedColormap([COLORS.get(c,'1') for c in np.unique(test_pred_mask)])\n",
        "        plt.imshow(test_pred_mask, interpolation='none',cmap=cmap)\n",
        "    return test_pred_mask"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "35da67e2-fe67-3099-8933-028f1ee4fa6d"
      },
      "outputs": [],
      "source": [
        "# Make predictions\n",
        "test_df = pd.DataFrame()\n",
        "preds = []\n",
        "for image_id in test_image_ids:\n",
        "    print (\"Predicting image \",image_id)\n",
        "    start = time.time()    \n",
        "    # Make predictions - pixels\n",
        "    pred_image = predict_image(image_id,clf,full_cols, plot_image = False)\n",
        "\n",
        "    image_shape = pred_image.shape[:2]\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",
        "\n",
        "    x_scaler, y_scaler = get_scalers(image_shape, x_max, y_min)   \n",
        "\n",
        "    ##Generate polygons from pixels\n",
        "    for cls in CLASSES:\n",
        "        polygons = pixels_to_poly(pred_image,pred_image==cls)\n",
        "        # Scale polygons\n",
        "        polygons = shapely.affinity.scale(polygons, xfact=1 / x_scaler, yfact=1 / y_scaler, origin=(0, 0, 0))\n",
        "        preds.append([image_id,cls,polygons])\n",
        "    print (\"Predictions for image finishend in %d seconds\" % (time.time()-start))\n",
        "        \n",
        "    "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "5dbf2706-4088-1fb5-8046-ba1d4897ff0f"
      },
      "source": [
        "## Make submission"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "575c4df4-98a6-4f03-98be-12bee1cebbb6"
      },
      "outputs": [],
      "source": [
        "preds_df=pd.DataFrame(preds, columns = ['ImageId','ClassType','MultipolygonWKT'])\n",
        "\n",
        "## Convert polygon precision - this is to reduce the volume of submission data as well as chance of errors for submission \n",
        "\n",
        "preds_df['MultipolygonWKT'] = preds_df['MultipolygonWKT'].\\\n",
        "apply(lambda x:x.simplify(0.0001, preserve_topology=False)).\\\n",
        "apply(lambda x:x if x.is_valid else x.buffer(0)).\\\n",
        "apply(lambda x:shapely.wkt.dumps(x,rounding_precision=5))   \n",
        "\n",
        "\n",
        "## This step is to ensure output will have the same sequence as sample submission\n",
        "output_df = pd.merge(sample_submission[['ImageId','ClassType']],preds_df, how = 'left', on = ['ImageId','ClassType'])\n",
        "\n",
        "## Final output\n",
        "output_df[['ImageId','ClassType','MultipolygonWKT']].to_csv(\"../output/submission.csv\", index=False)"
      ]
    }
  ],
  "metadata": {
    "_change_revision": 0,
    "_is_fork": false,
    "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
}