{"cells":[{"cell_type":"markdown","metadata":{"_cell_guid":"9940a12e-502e-630c-4f89-c3379a5de225"},"source":"**Polygon count can be misleading! The polygon count of trees is more but their area is not.**\n\nUpdate: Kernel now uses geojson data instead of the WKT data."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"1c2fa704-6eca-0309-7779-46a8ff2fb090"},"outputs":[],"source":"# import some useful packages\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom shapely.wkt import loads as wkt_loads\nimport matplotlib.pyplot as plt\nimport pylab\nimport json\nimport shapely.geometry\nimport os\nimport tensorflow as tf"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"486d2195-754c-3a6b-a54c-2cf949892148"},"outputs":[],"source":"# Input data files are available in the \"../input/\" directory.\ninDir = '../input'\ngeojsonDir = 'train_geojson_v3'\n\n# Give short names to object types\nCLASSES = {\n        1 : 'Bldg',\n        2 : 'Struct',\n        3 : 'Road',\n        4 : 'Track',\n        5 : 'Trees',\n        6 : 'Crops',\n        7 : 'Fast H20',\n        8 : 'Slow H20',\n        9 : 'Truck',\n        10 : 'Car',\n        }\nCLASSES_R = {}\nfor k in CLASSES:\n    CLASSES_R[CLASSES[k]] = k\nCOLORS = {\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\ndef load_tables():\n    # read the training data from train_wkt_v3.csv\n    df = pd.read_csv(inDir + '/train_wkt_v3.csv',\n            names=['SceneId', 'ClassType', 'MultipolygonWKT'], skiprows=1)\n\n    # grid size will also be needed later..\n    gs = pd.read_csv(inDir + '/grid_sizes.csv',\n            names=['SceneId', 'Xmax', 'Ymin'], skiprows=1)\n    return df, gs"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"f3d7c472-a0bc-5a72-5bb9-3f071bbec502"},"outputs":[],"source":"df, gs = load_tables()\ndf.head()\ngs.head()\nallSceneIds = gs.SceneId.unique()\ntrainSceneIds = np.sort(df.SceneId.unique())"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"1e66731f-9e08-a9db-5e1a-51b7eae9b503"},"outputs":[],"source":"filename_to_classType = {\n        '001_MM_L2_LARGE_BUILDING':1,\n        '001_MM_L3_RESIDENTIAL_BUILDING':1,\n        '001_MM_L3_NON_RESIDENTIAL_BUILDING':1,\n        '001_MM_L5_MISC_SMALL_STRUCTURE':2,\n        '002_TR_L3_GOOD_ROADS':3,\n        '002_TR_L4_POOR_DIRT_CART_TRACK':4,\n        '002_TR_L6_FOOTPATH_TRAIL':4,\n        '006_VEG_L2_WOODLAND':5,\n        '006_VEG_L3_HEDGEROWS':5,\n        '006_VEG_L5_GROUP_TREES':5,\n        '006_VEG_L5_STANDALONE_TREES':5,\n        '007_AGR_L2_CONTOUR_PLOUGHING_CROPLAND':6,\n        '007_AGR_L6_ROW_CROP':6, \n        '008_WTR_L3_WATERWAY':7,\n        '008_WTR_L2_STANDING_WATER':8,\n        '003_VH_L4_LARGE_VEHICLE':9,\n        '003_VH_L5_SMALL_VEHICLE':10,\n        '003_VH_L6_MOTORBIKE':10}\ncType_to_Filename = {}\nfor cType in range(1,11):\n    cType_to_Filename[cType] = [k for k in filename_to_classType if filename_to_classType[k] == cType]"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"5ba44e05-62e5-757f-dea3-00e1d91d80a5"},"outputs":[],"source":"# Basic functions to work with polygons\ndef get_grid_size(sceneId):\n    '''\n    Get the grid size of the scene\n\n    Parameters\n    ----------\n    sceneId : str\n        sceneId as used in grid_size.csv\n\n    Returns\n    -------\n    (xmax, ymin) : float\n    '''\n    xmax, ymin = gs[gs.SceneId == sceneId].iloc[0,1:].astype(float)\n    return xmax, ymin\n\n\ndef get_grid_area(sceneId):\n    '''\n    Get the area of the grid of the scene\n    '''\n    xmax, ymin = get_grid_size(sceneId)\n    return np.abs(xmax*ymin)\n\n\ndef get_polygons(sceneId):\n    '''\n    Parameters\n    ----------\n    sceneId : str\n        sceneId like \"6010_0_4\"\n\n    Returns\n    -------\n    polygonsList : dict\n        Keys are CLASSES\n        Values are shapely polygons\n    '''\n    # df_scene = df[df.SceneId == sceneId]\n    \n    polygonsList = {}\n    for cType in CLASSES.keys():\n        # WKT version:\n        # polygonsList[cType] = wkt_loads(df_scene[df_scene.ClassType == cType].MultipolygonWKT.values[0])\n        \n        # geojson version:\n        polygonsList[cType] = []\n        for filename in cType_to_Filename[cType]:\n            fullpath = '/'.join([inDir, geojsonDir, sceneId, filename + '.geojson'])\n            if not os.path.isfile(fullpath):\n                continue\n            with open(fullpath, 'r') as fp:\n                # print('DEBUG Opening file {}'.format(fullpath))\n                fcObj = json.load(fp)\n                fcList = fcObj['features']\n                for fc in fcList:\n                    geom = fc['geometry']\n                    polygonsList[cType].append(shapely.geometry.shape(geom))\n    return polygonsList"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"4fe8afef-12be-5222-0033-110f0661493b"},"outputs":[],"source":"# The main function for this kernel\ndef get_stats_polygons(polygonsList, image_area=1):\n    '''\n    Get stats from polygonsList\n    '''\n    count = {}\n    totalArea = {}\n    meanArea = {}\n    stdArea = {}\n    for cType in polygonsList:\n        count[cType] = len(polygonsList[cType])\n        if count[cType] > 0:\n            totalArea[cType] = np.sum ([polygon.area for polygon in polygonsList[cType]]) / image_area*100\n            meanArea [cType] = np.mean([polygon.area for polygon in polygonsList[cType]]) / image_area*100\n            stdArea  [cType] = np.std ([polygon.area for polygon in polygonsList[cType]]) / image_area*100\n        else:\n            totalArea[cType] = 0\n            meanArea [cType] = np.nan\n            stdArea  [cType] = np.nan\n\n    return pd.DataFrame({\n            'CLASS' : CLASSES,\n            'counts' : count,\n            'totalAreas' : totalArea,\n            'meanAreas' : meanArea,\n            'stdAreas' : stdArea,\n            })"},{"cell_type":"markdown","metadata":{"_cell_guid":"a5a2026c-b5e4-546d-2114-75455cd9a244"},"source":"Now, we are ready to collect some stats from the polygons of the training data."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"16d20f20-a52f-9fc4-a8ae-57176beaec35"},"outputs":[],"source":"def collect_stats():\n    pStatsList = []\n    for sceneId in trainSceneIds:\n        polyList = get_polygons(sceneId)\n        pStats = get_stats_polygons(polyList, image_area=get_grid_area(sceneId))\n        pStats['SceneId'] = sceneId\n        pStatsList.append(pStats)\n        print('Stats loaded for ', sceneId)\n        # print(pStats)\n    return pd.concat(pStatsList)\npStats = collect_stats()\npStats"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"c30647cc-b084-7eba-54be-dd6a0b173b3e"},"outputs":[],"source":"# Function to pivot and plot\ndef plot_stats(pStats, values, title):\n    pvt = pStats.pivot(index='CLASS', columns='SceneId', values=values)\n    fig, ax = plt.subplots(figsize=(10,4))\n    ax.set_aspect('equal')\n    plt.imshow(pvt, interpolation='nearest', cmap=plt.cm.plasma, extent=[0, 25, 10, 0])\n    plt.xticks(np.arange(0.5, 25.4, 1))\n    plt.yticks(np.arange(0.5, 10.4, 1))\n    ax.set_xticklabels(np.arange(1,26,1))\n    ax.set_yticklabels(pvt.index)\n    plt.xlabel('Image')\n    plt.ylabel('Class Type')\n    plt.title(title)\n    plt.colorbar()"},{"cell_type":"markdown","metadata":{"_cell_guid":"61f1e391-e41f-17c5-1afd-079a43e81b6e"},"source":"First, let's generate a plot similar to another [existing kernel](https://www.kaggle.com/torrinos/dstl-satellite-imagery-feature-detection/exploration-and-plotting/run/553107)."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"c5bdf517-4c0d-82d9-3990-8fd72f3129c3"},"outputs":[],"source":"plot_stats(pStats, 'counts', 'Number of polygons in image by type')"},{"cell_type":"markdown","metadata":{"_cell_guid":"a70352d7-36e9-2c94-b448-ee683233fdea"},"source":"As has been reported earlier, there are many trees in the training images."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"c0f324d1-7791-3537-8189-86e185f77fcd"},"outputs":[],"source":"# Now, let's plot percentage of area covered by each type\nplot_stats(pStats, 'totalAreas', 'Percentage of area in image by type')\n"},{"cell_type":"markdown","metadata":{"_cell_guid":"15e524a0-4212-59fd-45cb-48ebdb1ea43c"},"source":"While trees still cover a lot images,  some images have high coverage by crops."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"89bd84b6-4426-30db-1b31-1196bf837d20"},"outputs":[],"source":"# Since count is high but percentage is low for tree-type,\n# we expect its' meanArea per polygon is low.\n# Let's verify this.\nplot_stats(pStats, 'meanAreas', 'Average area covered by single polygon')"},{"cell_type":"markdown","metadata":{"_cell_guid":"e3a0bbc0-413c-3bb0-85d3-2f65d8e7be7d"},"source":"As claimed in other kernels, some images are of the same location (possibly different times).  Here's one way to find such pairs by looking at the polygon stats of the training set:"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"e4ea5461-8053-1e75-8a1f-e36c5357d9e1"},"outputs":[],"source":"pvt = pStats.pivot(index='CLASS', columns='SceneId', values='totalAreas')\npairImages = ['6110_1_2', '6140_1_2']\nprint(pvt[pairImages])\nfrom scipy.stats import pearsonr\nprint('Coverage for different classses in {} and {}: {:5.4f}'.format(\n    pairImages[0], pairImages[1], pearsonr(pvt[pairImages[0]],pvt[pairImages[1]])[0]))\n# Todo find more pairs"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"16e0ab03-9f22-09f6-300b-ffb3cdfcb7b0"},"outputs":[],"source":"percAreaCS = np.cumsum(pvt, axis=0)\nimport seaborn as sns\n\n#Set general plot properties\nsns.set_style(\"white\")\nsns.set_context({\"figure.figsize\": (12, 8)})\n\nfor i in range(1,11):\n    cTypeName = percAreaCS.index[-i]\n    cTypeId = CLASSES_R[cTypeName]\n    ax = sns.barplot(x = percAreaCS.columns, y = percAreaCS.iloc[-i],\n                              color = COLORS[cTypeId], label = cTypeName)\nl = plt.legend(loc=2)\n\nsns.despine(left=True)\nax.set_xlabel(\"sceneId\")\nax.set_ylabel(\"%age Area Covered\")    \n_ = ax.set_xticklabels(ax.xaxis.get_majorticklabels(), rotation=-60)\n"},{"cell_type":"markdown","metadata":{"_cell_guid":"69ec4118-b416-a1d7-966f-cfdb6f4d966f"},"source":"* Some bars go beyond 100% because some polygons may overlap."}],"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.5.2"}},"nbformat":4,"nbformat_minor":0}