{"cells":[{"cell_type":"markdown","metadata":{"_cell_guid":"a2c5f3be-61e9-1471-961c-6578f05b9fdc"},"source":"This is a quick demonstration of how to load a TIF, display it, and then overlay polygons of one of the feature types, `POOR_DIRT_CART_TRACK`.  \n\nI used tifffile to load the image because it's easier to install than GDAL, and the files don't have headers anyway so there's no advantage to GDAL.  I used matplotlib to handle polygons, although I should mention that some of my other code uses openCV for some polygon operations."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"45757e81-8479-7383-4ae2-9b50d5234b96"},"outputs":[],"source":"import os\nimport numpy as np\nimport tifffile\nimport json\nfrom PIL import Image\nimport matplotlib\nimport matplotlib.pyplot as plt\nimport pandas as pd\nfrom matplotlib.collections import PatchCollection\n\n%matplotlib inline\n\ndata_dir = '../input'\ngrid_name = '6010_4_2'\n\n# Load grid CSV\ngrid_sizes = pd.read_csv(os.path.join(data_dir, 'grid_sizes.csv'), index_col=0)\ngrid_sizes.ix[grid_name]\n    \n# Load JSON of image overlays\nsh_fname = os.path.join(data_dir, 'train_geojson_v3/%s/002_TR_L4_POOR_DIRT_CART_TRACK.geojson'%grid_name)\nwith open(sh_fname, 'r') as f:\n    sh_json = json.load(f)\n    \n# Load the tif file\nim_fname = os.path.join(data_dir, 'three_band','%s.tif'%grid_name)\ntif_data = tifffile.imread(im_fname).transpose([1,2,0])"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"b593d958-47d4-396d-0fb4-10d03284e54b"},"outputs":[],"source":"def scale_coords(tif_data, grid_name, point):\n    \"\"\"Scale the coordinates of a polygon into the image coordinates for a grid cell\"\"\"\n    w,h,_ = tif_data.shape\n    Xmax, Ymin = grid_sizes.ix[grid_name][['Xmax', 'Ymin']]\n    x,y = point[:,0], point[:,1]\n\n    wp = float(w**2)/(w+1)\n    xp = x/Xmax*wp\n\n    hp = float(h**2)/(h+1)\n    yp = y/Ymin*hp\n\n    return np.concatenate([xp[:,None],yp[:,None]], axis=1)\n\ndef scale_percentile(matrix):\n    \"\"\"Fixes the pixel value range to 2%-98% original distribution of values\"\"\"\n    orig_shape = matrix.shape\n    matrix = np.reshape(matrix, [matrix.shape[0]*matrix.shape[1], 3]).astype(float)\n    \n    # Get 2nd and 98th percentile\n    mins = np.percentile(matrix, 1, axis=0)\n    maxs = np.percentile(matrix, 99, axis=0) - mins\n    \n    matrix = (matrix - mins[None,:])/maxs[None,:]\n    matrix = np.reshape(matrix, orig_shape)\n    matrix = matrix.clip(0,1)\n    return matrix"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"df46e4ec-830d-4056-b2a3-d7d2cdf236d7"},"outputs":[],"source":"# Show the image with the values scaled from 2-98 percentile to make them visible\nfixed_im = scale_percentile(tif_data)\nplt.imshow(fixed_im)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"78115b1c-4629-cbe4-eb87-d80cb8b077af"},"outputs":[],"source":"# Load JSON of image overlays, and convert it into image coordinates\ndef load_overlays(tile_name):\n    \"\"\"Get all of the polygon overlays for a tile.\n    Returns a dict: {LABEL: POLYGON}\"\"\"\n    dirname = os.path.join(data_dir, 'train_geojson_v3/%s/'%tile_name)\n    fnames = [os.path.join(dirname, fname) for fname in \n              os.listdir(dirname) \n              if fname.endswith('.geojson') and not fname.startswith('Grid')]\n    \n    overlays = dict()\n    for fname in fnames:\n        with open(fname, 'r') as f:\n            sh_json = json.load(f)\n        label = sh_json['features'][0]['properties']['LABEL']\n        print(label)\n        \n        polygons = []\n        for sh in sh_json['features']:\n            pts = scale_coords(tif_data, grid_name, np.array(sh['geometry']['coordinates'][0])).squeeze()\n            \n            # Remove badly formatted polygons\n            if not ((len(pts.shape)==2) and (pts.shape[1]==2) and (pts.shape[0] > 2)):\n                continue\n            polygons.append(pts)\n            \n        overlays[label] = polygons\n\n    return overlays\n    \noverlays_poly = load_overlays(grid_name)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"ea638582-fd7b-eca7-401c-ccb2095635a1"},"outputs":[],"source":"# Put polygons on top of the image\n\nfig, ax = plt.subplots(figsize=(20,20))\npatches = []\nfor pts in overlays_poly['GROUP_TREES']:\n    \n    poly = matplotlib.patches.Polygon(pts)\n    patches.append(poly)\n\np = PatchCollection(patches, cmap=matplotlib.cm.jet, alpha=0.4)\n\ncolors = 100*np.random.rand(len(patches))\np.set_array(np.array(colors))\n\nax.imshow(fixed_im)\nax.add_collection(p)\n\nplt.show()"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"4ceddf90-25c1-3238-cbd1-e9489c72814d"},"outputs":[],"source":""}],"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}