{"nbformat":4,"nbformat_minor":1,"metadata":{"kernelspec":{"name":"conda-root-py","language":"python","display_name":"Python [conda root]"},"language_info":{"mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"version":"3.5.4"}},"cells":[{"cell_type":"markdown","source":"

[dstl] Satellite Imagery Feature Detection Project

","metadata":{}},{"cell_type":"markdown","source":"Project is done as part of the [dstl] Satellite Imagery Feature Detection Kaggle Competition. This work seeks to better understand how to do semantic segmentation (a category of Computer Vision) in Python. Work is inspired by the great kernels from participants of this competition. See the Github \"ReadMe\" for a list of kernels inspired from.\n

Data

\nAs the image dataset was given in GeoTiff format, we required the use of tifffile and matplotlib to visualize them. In the same way, As the training WKT gave the information as MultipolygonWKT, we needed shapely (or geopandas) to keep track of the geometric form of the polygon classes.","metadata":{}},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"import csv\nimport sys\nimport cv2\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport matplotlib.image as mpimg\nimport matplotlib\nfrom matplotlib.patches import Polygon as polyg"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"import tifffile as tiff\nfrom shapely.wkt import loads\nimport shapely\nfrom shapely.geometry import MultiPolygon, Polygon\nimport shapely.affinity\nimport shapely.wkt\nfrom collections import defaultdict\n#from PIL import Image"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"from sklearn.linear_model import SGDClassifier\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.metrics import average_precision_score\nimport skimage\nfrom skimage.feature import greycomatrix, greycoprops\nfrom skimage.filters import sobel"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"import warnings\nwarnings.filterwarnings('ignore')"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"%matplotlib inline"},{"cell_type":"markdown","source":"

Lists and variables

","metadata":{}},{"cell_type":"markdown","source":"Because the dataset we are working with in this project is fairly large, the following cell is needed to not get a \"field larger than field limit\" Error when using the csv module.","metadata":{}},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"#to load large files\nmaxInt = sys.maxsize\ndecrement = True\n\nwhile decrement:\n # decrease the maxInt value by factor 10 \n # as long as the OverflowError occurs.\n decrement = False\n try:\n csv.field_size_limit(maxInt)\n except OverflowError:\n maxInt = int(maxInt/10)\n decrement = True"},{"cell_type":"markdown","source":"The following dictionaries map the class numerical value to its filename. ClassType can be found in the train_wkt.csv document.","metadata":{}},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"CLASSES = {\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}\nfilename_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\n}"},{"cell_type":"markdown","source":"

Loading the dataset

","metadata":{}},{"cell_type":"markdown","source":"Here, we are loading our datasets (both the training and the grid sizes files as pandas dataframes). The column name is changed such that ImageId becomes the common key between the two tables.","metadata":{}},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{},"source":"grid_sizes = pd.read_csv('./input/grid_sizes.csv')\ngrid_sizes.rename(columns={'Unnamed: 0':'ImageId'}, inplace=True)\ngrid_sizes.head()"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"scrolled":true},"source":"wkt_v4 = pd.read_csv('./input/train_wkt_v4.csv')\nwkt_df = pd.read_csv('./input/train_wkt_v4.csv')\nwkt_df.loc[wkt_df.MultipolygonWKT =='MULTIPOLYGON EMPTY','MultipolygonWKT'] = np.nan\nwkt_df = wkt_df[wkt_df.MultipolygonWKT != 'MULTIPOLYGON EMPTY']\nwkt_df.head()"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"#Define the variables we care about in this study\nallImageIds = grid_sizes.ImageId.unique()\nallClassTypes = wkt_df.ClassType.unique()\nunique_set_ImageIds = grid_sizes.ImageId.astype(str).str[0:4].unique()\nImageId = '6120_2_2'\nPoly_class = 1 #np.asscalar(allClassTypes[0])"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"#this is a common function that will be called often later to redeclare the variables when needed\ndef variables(ImageId, Poly_class):\n global rgb, im_rgb, x_max, y_min, train_polygons\n rgb = tiff.imread('./data/three_band/{}.tif'.format(ImageId))\n im_rgb = tiff.imread('./data/three_band/{}.tif'.format(ImageId)).transpose([1, 2, 0])\n x_max = y_min = None\n for i, row in grid_sizes.iterrows():\n if row['ImageId'] == '6120_2_2':\n x_max, y_min = float(row['Xmax']), float(row['Ymin'])\n break\n train_polygons = None\n for i, row in wkt_v4.iterrows():\n if row['ImageId'] == '6120_2_2' and row['ClassType'] == Poly_class:\n train_polygons = loads(row['MultipolygonWKT'])\n break"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"variables(ImageId, Poly_class)"},{"cell_type":"markdown","source":"

Visualizing

","metadata":{}},{"cell_type":"markdown","source":"The most classic way to visualize a satellite image is to simply plot it using the tifffine module.","metadata":{}},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{},"source":"tiff.imshow(im_rgb)"},{"cell_type":"markdown","source":"However a way better way to visualize the satellite image is by changing the pixel value of images into 5 to 95%, which is a robust way to eliminate the outlier pixel coloring.","metadata":{}},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{},"source":"def scale_percentile(matrix):\n w, h, d = matrix.shape\n matrix = np.reshape(matrix, [w * h, d]).astype(np.float64)\n # Get 5nd and 95th percentile\n mins = np.percentile(matrix, 6, axis=0)\n maxs = np.percentile(matrix, 94, axis=0) - mins\n matrix = (matrix - mins[None, :]) / maxs[None, :]\n matrix = np.reshape(matrix, [w, h, d])\n matrix = matrix.clip(0, 1)\n return matrix\ntiff.imshow(255*scale_percentile(im_rgb)); #choosing size"},{"cell_type":"markdown","source":"

Some Data Analytics on the Dataset\"

","metadata":{}},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"# PolygonList is the list of all classes inside the image with its MultipolygonWKT values\npolygonsList = {}\nimage = wkt_v4[wkt_df.ImageId == ImageId]\nfor cType in image.ClassType.unique():\n polygonsList[cType] = loads(image[image.ClassType == cType].MultipolygonWKT.values[0])"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"scrolled":true},"source":"# number of objects on the image by type\nfor p in polygonsList:\n print(\"Class type: [{}] {} -> Number of objects: {}\".format(p, CLASSES[p],len(polygonsList[p].geoms)))"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"scrolled":false},"source":"# Do the same as previously but for all existing ImageId identified by the training WKT\n# Below is a pivot table where the columns are class types and the rows are for each ImageId\nwkt_v4['polygons'] = wkt_v4.apply(lambda row: loads(row.MultipolygonWKT),axis=1)\nwkt_v4['nPolygons'] = wkt_v4.apply(lambda row: len(row['polygons'].geoms),axis=1)\npvt = wkt_v4.pivot(index='ImageId', columns='ClassType', values='nPolygons')\npvt"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"scrolled":true},"source":"#Heatmap for the frequency of class types\nfig, ax = plt.subplots(figsize=(10, 4))\nax.set_aspect('equal')\nplt.imshow(pvt.T, interpolation='nearest', cmap=plt.cm.Blues, extent=[0,22,10,1])\nplt.yticks(np.arange(1, 11, 1.0))\nplt.title('Number of objects by type')\nplt.ylabel('Class Type')\nplt.xlabel('Image')\nplt.colorbar()\nplt.show()"},{"cell_type":"markdown","source":"For this dataset, we see that the trees (having classtype 5) is the most identified object type for the labeled images. The dataset we got represents satellite images of 1km X 1km rural and urban areas. If you go check our second kernel where we visualized 7 different channels, we can see that trees and natural elements are dominating in most satellite images. The other major class types represented in the satellite images are buildings and vehicles. Waterways and routes are present at much smaller scale.","metadata":{}},{"cell_type":"markdown","source":"

Show Polygon Masks

\nTo help with the training of the dataset, some of the satellite has their polygon classes labelled in shapely multipolygon format. This implies that for specific satellite images, we can visualize it using Polygon masks. Three different methods are presented below.","metadata":{}},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"def show_mask(m):\n # hack for nice display\n tiff.imshow(255 * np.stack([m, m, m]));"},{"cell_type":"markdown","source":"

Visualizing the Polygons Method 1: Using matplotlib.patches.Polygon

\nThe matplotlib.patches.Polygon returns a general polygon patch which can be plotted.","metadata":{}},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{},"source":"# plotting by class type (p variable)\nfig, ax = plt.subplots(2, 5, figsize=(10, 5))\nfor i, p in enumerate(polygonsList):\n if i<5:\n for polygon in polygonsList[p]:\n #color=plt.cm.Set1(p*10)\n mpl_poly = polyg(np.array(polygon.exterior), lw=5, alpha=0.5)\n ax[0, i].add_patch(mpl_poly)\n ax[0, i].relim()\n ax[0, i].autoscale_view()\n ax[0, i].axis('off')\n else:\n for polygon in polygonsList[p]:\n #color=plt.cm.Set1(p*10)\n mpl_poly = polyg(np.array(polygon.exterior), lw=5, alpha=0.5)\n ax[1, i-5].add_patch(mpl_poly)\n ax[1, i-5].relim()\n ax[1, i-5].autoscale_view()\n ax[1, i-5].axis('off')\nplt.tight_layout()\nfig.subplots_adjust(wspace=0, hspace=0)"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{},"source":"# plot all the polygons on top of each other using matplotlib\nfig, ax = plt.subplots(figsize=(10, 10))\n# plotting, color by class type (p variable)\nfor i, p in enumerate(polygonsList):\n for polygon in polygonsList[p]:\n #color=plt.cm.Set1(p*10)\n mpl_poly = polyg(np.array(polygon.exterior), lw=0, alpha=0.3)\n ax.add_patch(mpl_poly)\n ax.relim()\n ax.autoscale_view()\nfig.subplots_adjust(wspace=0, hspace=0)"},{"cell_type":"markdown","source":"

Visualizing the Polygons Method 2: Using OpenCV.fillPoly to Fill

","metadata":{}},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"def _convert_coordinates_to_raster(coords, img_size, xymax):\n Xmax,Ymax = xymax\n H,W = img_size\n W1 = 1.0*W*W/(W+1)\n H1 = 1.0*H*H/(H+1)\n xf = W1/Xmax\n yf = H1/Ymax\n coords[:,1] *= yf\n coords[:,0] *= xf\n coords_int = np.round(coords).astype(np.int32)\n return coords_int\n\ndef _get_xmax_ymin(grid_sizes_panda, imageId):\n xmax, ymin = grid_sizes_panda[grid_sizes_panda.ImageId == imageId].iloc[0,1:].astype(float)\n return (xmax,ymin)\n\ndef _get_polygon_list(wkt_list_pandas, imageId, cType):\n df_image = wkt_list_pandas[wkt_list_pandas.ImageId == imageId]\n multipoly_def = df_image[df_image.ClassType == cType].MultipolygonWKT\n polygonList = None\n if len(multipoly_def) > 0:\n assert len(multipoly_def) == 1\n polygonList = loads(multipoly_def.values[0])\n return polygonList\n\ndef _get_and_convert_contours(polygonList, raster_img_size, xymax):\n perim_list = []\n interior_list = []\n if polygonList is None:\n return None\n for k in range(len(polygonList)):\n poly = polygonList[k]\n perim = np.array(list(poly.exterior.coords))\n perim_c = _convert_coordinates_to_raster(perim, raster_img_size, xymax)\n perim_list.append(perim_c)\n for pi in poly.interiors:\n interior = np.array(list(pi.coords))\n interior_c = _convert_coordinates_to_raster(interior, raster_img_size, xymax)\n interior_list.append(interior_c)\n return perim_list,interior_list\n\ndef _plot_mask_from_contours(raster_img_size, contours, class_value = 1):\n img_mask = np.zeros(raster_img_size, np.uint8)\n if contours is None:\n return img_mask\n perim_list,interior_list = contours\n cv2.fillPoly(img_mask, perim_list, class_value)\n cv2.fillPoly(img_mask, interior_list, 0)\n return img_mask\n\ndef generate_mask_for_image_and_class(raster_size, imageId, class_type, grid_sizes_panda, wkt_list_pandas):\n xymax = _get_xmax_ymin(grid_sizes_panda,imageId)\n polygon_list = _get_polygon_list(wkt_list_pandas,imageId,class_type)\n contours = _get_and_convert_contours(polygon_list,raster_size,xymax)\n mask = _plot_mask_from_contours(raster_size,contours,1)\n return mask"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"scrolled":false},"source":"# Visualizing for one polygon class, notice that 1 is for buildings\nmask = generate_mask_for_image_and_class((500,500), ImageId, 1, grid_sizes, wkt_v4)\nshow_mask(mask)"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"scrolled":false},"source":"#combine all the polygon masks together\nset_of_mask = dict()\nmask_test = np.zeros((500, 500))\nfor i in range(0,9):\n mask = generate_mask_for_image_and_class((500,500), ImageId, i , grid_sizes, wkt_v4)\n set_of_mask[i] = mask*255/9*i\n mask_test = mask_test + mask*255/9*i\nx = tiff.imshow(255 * mask_test)"},{"cell_type":"markdown","source":"

Visualizing the Polygons Method 3: Using OpenCV.polylines to Render

","metadata":{}},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"def adjust_contrast(x): \n for i in range(3):\n x[:,:,i] = stretch2(x[:,:,i])\n return x.astype(np.uint8)\n\ndef stretch2(band, lower_percent=2, higher_percent=98):\n a = 0 #np.min(band)\n b = 255 #np.max(band)\n c = np.percentile(band, lower_percent)\n d = np.percentile(band, higher_percent) \n out = a + (band - c) * (b - a) / (d - c) \n out[outb] = b\n return out"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"def truth_polys(image_id, class_id, W, H):\n rows = wkt_df.loc[(wkt_df.ImageId==image_id) & (wkt_df.ClassType==class_id), 'MultipolygonWKT']\n mp = loads(rows.values[0])\n xmax, ymin = grid_sizes[grid_sizes.ImageId == ImageId].iloc[0,1:].astype(float) \n W_ = W * (W/(W+1.))\n H_ = H * (H/(H+1.))\n x_scaler = W_ / xmax\n y_scaler = H_ / ymin\n return shapely.affinity.scale(mp, xfact = x_scaler, yfact= y_scaler, origin=(0,0,0)) "},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"scrolled":false},"source":"rgb = tiff.imread('./data/three_band/{}.tif'.format(ImageId))\nrgb = np.rollaxis(rgb, 0, 3)\nx = adjust_contrast(rgb).copy()\nH=len(x); W=len(x[0])\n#Read Polygons\npolys = truth_polys(ImageId, Poly_class, W, H)\nint_vertices=lambda x: np.array(x).round().astype(np.int32)\nfor poly_id, poly in enumerate(polys):\n xys=int_vertices(poly.exterior.coords)\n cv2.polylines(x,[xys],True,(255,0,0),3)\n for pi in poly.interiors:\n ixys=int_vertices(pi.coords)\n cv2.polylines(x,[ixys],True,(255,0,0),3) \n#Plot\nfig, ax = plt.subplots(figsize=(9,9))\nax.imshow(x)"},{"cell_type":"markdown","source":"

Bonus: Visualize the whole satellite channel together!

\nAs a bonus, we defined functions to visualize the whole channel's 25 satellite images together. We also wrote a different Jupyter notebook doing exactly that for many different channels.","metadata":{}},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"def plot_5_by_5(list_of_id):\n fig, ax = plt.subplots(5, 5, figsize=(20, 20))\n for k in list_of_id:\n i = int(k[5])\n j = int(k[7])\n temp = tiff.imread('./data/three_band/{}.tif'.format(k))\n temp = np.rollaxis(temp, 0, 3)\n y = adjust_contrast(temp).copy()\n ax[i, j].imshow(y)\n ax[i, j].axis('off')\n fig.subplots_adjust(wspace=0, hspace=0)"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{},"source":"%%time\n#plot all images: plot 5 images (top and bottom of the image we are studying)\ntemp_set_id = ImageId[0:4]\nlist_of_id = []\nfor i in allImageIds:\n if i[0:4] == temp_set_id:\n list_of_id.append(i)\nplot_5_by_5(list_of_id)"},{"cell_type":"markdown","source":"

Prediction: Logistic regression classifier

\nIn the context of this project, we have labeled training set data. It hence makes sense to use a supervised machine learning algorithm to do the prediction. Although Computer Vision has benefitted extensively from the advancement of Deep Learning and its libraries (TensorFlow, Keras, CNN, etc), we decided to follow the simpler method introduced on the Kaggle competition kernel. Hence a logistic regression classifier that takes as input the values of the multipolygon from the training WKT to was used to do the prediction.\n

Train dataset

","metadata":{}},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"#scale the image\ndef get_scalers():\n h, w = im_rgb.shape[:2] # 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"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"#create a mask from polygons\ndef mask_for_polygons(polygons):\n img_mask = np.zeros(im_rgb.shape[:2], np.uint8)\n if not polygons:\n return img_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(img_mask, exteriors, 1)\n cv2.fillPoly(img_mask, interiors, 0)\n return img_mask"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"x_scaler, y_scaler = get_scalers()\ntrain_polygons_scaled = shapely.affinity.scale(train_polygons, xfact=x_scaler, yfact=y_scaler, origin=(0, 0, 0))\ntrain_mask = mask_for_polygons(train_polygons_scaled)"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{},"source":"#Get the binary image\nshow_mask(train_mask)"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"scrolled":true},"source":"xs = im_rgb.reshape(-1, 3).astype(np.float32)\nys = train_mask.reshape(-1)\npipeline = make_pipeline(StandardScaler(), SGDClassifier(loss='log'))\n\nprint('training...')\n# do not care about overfitting here\npipeline.fit(xs, ys)\npred_ys = pipeline.predict_proba(xs)[:, 1]\nprint('average precision', average_precision_score(ys, pred_ys))\npred_mask = pred_ys.reshape(train_mask.shape)"},{"cell_type":"markdown","source":"

Test Images

","metadata":{}},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{},"source":"# Predicted mask\nshow_mask(pred_mask)\n# Predicted mask with threshold\nthreshold = 0.3\npred_binary_mask = pred_mask >= threshold\nshow_mask(pred_binary_mask)"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{},"source":"#create polygons from bit masks\ndef mask_to_polygons(mask, epsilon=10., min_area=10.):\n image, contours, hierarchy = cv2.findContours(\n ((mask == 1) * 255).astype(np.uint8),\n cv2.RETR_CCOMP, cv2.CHAIN_APPROX_TC89_KCOS)\n # create approximate contours to have reasonable submission size\n approx_contours = [cv2.approxPolyDP(cnt, epsilon, True)\n for cnt in contours]\n if not contours:\n return MultiPolygon()\n # associate parent and child contours\n cnt_children = defaultdict(list)\n child_contours = set()\n assert hierarchy.shape[0] == 1\n for idx, (_, _, _, parent_idx) in enumerate(hierarchy[0]):\n if parent_idx != -1:\n child_contours.add(idx)\n cnt_children[parent_idx].append(approx_contours[idx])\n # create actual polygons filtering by area\n all_polygons = []\n for idx, cnt in enumerate(approx_contours):\n if idx not in child_contours and cv2.contourArea(cnt) >= min_area:\n assert cnt.shape[1] == 1\n poly = Polygon(\n shell=cnt[:, 0, :],\n holes=[c[:, 0, :] for c in cnt_children.get(idx, [])\n if cv2.contourArea(c) >= min_area])\n all_polygons.append(poly)\n # approximating polygons might have created invalid ones, fix them\n all_polygons = MultiPolygon(all_polygons)\n if not all_polygons.is_valid:\n all_polygons = all_polygons.buffer(-1)\n if all_polygons.type == 'Polygon':\n all_polygons = MultiPolygon([all_polygons])\n return all_polygons\n\n#Rendering Shapes: Jaccard on the pixel level\ntp, fp, fn = (( pred_binary_mask & train_mask).sum(),\n ( pred_binary_mask & ~train_mask).sum(),\n (~pred_binary_mask & train_mask).sum())\nprint('Pixel jaccard', tp / (tp + fp + fn))\n#prediction to polygon and back to mask! This shows that the shapes has been rendered!\npred_polygons = mask_to_polygons(pred_binary_mask)\npred_poly_mask = mask_for_polygons(pred_polygons)\nshow_mask(pred_poly_mask)"},{"cell_type":"markdown","source":"

Visualize the whole channel's prediction

\nThe excellent result can be explained by the allowance of overfitting. Indeed, in this case, we were looking at the file \"6120_2_2\" which is a labeled geotiff file. It hence becomes important to retest the algorithm on cases where the multipolygon values were not defined. This is done for all files in the channel \"6120\".

\nAs we can see below, the prediction values become a lot lower, averaging 25%. Although disappoint, it makes sense because the dataset given to do the training had most of the images unlabeled and each satellite image depends on multiple geographical factors and the 2D rendered masks changes accordingly. In fact, out of 450 satellite images in our dataset, only 25 of them were properly labeled for training. Maybe using CNN could result in better result.","metadata":{}},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"def plot_5_by_5_pred(list_of_id, list_of_img):\n fig, ax = plt.subplots(5, 5, figsize=(20, 20))\n n = len(list_of_id)\n for i in range(0, n):\n ax[int(list_of_id[i][5]), int(list_of_id[i][7])].imshow(list_of_img[i])\n ax[int(list_of_id[i][5]), int(list_of_id[i][7])].axis('off')\n fig.subplots_adjust(wspace=0, hspace=0)"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":"def predicting_values(im_ID):\n variables(im_ID, Poly_class)\n x_scaler, y_scaler = get_scalers()\n train_polygons_scaled = shapely.affinity.scale(train_polygons, xfact=x_scaler, yfact=y_scaler, origin=(0, 0, 0))\n train_mask = mask_for_polygons(train_polygons_scaled)\n xs = im_rgb.reshape(-1, 3).astype(np.float32)\n ys = train_mask.reshape(-1)\n pipeline = make_pipeline(StandardScaler(), SGDClassifier(loss='log'))\n # overfitting is not considered\n pipeline.fit(xs, ys)\n pred_ys = pipeline.predict_proba(xs)[:, 1]\n print('training average precision for {} is {}'.format(im_ID, average_precision_score(ys, pred_ys)))\n pred_mask = pred_ys.reshape(im_rgb.shape[:2])\n # Predicted mask with threshold into binary image\n threshold = 0.3\n pred_binary_mask = pred_mask >= threshold\n return pred_binary_mask"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{},"source":"%%time\n#plot all images: plot 5 images (top and bottom of the image we are studying)\ntemp_set_id = '6120'\nlist_of_id = []\nlist_of_img = []\nfor i in allImageIds:\n if i[0:4] == temp_set_id:\n x = predicting_values(i)\n list_of_id.append(i)\n list_of_img.append(x)"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{},"source":"%%time\nplot_5_by_5_pred(list_of_id, list_of_img)"},{"cell_type":"code","outputs":[],"execution_count":null,"metadata":{"collapsed":true},"source":""}]}