<jupyter_start><jupyter_text>Open Cities AI Challenge: Segmenting Buildings for Disaster Resilience This is my solution to [Open Cities AI Challenge: Segmenting Buildings for Disaster Resilience](https://www.drivendata.org/competitions/60/building-segmentation-disaster-resilience/page/150/). Sources:1) How to prepare [patches of data for training](https://medium.com/@anthropoco/how-to-segment-buildings-on-drone-imagery-with-fast-ai-cloud-native-geodata-tools-ae249612c321).2) [Getting started guide](https://colab.research.google.com/drive/1Fv-80b1m-O-0p1g59NDzD82XdgurWlwa) from [johnowhitaker](https://community.drivendata.org/u/johnowhitaker/summary).3) [Fast.AI](https://docs.fast.ai/)Google colab notebook was bootstraped using [colab-bootstrap](https://github.com/zaitra/colab-bootstrap) by [Zaitra](https://zaitra.io).  Install python packages<jupyter_code>!add-apt-repository ppa:ubuntugis/ubuntugis-unstable -y
!apt-get update
!apt-get install python-numpy gdal-bin libgdal-dev python3-rtree
!pip3 install geopandas pillow rasterio descartes solaris rio-tiler pystac<jupyter_output><empty_output><jupyter_text>Imports<jupyter_code># Jupyter notebook related
%reload_ext autoreload
%autoreload 2
%matplotlib inline

# Built-in modules
import sys
import os
import datetime
from enum import Enum
from pathlib import Path


#import solaris as sol
# Basics of Python data handling and visualization
import skimage
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd

import rasterio
from rasterio.transform import from_bounds
from rasterio.windows import Window
from rio_tiler import main as rt_main

from shapely.geometry import Polygon
from shapely.ops import cascaded_union

from pystac import Catalog, CatalogType, Item, Asset, LabelItem, Collection
from tqdm import tqdm<jupyter_output><empty_output><jupyter_text>Google driveMount google drive folder with data to machine provided by google collab.<jupyter_code>import os

from google.colab import drive
from pathlib import Path

drive.mount('/content/gdrive')

GOOGLE_DRIVE_PATH = Path('/content/gdrive/My Drive')
WORKDIR = 'segmentation' # specify the path to folder which you are intend to work with in this notebook
WORKDIR_PATH = GOOGLE_DRIVE_PATH / WORKDIR

if not os.path.exists(WORKDIR_PATH):
          os.mkdir(WORKDIR_PATH)<jupyter_output>Go to this URL in a browser: https://accounts.google.com/o/oauth2/auth?client_id=947318989803-6bn6qk8qdgf4n4g3pfee6491hc0brc4i.apps.googleusercontent.com&redirect_uri=urn%3aietf%3awg%3aoauth%3a2.0%3aoob&response_type=code&scope=email%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdocs.test%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdrive%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdrive.photos.readonly%20https%3a%2f%2fwww.googleapis.com%2fauth%2fpeopleapi.readonly

Enter your authorization code:
··········
Mounted at /content/gdrive<jupyter_text>Download and unzip the dataset. Training data<jupyter_code>!curl https://drivendata-public-assets.s3.amazonaws.com/train_tier_1.tgz | tar -xz  -C "/content/gdrive/My Drive/segmentation/data"<jupyter_output><empty_output><jupyter_text>Testing data<jupyter_code>!curl https://drivendata-public-assets.s3.amazonaws.com/test.tgz | tar -xz  -C "/content/gdrive/My Drive/segmentation/data"<jupyter_output><empty_output><jupyter_text>Set the paths<jupyter_code>IMG_PATH = WORKDIR_PATH /'images-256'
MASK_PATH = WORKDIR_PATH /'masks-256'
IMG_PATH.mkdir(exist_ok=True)
MASK_PATH.mkdir(exist_ok=True)<jupyter_output><empty_output><jupyter_text>Load our training catalog<jupyter_code>train1_catalog = Catalog.from_file('/content/gdrive/My Drive/segmentation/data/train_tier_1/catalog.json')
cols = {cols.id:cols for cols in train1_catalog.get_children()}
cols<jupyter_output><empty_output><jupyter_text>Get a list of the possible areas(scenes) and ids<jupyter_code>areas = []
for c in cols:
  itms = [x for x in cols[c].get_all_items()]
  for i, id in enumerate(itms):
    if i % 2 == 0 and i+1 < len(itms):
      areas.append((c, itms[i].id, itms[i+1].id))
areas<jupyter_output><empty_output><jupyter_text>Split the STAC into tiles.In this step we want to split the STAC file into images and masks of maller size (256px) Define function for saving images<jupyter_code>def save_tile_img(tif_url, xyz, tile_size, save_path='', prefix=''):
  x,y,z = xyz
  tile, mask = rt_main.tile(tif_url, x,y,z, tilesize=tile_size)
  skimage.io.imsave(f'{save_path}/{prefix}{z}_{x}_{y}.png', np.moveaxis(tile,0,2), check_contrast=False)<jupyter_output><empty_output><jupyter_text>Define function for saving masks.<jupyter_code>def save_tile_mask(labels_poly, tile_poly, xyz, tile_size, save_path='', prefix='', greyscale=True):
  x,y,z = xyz
  tfm = from_bounds(*tile_poly.bounds, tile_size, tile_size) 
  
  cropped_polys = [poly for poly in labels_poly if poly.intersects(tile_poly)]
  cropped_polys_gdf = gpd.GeoDataFrame(geometry=cropped_polys, crs='epsg:4326')
  
  fbc_mask = sol.vector.mask.df_to_px_mask(df=cropped_polys_gdf,
                                         channels=['footprint', 'boundary', 'contact'],
                                         affine_obj=tfm, shape=(tile_size,tile_size),
                                         boundary_width=5, boundary_type='inner', contact_spacing=5, meters=True)
  
  if greyscale:
    skimage.io.imsave(f'{save_path}/{prefix}{z}_{x}_{y}_mask.png',fbc_mask[:,:,0], check_contrast=False) 
  else:
    skimage.io.imsave(f'{save_path}/{prefix}{z}_{x}_{y}_mask.png',fbc_mask, check_contrast=False)<jupyter_output><empty_output><jupyter_text>Define function for splitting the STAC for whole area<jupyter_code>def save_area_id_images(area, img_id, label_id, zoom_level = 19, tile_size = 256):

  # The item
  one_item = cols[area].get_item(id=img_id)

  # Load labels shapefile
  lab = cols[area].get_item(id=label_id)
  gdf = gpd.read_file(lab.make_asset_hrefs_absolute().assets['labels'].href)
  # get the geometries from the geodataframe
  all_polys = gdf.geometry

  # Get outlines as polygons
  polygon_geom = Polygon(one_item.to_dict()['geometry']['coordinates'][0])
  polygon = gpd.GeoDataFrame(index=[0], crs=gdf.crs, geometry=[polygon_geom])   

  # Tile at zoom_level
  polygon['geometry'].to_file(img_id+'.geojson', driver='GeoJSON')
  !cat {img_id}.geojson | supermercado burn {zoom_level} | mercantile shapes | fio collect > {img_id}{zoom_level}tiles.geojson

  # Load tiles
  tiles = gpd.read_file(f'{img_id}{zoom_level}tiles.geojson')

  # Add a convenience column
  tiles['xyz'] = tiles.id.apply(lambda x: x.lstrip('(,)').rstrip('(,)').split(','))
  tiles['xyz'] = [[int(q) for q in p] for p in tiles['xyz']]

  # IMG URL
  tif_url = one_item.assets['image'].href

  # Sometimes it's just ./id.tif - add full path (should maybe use make_asset_hrefs_absolute instead!!)
  if tif_url.startswith("./"):
    tif_url = '/'.join(one_item.to_dict()['links'][1]['href'].split("/")[:-1])+tif_url[1:]

  print("TIF URL:", tif_url)

  print("Number of tiles:", len(tiles))

  # Loop through tiles, downloading and saving
  for idx in range(len(tiles)):
    tile, mask = rt_main.tile(tif_url, *tiles.iloc[idx]['xyz'], tilesize=tile_size)

    tile_poly = tiles.iloc[idx]['geometry']

    # get affine transformation matrix for this tile using rasterio.transform.from_bounds: https://rasterio.readthedocs.io/en/stable/api/rasterio.transform.html#rasterio.transform.from_bounds
    tfm = from_bounds(*tile_poly.bounds, tile_size, tile_size) 

    # crop geometries to what overlaps our tile polygon bounds
    cropped_polys = [poly for poly in all_polys if poly.intersects(tile_poly)]
    cropped_polys_gdf = gpd.GeoDataFrame(geometry=cropped_polys, crs='epsg:4326')

    # burn a footprint/boundary/contact 3-channel mask with solaris: https://solaris.readthedocs.io/en/latest/tutorials/notebooks/api_masks_tutorial.html
    fbc_mask = sol.vector.mask.df_to_px_mask(df=cropped_polys_gdf,
                                            channels=['footprint', 'boundary', 'contact'],
                                            affine_obj=tfm, shape=(tile_size,tile_size),
                                            boundary_width=5, boundary_type='inner', contact_spacing=5, meters=True)

    save_tile_img(tif_url, tiles.iloc[idx]['xyz'], tile_size, save_path=IMG_PATH, prefix=f'{area}_{img_id}_{idx}_')
    save_tile_mask(all_polys, tile_poly, tiles.iloc[idx]['xyz'], tile_size, save_path=MASK_PATH,prefix=f'{area}_{img_id}_{idx}_')
    print("Saved", f'{area}_{img_id}_{idx}_')<jupyter_output><empty_output><jupyter_text>Actually create tiles<jupyter_code>for a in areas:
  save_area_id_images(area = a[0], img_id = a[1], label_id=a[2], zoom_level = 19, tile_size = 256)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1254722381224957, 'fraction_numerical': 0.016370921579851174, 'mean_word_length': 5.293948126801153, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 49, 'https://': 12, 'lorem ipsum': 0, 'www.': 5, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '577889', 'license': 'MIT', 'n_tokens_mistral': 3274, 'n_tokens_neox': 3101, 'n_words': 588, 'path': 'preprocessing.ipynb', 'repo_name': 'rpitonak/building-segmentation', 'size': '8735'}