UNET_G4G_2019_Parking.ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"UNET_G4G_2019_Parking.ipynb","provenance":[],"private_outputs":true,"collapsed_sections":[],"toc_visible":true,"machine_shape":"hm"},"kernelspec":{"name":"python3","display_name":"Python 3"},"accelerator":"GPU"},"cells":[{"cell_type":"markdown","metadata":{"id":"view-in-github","colab_type":"text"},"source":["<a href=\"https://colab.research.google.com/github/mjevans26/Satellite_ComputerVision/blob/master/UNET_G4G_2019_Parking.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"]},{"cell_type":"code","metadata":{"id":"esIMGVxhDI0f","colab_type":"code","colab":{}},"source":["#@title Copyright 2019 Google LLC. { display-mode: \"form\" }\n","# Licensed under the Apache License, Version 2.0 (the \"License\");\n","# you may not use this file except in compliance with the License.\n","# You may obtain a copy of the License at\n","#\n","# https://www.apache.org/licenses/LICENSE-2.0\n","#\n","# Unless required by applicable law or agreed to in writing, software\n","# distributed under the License is distributed on an \"AS IS\" BASIS,\n","# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n","# See the License for the specific language governing permissions and\n","# limitations under the License."],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"_SHAc5qbiR8l","colab_type":"text"},"source":["# Introduction\n","\n","This is a Google Colab notebook demonstrating the process used to export training, evaluation, and prediction data from Google Earth Engine used to develop a [Deeplab V3](https://arxiv.org/abs/1706.05587) convolutional neural network that delineates ground mounted solar arrays in [NAIP](https://www.fsa.usda.gov/programs-and-services/aerial-photography/imagery-programs/naip-imagery/).  Model training and predictions are accomplished in a separate notebook."]},{"cell_type":"markdown","metadata":{"id":"_MJ4kW1pEhwP","colab_type":"text"},"source":["# Setup software libraries\n","\n","Install needed libraries to the notebook VM.  Authenticate as necessary."]},{"cell_type":"code","metadata":{"id":"neIa46CpciXq","colab_type":"code","colab":{}},"source":["# Cloud authentication.\n","from google.colab import auth\n","auth.authenticate_user()"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"4D6ArFWrckmS","colab_type":"code","colab":{}},"source":["# Earth Engine install to notebook VM.\n","!pip install earthengine-api"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"jat01FEoUMqg","colab_type":"code","colab":{}},"source":["# Import, authenticate and initialize the Earth Engine library.\n","import ee\n","ee.Authenticate()\n","ee.Initialize()"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"8RnZzcYhcpsQ","colab_type":"code","colab":{}},"source":["# Tensorflow setup.\n","import tensorflow as tf\n","\n","tf.enable_eager_execution()\n","print(tf.__version__)\n","\n","%load_ext tensorboard"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"n1hFdpBQfyhN","colab_type":"code","colab":{}},"source":["# Folium setup.\n","import folium\n","print(folium.__version__)\n","\n","# Define the URL format used for Earth Engine generated map tiles.\n","EE_TILES = 'https://earthengine.googleapis.com/map/{mapid}/{{z}}/{{x}}/{{y}}?token={token}'"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"WjUgYcsAs9Ed","colab_type":"text"},"source":["##Mount Google Drive"]},{"cell_type":"code","metadata":{"id":"JKDKpX4FtQA1","colab_type":"code","colab":{}},"source":["# Attach specified google drive directory to this notebook\n","from google.colab import drive\n","drive.mount('/content/drive')"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"M6pVAfdDIJ-a","colab_type":"code","colab":{}},"source":["%cd '/content/drive/My Drive/repos/Satellite_ComputerVision'"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"iT8ycmzClYwf","colab_type":"text"},"source":["# Variables\n","\n","Declare the variables that will be in use throughout the notebook."]},{"cell_type":"markdown","metadata":{"id":"qKs6HuxOzjMl","colab_type":"text"},"source":["## Specify your Cloud Storage Bucket\n","You must have write access to a bucket to run this demo!  To run it read-only, use the demo bucket below, but note that writes to this bucket will not work."]},{"cell_type":"code","metadata":{"id":"obDDH1eDzsch","colab_type":"code","colab":{}},"source":["# This is read-only:\n","BUCKET = 'cvod-203614-mlengine'"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"wmfKLl9XcnGJ","colab_type":"text"},"source":["## Set other global variables"]},{"cell_type":"code","metadata":{"id":"psz7wJKalaoj","colab_type":"code","colab":{}},"source":["# Specify names locations for outputs in Cloud Storage. \n","FOLDER = 'LI_parking'\n","PRED_BASE = 'data/predict'\n","TRAINING_BASE = 'data/training'\n","EVAL_BASE = 'data/eval'\n","MODEL_BASE = 'models'\n","log_dir = 'drive/My Drive/Tensorflow/models/UNET256'\n","\n","# Specify inputs (Landsat bands) to the model and the response variable.\n","opticalBands = ['R', 'G', 'B']\n","thermalBands = ['B8', 'B11', 'B12']\n","pcaBands = ['pc1', 'pc2', 'pc3']\n","BANDS = opticalBands# + thermalBands# + pcaBands\n","RESPONSE = 'landcover'\n","FEATURES = BANDS + [RESPONSE]\n","\n","# Specify the size and shape of patches expected by the model.\n","KERNEL_SIZE = 512\n","KERNEL_SHAPE = [KERNEL_SIZE, KERNEL_SIZE]\n","COLUMNS = [\n","  tf.io.FixedLenFeature(shape=KERNEL_SHAPE, dtype=tf.float32) for k in FEATURES\n","]\n","FEATURES_DICT = dict(zip(FEATURES, COLUMNS))\n","\n","# Sizes of the training and evaluation datasets.\n","TRAIN_SIZE = 8000\n","EVAL_SIZE = 5000\n","\n","# Specify model training parameters.\n","BATCH_SIZE = 16\n","EPOCHS = 20\n","BUFFER_SIZE = 8000\n","OPTIMIZER = tf.keras.optimizers.Adam(learning_rate=0.001, beta_1=0.9, beta_2=0.999)\n","LOSS = 'binary_crossentropy'\n","METRICS = [tf.keras.metrics.categorical_accuracy, tf.keras.metrics.MeanIoU(num_classes=2)]"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"hgoDc7Hilfc4","colab_type":"text"},"source":["# Imagery\n","\n","Gather and setup the imagery to use for inputs (predictors).  This is a three-year, cloud-free, Landsat 8 composite.  Display it in the notebook for a sanity check."]},{"cell_type":"code","metadata":{"id":"-IlgXu-vcUEY","colab_type":"code","colab":{}},"source":["# Use Landsat 8 surface reflectance data.\n","NAIP = ee.ImageCollection(\"USDA/NAIP/DOQQ\")\n","towns = ee.FeatureCollection(\"users/defendersofwildlifeGIS/LongIsland/towns\")\n","\n","begin = '2017-01-01'\n","end = '2017-12-30'\n","\n","# The image input data is a cloud-masked median composite.\n","image = NAIP.filterDate(begin, end)\\\n",".filterBounds(towns)\\\n",".filterDate(begin, end)\\\n",".median()\\\n",".select(BANDS)\\\n",".clip(towns)\n","\n","# Use folium to visualize the imagery.\n","mapid = image.getMapId({'bands': ['R', 'G', 'B'], 'min': 0, 'max': 256})\n","map = folium.Map(location=[40.8175, -73.195])\n","folium.TileLayer(\n","    tiles=EE_TILES.format(**mapid),\n","    attr='Google Earth Engine',\n","    overlay=True,\n","    name='median composite',\n","  ).add_to(map)\n","\n","map.add_child(folium.LayerControl())\n","map"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"gHznnctkJsZJ","colab_type":"text"},"source":["Prepare the response (what we want to predict).  This is impervious surface area (in fraction of a pixel) from the 2016 NLCD dataset.  Display to check."]},{"cell_type":"code","metadata":{"id":"5Wxz9BPYHBwh","colab_type":"code","colab":{}},"source":["def set_landcover(ft):\n","  return ft.set('label', 1)\n","\n","nassauParkingFootprints = ee.FeatureCollection(\"users/defendersofwildlifeGIS/LongIsland/NassauParking\")\n","suffolkParkingFootprints = ee.FeatureCollection('users/defendersofwildlifeGIS/LongIsland/SuffolkParking')\n","parkingFootprints = nassauParkingFootprints.merge(suffolkParkingFootprints)\n","parking = parkingFootprints.map(set_landcover)\n","blankimg = ee.Image.constant(0)\n","parking_footprint = parking.reduceToImage(['label'], ee.Reducer.first())\n","labelimg = blankimg.where(parking_footprint, parking_footprint).rename('landcover')\n","\n","mapid = labelimg.getMapId({'bands': 'landcover', 'min':0, 'max': 1})\n","print(mapid)\n","map = folium.Map(location = [40.8175, -73.195])\n","folium.TileLayer(\n","    tiles = EE_TILES.format(**mapid),\n","    attr='Google Earth Engine',\n","    overlay = True,\n","    name = 'parking lots',\n",").add_to(map)\n","map.add_child(folium.LayerControl())\n","map"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"CTS7_ZzPDhhg","colab_type":"text"},"source":["Stack the 2D images (Landsat composite and NLCD impervious surface) to create a single image from which samples can be taken.  Convert the image into an array image in which each pixel stores 256x256 patches of pixels for each band.  This is a key step that bears emphasis: to export training patches, convert a multi-band image to [an array image](https://developers.google.com/earth-engine/arrays_array_images#array-images) using [`neighborhoodToArray()`](https://developers.google.com/earth-engine/api_docs#eeimageneighborhoodtoarray), then sample the image at points."]},{"cell_type":"code","metadata":{"id":"eGHYsdAOipa4","colab_type":"code","colab":{}},"source":["featureStack = ee.Image.cat([\n","  image.select(BANDS),\n","  labelimg.select(RESPONSE)\n","]).float()\n","\n","print(featureStack.bandNames().getInfo())\n","\n","list = ee.List.repeat(1, KERNEL_SIZE)\n","lists = ee.List.repeat(list, KERNEL_SIZE)\n","kernel = ee.Kernel.fixed(KERNEL_SIZE, KERNEL_SIZE, lists)\n","\n","arrays = featureStack.neighborhoodToArray(kernel)"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"F4djSxBRG2el","colab_type":"text"},"source":["Use some pre-made geometries to sample the stack in strategic locations.  Specifically, these are hand-made polygons in which to take the 256x256 samples.  Display the sampling polygons on a map, red for training polygons, blue for evaluation."]},{"cell_type":"code","metadata":{"id":"ure_WaD0itQY","colab_type":"code","colab":{}},"source":["import re\n","towns = towns.randomColumn('random', 52.0)\n","townList = ee.List(towns.aggregate_array('TOWN')).distinct()\n","townList = [townList.get(town).getInfo() for town in range(townList.size().getInfo())]\n","townList = [town for town in townList if not re.search(r\"City|Water|Indian\", town)]\n","trainList = townList[0:(len(townList)//10) * 8]\n","evalList = townList[(len(townList)//10) * 8:]\n","\n","trainFilter = ee.Filter.inList('TOWN', ee.List(trainList))\n","evalFilter = ee.Filter.inList(\"TOWN\", ee.List(evalList))\n","\n","trainingPolys = towns.filter(trainFilter)\n","print('training size', len(trainList))\n","\n","evalPolys = towns.filter(evalFilter)\n","print('eval size', len(evalList))\n","\n","polyImage = ee.Image(0).byte().paint(trainingPolys, 1).paint(evalPolys, 2)\n","polyImage = polyImage.updateMask(polyImage)\n","\n","mapid = polyImage.getMapId({'min': 1, 'max': 2, 'palette': ['red', 'blue']})\n","map = folium.Map(location=[40.8175, -73.195], zoom_start=8)\n","folium.TileLayer(\n","    tiles=EE_TILES.format(**mapid),\n","    attr='Google Earth Engine',\n","    overlay=True,\n","    name='training polygons',\n","  ).add_to(map)\n","map.add_child(folium.LayerControl())\n","map"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"ZV890gPHeZqz","colab_type":"text"},"source":["# Sampling\n","\n","The mapped data look reasonable so take a sample from each polygon and merge the results into a single export.  The key step is sampling the array image at points, to get all the pixels in a 256x256 neighborhood at each point.  It's worth noting that to build the training and testing data for the FCNN, you export a single TFRecord file that contains patches of pixel values in each record.  You do NOT need to export each training/testing patch to a different image.  Since each record potentially contains a lot of data (especially with big patches or many input bands), some manual sharding of the computation is necessary to avoid the `computed value too large` error.  Specifically, the following code takes multiple (smaller) samples within each geometry, merging the results to get a single export."]},{"cell_type":"code","metadata":{"id":"FyRpvwENxE-A","colab_type":"code","cellView":"code","colab":{}},"source":["#@title Don't run\n","# Convert the feature collections to lists for iteration.\n","#trainingPolysList = trainingPolys.toList(trainingPolys.size())\n","#evalPolysList = trainingPolys.toList(trainingPolys.size())\n","\n","# These numbers determined experimentally.\n","n = 100 # Number of shards in each town.\n","N = 1000 # Total sample size in each town.\n","\n","for town in trainList:\n","  geomSample = ee.FeatureCollection([])\n","  for i in range (n):\n","    sample = arrays.sample(\n","      region = trainingPolys.filterMetadata('TOWN', 'equals', town),\n","      scale = 1,\n","      numPixels = N/n,\n","      seed = i,\n","      tileScale = 8\n","    )\n","    geomSample = geomSample.merge(sample)\n","    \n","  desc = 'DeepLab_' + str(KERNEL_SIZE) + '_NAIP_' + town\n","  task = ee.batch.Export.table.toCloudStorage(\n","      collection = geomSample,\n","      description = desc,\n","      bucket = BUCKET,\n","      fileNamePrefix = FOLDER + '/' + TRAINING_BASE + '/' + desc,\n","      fileFormat = 'TFRecord',\n","      selectors = BANDS + [RESPONSE]\n","  )\n","  task.start()\n","    \n","for town in evalList:\n","  geomSample = ee.FeatureCollection([])\n","  for i in range(n):\n","    sample = arrays.sample(\n","        region = evalPolys.filterMetadata('TOWN', 'equals', town),\n","        scale = 1,\n","        numPixels = N/n,\n","        seed = i,\n","        tileScale = 8\n","    )\n","    geomSample = geomSample.merge(sample)\n","    \n","  desc = 'DeepLab_' + str(KERNEL_SIZE) + 'NAIP_' + town\n","  task = ee.batch.Export.table.toCloudStorage(\n","      collection = geomSample,\n","      description = desc,\n","      bucket = BUCKET,\n","      fileNamePrefix = FOLDER + '/' + EVAL_BASE + '/' + desc,\n","      fileFormat = 'TFRecord',\n","      selectors = BANDS + [RESPONSE]\n","  )\n","  task.start()\n","                            \n","#Export all the training data (in many pieces), with one task \n","#per geometry.\n","# for g in range(trainingPolys.size().getInfo()):\n","#   geomSample = ee.FeatureCollection([])\n","#   for i in range(n):\n","#     sample = arrays.sample(\n","#       region = ee.Feature(trainingPolysList.get(g)).geometry(), \n","#       scale = 30, \n","#       numPixels = N / n, # Size of the shard.\n","#       seed = i,\n","#       tileScale = 8\n","#     )\n","#     geomSample = geomSample.merge(sample)\n","  \n","#   desc = TRAINING_BASE + '_g' + str(g)\n","#   task = ee.batch.Export.table.toCloudStorage(\n","#     collection = geomSample,\n","#     description = desc, \n","#     bucket = BUCKET, \n","#     fileNamePrefix = FOLDER + '/' + desc,\n","#     fileFormat = 'TFRecord',\n","#     selectors = BANDS + [RESPONSE]\n","#   )\n","#   task.start()\n","\n","# # Export all the evaluation data.\n","# for g in range(evalPolys.size().getInfo()):\n","#   geomSample = ee.FeatureCollection([])\n","#   for i in range(n):\n","#     sample = arrays.sample(\n","#       region = ee.Feature(evalPolysList.get(g)).geometry(), \n","#       scale = 30, \n","#       numPixels = N / n,\n","#       seed = i,\n","#       tileScale = 8\n","#     )\n","#     geomSample = geomSample.merge(sample)\n","  \n","#   desc = EVAL_BASE + '_g' + str(g)\n","#   task = ee.batch.Export.table.toCloudStorage(\n","#     collection = geomSample,\n","#     description = desc, \n","#     bucket = BUCKET, \n","#     fileNamePrefix = FOLDER + '/' + desc,\n","#     fileFormat = 'TFRecord',\n","#     selectors = BANDS + [RESPONSE]\n","#   )\n","#   task.start()"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"dk51-l7MH2Sa","colab_type":"text"},"source":["##Preprocessing\n","Define functions that apply random manipulations to our training data"]},{"cell_type":"code","metadata":{"id":"ajyp48-vINuy","colab_type":"code","colab":{}},"source":["def augColor(x):\n","    \"\"\"Color augmentation\n","\n","    Args:\n","        x: Image\n","\n","    Returns:\n","        Augmented image\n","    \"\"\"\n","    x = tf.image.random_hue(x, 0.08)\n","    x = tf.image.random_saturation(x, 0.6, 1.6)\n","    x = tf.image.random_brightness(x, 0.05)\n","    x = tf.image.random_contrast(x, 0.7, 1.3)\n","    return x\n","  \n","  \n","def augImg(img):\n","  outDims = tf.shape(img)[0:1]\n","  x = tf.image.random_flip_left_right(img)\n","  x = tf.image.random_flip_up_down(x)\n","  x = rotated = tf.image.rot90(x, tf.random_uniform(shape=[], minval=0, maxval=4, dtype=tf.int32))\n","  #x = zoom(x, outDims)\n","  #since were gonna map_fn this on a 4d image, output must be 3d, so squeeze the artificial 'sample' dimension\n","  return tf.squeeze(x)\n","\n","def preprocess(img, labels):\n","  dims = tf.shape(img)\n","  #need to combine labels and bands for morphological transformations\n","  comb = tf.concat([img, tf.expand_dims(labels, axis = 2)], axis = 2)\n","  aug = aug_img(comb)\n","  #aug = tf.map_fn(fn = aug_img, elems = comb)\n","  labels = tf.squeeze(aug[:, :, -1])\n","  band_stack = color(aug[:, :, 0:dims[2]])\n","  return band_stack, labels"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"rWXrvBE4607G","colab_type":"text"},"source":["# Training data\n","\n","Load the data exported from Earth Engine into a `tf.data.Dataset`.  The following are helper functions for that."]},{"cell_type":"code","metadata":{"id":"WWZ0UXCVMyJP","colab_type":"code","colab":{}},"source":["def parse_tfrecord(example_proto):\n","  \"\"\"The parsing function.\n","  Read a serialized example into the structure defined by FEATURES_DICT.\n","  Args:\n","    example_proto: a serialized Example.\n","  Returns: \n","    A dictionary of tensors, keyed by feature name.\n","  \"\"\"\n","  return tf.io.parse_single_example(example_proto, FEATURES_DICT)\n","\n","\n","def to_tuple(inputs):\n","  \"\"\"Function to convert a dictionary of tensors to a tuple of (inputs, outputs).\n","  Turn the tensors returned by parse_tfrecord into a stack in HWC shape.\n","  Args:\n","    inputs: A dictionary of tensors, keyed by feature name.\n","  Returns: \n","    A dtuple of (inputs, outputs).\n","  \"\"\"\n","  inputsList = [inputs.get(key) for key in FEATURES]\n","  stacked = tf.stack(inputsList, axis=0)\n","  # Convert from CHW to HWC\n","  stacked = tf.transpose(stacked, [1, 2, 0])\n","  stacked = augImg(stacked)\n","  return stacked[:,:,:len(BANDS)], stacked[:,:,len(BANDS):]\n","\n","\n","def get_dataset(pattern):\n","  \"\"\"Function to read, parse and format to tuple a set of input tfrecord files.\n","  Get all the files matching the pattern, parse and convert to tuple.\n","  Args:\n","    pattern: A file pattern to match in a Cloud Storage bucket.\n","  Returns: \n","    A tf.data.Dataset\n","  \"\"\"\n","  glob = tf.gfile.Glob(pattern)\n","  dataset = tf.data.TFRecordDataset(glob, compression_type='GZIP')\n","  dataset = dataset.map(parse_tfrecord, num_parallel_calls=5)\n","  dataset = dataset.map(to_tuple, num_parallel_calls=5)\n","  return dataset"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"Xg1fa18336D2","colab_type":"text"},"source":["Use the helpers to read in the training dataset.  Print the first record to check."]},{"cell_type":"code","metadata":{"id":"rm0qRF0fAYcC","colab_type":"code","colab":{}},"source":["def get_training_dataset():\n","\t\"\"\"Get the preprocessed training dataset\n","  Returns: \n","    A tf.data.Dataset of training data.\n","  \"\"\"\n","\tglob = 'gs://' + BUCKET + '/' + FOLDER + '/' + TRAINING_BASE + '/*'\n","\tprint(glob)\n","\tdataset = get_dataset(glob)\n","\tdataset = dataset.shuffle(8000).batch(BATCH_SIZE).repeat()\n","\treturn dataset\n","\n","training = get_training_dataset()\n","\n","print(iter(training.take(1)).next())"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"7CRGG26bYWQZ","colab_type":"code","colab":{}},"source":["print(iter(training.take(1)).next())"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"j-cQO5RL6vob","colab_type":"text"},"source":["# Evaluation data\n","\n","Now do the same thing to get an evaluation dataset.  Note that unlike the training dataset, the evaluation dataset has a batch size of 1, is not repeated and is not shuffled."]},{"cell_type":"code","metadata":{"id":"ieKTCGiJ6xzo","colab_type":"code","colab":{}},"source":["def get_eval_dataset():\n","\t\"\"\"Get the preprocessed evaluation dataset\n","  Returns: \n","    A tf.data.Dataset of evaluation data.\n","  \"\"\"\n","\tglob = 'gs://' + BUCKET + '/' + FOLDER + '/' + EVAL_BASE + '/*'\n","\tprint(glob)\n","\tdataset = get_dataset(glob)\n","\tdataset = dataset.batch(1).repeat()\n","\treturn dataset\n","\n","evaluation = get_eval_dataset()"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"keoalUvBbSkh","colab_type":"code","colab":{}},"source":["print(iter(evaluation.take(1)).next())"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"9JIE7Yl87lgU","colab_type":"text"},"source":["# Model\n","\n","Here we use the Keras implementation of the U-Net model as found [in the TensorFlow examples](https://github.com/tensorflow/models/blob/master/samples/outreach/blogs/segmentation_blogpost/image_segmentation.ipynb).  The U-Net model takes 256x256 pixel patches as input and outputs per-pixel class probability, label or a continuous output.  We can implement the model essentially unmodified, but will use mean squared error loss on the sigmoidal output since we are treating this as a regression problem, rather than a classification problem.  Since impervious surface fraction is constrained to [0,1], with many values close to zero or one, a saturating activation function is suitable here."]},{"cell_type":"markdown","metadata":{"id":"Xh2EZyyPu84H","colab_type":"text"},"source":["##Metrics"]},{"cell_type":"code","metadata":{"id":"mISCOXUHu7G_","colab_type":"code","colab":{}},"source":["def weighted_bce(y_true, y_pred):\n","  bce = tf.nn.weighted_cross_entropy_with_logits(labels = y_true, logits = y_pred, pos_weight = 20)\n","  return tf.reduce_mean(bce)\n","\n","def iou(true, pred):\n","\n","    intersection = true * pred\n","\n","    notTrue = 1 - true\n","    union = true + (notTrue * pred)\n","\n","    return tf.reduce_sum(intersection)/tf.reduce_sum(union)\n","\n"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"wsnnnz56yS3l","colab_type":"code","colab":{}},"source":["from tensorflow.python.keras import layers\n","from tensorflow.python.keras import losses\n","from tensorflow.python.keras import models\n","from tensorflow.python.keras import metrics\n","from tensorflow.python.keras import optimizers\n","\n","def conv_block(input_tensor, num_filters):\n","\tencoder = layers.Conv2D(num_filters, (3, 3), padding='same')(input_tensor)\n","\tencoder = layers.BatchNormalization()(encoder)\n","\tencoder = layers.Activation('relu')(encoder)\n","\tencoder = layers.Conv2D(num_filters, (3, 3), padding='same')(encoder)\n","\tencoder = layers.BatchNormalization()(encoder)\n","\tencoder = layers.Activation('relu')(encoder)\n","\treturn encoder\n","\n","def encoder_block(input_tensor, num_filters):\n","\tencoder = conv_block(input_tensor, num_filters)\n","\tencoder_pool = layers.MaxPooling2D((2, 2), strides=(2, 2))(encoder)\n","\treturn encoder_pool, encoder\n","\n","def decoder_block(input_tensor, concat_tensor, num_filters):\n","\tdecoder = layers.Conv2DTranspose(num_filters, (2, 2), strides=(2, 2), padding='same')(input_tensor)\n","\tdecoder = layers.concatenate([concat_tensor, decoder], axis=-1)\n","\tdecoder = layers.BatchNormalization()(decoder)\n","\tdecoder = layers.Activation('relu')(decoder)\n","\tdecoder = layers.Conv2D(num_filters, (3, 3), padding='same')(decoder)\n","\tdecoder = layers.BatchNormalization()(decoder)\n","\tdecoder = layers.Activation('relu')(decoder)\n","\tdecoder = layers.Conv2D(num_filters, (3, 3), padding='same')(decoder)\n","\tdecoder = layers.BatchNormalization()(decoder)\n","\tdecoder = layers.Activation('relu')(decoder)\n","\treturn decoder\n","\n","def get_model():\n","\tinputs = layers.Input(shape=[None, None, len(BANDS)]) # 256\n","\tencoder0_pool, encoder0 = encoder_block(inputs, 32) # 128\n","\tencoder1_pool, encoder1 = encoder_block(encoder0_pool, 64) # 64\n","\tencoder2_pool, encoder2 = encoder_block(encoder1_pool, 128) # 32\n","\tencoder3_pool, encoder3 = encoder_block(encoder2_pool, 256) # 16\n","\tencoder4_pool, encoder4 = encoder_block(encoder3_pool, 512) # 8\n","\tcenter = conv_block(encoder4_pool, 1024) # center\n","\tdecoder4 = decoder_block(center, encoder4, 512) # 16\n","\tdecoder3 = decoder_block(decoder4, encoder3, 256) # 32\n","\tdecoder2 = decoder_block(decoder3, encoder2, 128) # 64\n","\tdecoder1 = decoder_block(decoder2, encoder1, 64) # 128\n","\tdecoder0 = decoder_block(decoder1, encoder0, 32) # 256\n","\toutputs = layers.Conv2D(1, (1, 1), activation='sigmoid')(decoder0)\n","\n","\tmodel = models.Model(inputs=[inputs], outputs=[outputs])\n","\n","\tmodel.compile(\n","\t\toptimizer=OPTIMIZER, \n","    loss = weighted_bce,\n","\t\t#loss=losses.get(LOSS),\n","\t\tmetrics=[metrics.get(metric) for metric in METRICS])\n","\n","\treturn model\n","\n","\n","log_dir = 'drive/My Drive/Tensorflow/models/UNET256'\n","\n","checkpoint = tf.keras.callbacks.ModelCheckpoint(\n","    log_dir+'best_weights.hdf5',\n","    monitor='val_mean_io_u',\n","    verbose=1,\n","    save_best_only=True,\n","    mode='max'\n","    )"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"uu_E7OTDBCoS","colab_type":"text"},"source":["# Training the model\n","\n","You train a Keras model by calling `.fit()` on it.  Here we're going to train for 10 epochs, which is suitable for demonstration purposes.  For production use, you probably want to optimize this parameter, for example through [hyperparamter tuning](https://cloud.google.com/ml-engine/docs/tensorflow/using-hyperparameter-tuning)."]},{"cell_type":"code","metadata":{"id":"NzzaWxOhSxBy","colab_type":"code","colab":{}},"source":["m = get_model()\n","\n","\n","m.fit(\n","    x=training, \n","    epochs=EPOCHS, \n","    steps_per_epoch=int(TRAIN_SIZE / BATCH_SIZE), \n","    validation_data=evaluation,\n","    validation_steps=EVAL_SIZE/BATCH_SIZE,\n","    callbacks = [checkpoint]\n","    )\n","\n","m.save('drive/My Drive/Tensorflow/models/UNET256/UNET256.h5')\n","\n","#!gsutil cp best_weights.hdf5 gs://cvod-203614-mlengine/NC_solar/models/UNET256/best_weights.hdf5\n","#!gsutil cp UNET256.h5 gs://cvod-203614-mlengine/NC_solar/models/UNET256/UNET256.h5"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"U2XrwZHp66j4","colab_type":"text"},"source":["Note that the notebook VM is sometimes not heavy-duty enough to get through a whole training job, especially if you have a large buffer size or a large number of epochs.  You can still use this notebook for training, but may need to set up an alternative VM ([learn more](https://research.google.com/colaboratory/local-runtimes.html)) for production use.  Alternatively, you can package your code for running large training jobs on Google's AI Platform [as described here](https://cloud.google.com/ml-engine/docs/tensorflow/trainer-considerations).  The following code loads a pre-trained model, which you can use for predictions right away."]},{"cell_type":"markdown","metadata":{"id":"zvIqqpNXqJSE","colab_type":"text"},"source":["##Load model and resume training"]},{"cell_type":"code","metadata":{"id":"q0xgBhsaqInV","colab_type":"code","colab":{}},"source":["#bring in the architecture and best weights from GCS\n","m = models.load_model('drive/My Drive/Tensorflow/models/UNET256/UNET256.h5', custom_objects={'weighted_bce': weighted_bce})\n","m.load_weights('drive/My Drive/Tensorflow/models/UNET256/best_weights.hdf5') \n","\n","#lets see where were at\n","evalMetrics = m.evaluate(x=evaluation, steps = EVAL_SIZE, verbose = 1)"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"xlsFciElxOUA","colab_type":"code","colab":{}},"source":["#set the monitored value (val_mean_io_u) to current evaluation output\n","checkpoint = tf.keras.callbacks.ModelCheckpoint(\n","    log_dir+'best_weights.hdf5',\n","    monitor='val_mean_io_u',\n","    verbose=1,\n","    save_best_only=True,\n","    mode='max'\n","    )\n","\n","checkpoint.best = evalMetrics[2]\n","print(checkpoint.__dict__)\n","print(checkpoint.best)"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"7eq0aLlw864A","colab_type":"text"},"source":["## Set up tensorboard"]},{"cell_type":"code","metadata":{"id":"PA2gJENE8-J1","colab_type":"code","colab":{}},"source":["tensorboard = tf.keras.callbacks.TensorBoard(log_dir= 'drive/My Drive/Tensorflow/models/UNET256')"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"Ty8wCxDtqWBM","colab_type":"code","colab":{}},"source":["#Now keep training!\n","m.fit(\n","    x=training, \n","    epochs= 10, \n","    steps_per_epoch=int(TRAIN_SIZE / BATCH_SIZE), \n","    validation_data=evaluation,\n","    validation_steps=EVAL_SIZE/BATCH_SIZE,\n","    callbacks = [checkpoint, tensorboard]\n","    )\n","#m.save('drive/My Drive/Tensorflow/models/UNET256/UNET256.h5')"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"tyhWcGHJ82e8","colab_type":"code","colab":{}},"source":["m.save('drive/My Drive/Tensorflow/models/UNET256/UNET256.h5')"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"i9OM5BiS1xYQ","colab_type":"code","colab":{}},"source":["%tensorboard --logdir 'drive/My Drive/Tensorflow/models/UNET256'"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"-RJpNfEUS1qp","colab_type":"code","colab":{}},"source":["# Load a trained model. 50 epochs. 25 hours. Final RMSE ~0.08.\n","MODEL_DIR = BUCKET + '/' + FOLDER + '/' + 'models/UNET256'\n","m = tf.contrib.saved_model.load_keras_model(MODEL_DIR)\n","m.summary()"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"J1ySNup0xCqN","colab_type":"text"},"source":["# Prediction\n","\n","The prediction pipeline is:\n","\n","1.  Export imagery on which to do predictions from Earth Engine in TFRecord format to a Cloud Storge bucket.\n","2.  Use the trained model to make the predictions.\n","3.  Write the predictions to a TFRecord file in a Cloud Storage.\n","4.  Upload the predictions TFRecord file to Earth Engine.\n","\n","The following functions handle this process.  It's useful to separate the export from the predictions so that you can experiment with different models without running the export every time."]},{"cell_type":"code","metadata":{"id":"lv6nb0ShH4_T","colab_type":"code","colab":{}},"source":["#Inspect the prediction outputs\n","predictions = m.predict(evaluation, steps=1, verbose=1)\n","for prediction in predictions:\n","  print(predictions)"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"M3WDAa-RUpXP","colab_type":"code","colab":{}},"source":["def doExport(image, out_image_base, kernel_buffer, region):\n","  \"\"\"Run the image export task.  Block until complete.\n","  \"\"\"\n","  task = ee.batch.Export.image.toCloudStorage(\n","    image = image.select(BANDS+[RESPONSE]), \n","    description = out_image_base, \n","    bucket = BUCKET, \n","    fileNamePrefix = FOLDER + '/' + PRED_BASE + '/' + out_image_base, \n","    region = region.getInfo()['coordinates'], \n","    scale = 1, \n","    fileFormat = 'TFRecord', \n","    maxPixels = 1e10,\n","    formatOptions = { \n","      'patchDimensions': KERNEL_SHAPE,\n","      'kernelSize': kernel_buffer,\n","      'compressed': True,\n","      'maxFileSize': 104857600\n","    }\n","  )\n","  task.start()\n","\n","  # Block until the task completes.\n","  print('Running image export to Cloud Storage...')\n","  import time\n","  while task.active():\n","    time.sleep(30)\n","\n","  # Error condition\n","  if task.status()['state'] != 'COMPLETED':\n","    print('Error with image export.')\n","  else:\n","    print('Image export completed.')"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"zb_9_FflygVw","colab_type":"code","colab":{}},"source":["def doPrediction(out_image_base, user_folder, kernel_buffer, region):\n","  \"\"\"Perform inference on exported imagery, upload to Earth Engine.\n","  \"\"\"\n","\n","  print('Looking for TFRecord files...')\n","  \n","  # Get a list of all the files in the output bucket.\n","  filesList = !gsutil ls 'gs://'{BUCKET}'/'{FOLDER}'/'{PRED_BASE}\n","  # Get only the files generated by the image export.\n","  exportFilesList = [s for s in filesList if out_image_base in s]\n","\n","  # Get the list of image files and the JSON mixer file.\n","  imageFilesList = []\n","  jsonFile = None\n","  for f in exportFilesList:\n","    if f.endswith('.tfrecord.gz'):\n","      imageFilesList.append(f)\n","    elif f.endswith('.json'):\n","      jsonFile = f\n","\n","  # Make sure the files are in the right order.\n","  imageFilesList.sort()\n","\n","  from pprint import pprint\n","  pprint(imageFilesList)\n","  print(jsonFile)\n","  \n","  import json\n","  # Load the contents of the mixer file to a JSON object.\n","  jsonText = !gsutil cat {jsonFile}\n","  # Get a single string w/ newlines from the IPython.utils.text.SList\n","  mixer = json.loads(jsonText.nlstr)\n","  pprint(mixer)\n","  patches = mixer['totalPatches']\n","  \n","  # Get set up for prediction.\n","  x_buffer = int(kernel_buffer[0] / 2)\n","  y_buffer = int(kernel_buffer[1] / 2)\n","\n","  buffered_shape = [\n","      KERNEL_SHAPE[0] + kernel_buffer[0],\n","      KERNEL_SHAPE[1] + kernel_buffer[1]]\n","\n","  imageColumns = [\n","    tf.FixedLenFeature(shape=buffered_shape, dtype=tf.float32) \n","      for k in BANDS\n","  ]\n","\n","  imageFeaturesDict = dict(zip(BANDS, imageColumns))\n","\n","  def parse_image(example_proto):\n","    return tf.parse_single_example(example_proto, imageFeaturesDict)\n","\n","  def toTupleImage(dict):\n","    inputsList = [dict.get(key) for key in BANDS]\n","    stacked = tf.stack(inputsList, axis=0)\n","    stacked = tf.transpose(stacked, [1, 2, 0])\n","    return stacked\n","  \n","   # Create a dataset from the TFRecord file(s) in Cloud Storage.\n","  imageDataset = tf.data.TFRecordDataset(imageFilesList, compression_type='GZIP')\n","  imageDataset = imageDataset.map(parse_image, num_parallel_calls=5)\n","  imageDataset = imageDataset.map(toTupleImage).batch(1)\n","  \n","  # Perform inference.\n","  print('Running predictions...')\n","  predictions = m.predict(imageDataset, steps=patches, verbose=1)\n","  # print(predictions[0])\n","\n","  print('Writing predictions...')\n","  out_image_file = 'gs://' + BUCKET + '/' + FOLDER + '/' + PRED_BASE + '/outputs/' + out_image_base + '.TFRecord'\n","  writer = tf.python_io.TFRecordWriter(out_image_file)\n","  patches = 0\n","  for predictionPatch in predictions:\n","    print('Writing patch ' + str(patches) + '...')\n","    predictionPatch = predictionPatch[\n","        x_buffer:x_buffer+KERNEL_SIZE, y_buffer:y_buffer+KERNEL_SIZE]\n","\n","    # Create an example.\n","    example = tf.train.Example(\n","      features=tf.train.Features(\n","        feature={\n","          'probability': tf.train.Feature(\n","              float_list=tf.train.FloatList(\n","                  value=predictionPatch.flatten()))\n","        }\n","      )\n","    )\n","    # Write the example.\n","    writer.write(example.SerializeToString())\n","    patches += 1\n","\n","  writer.close()\n","\n","  # Start the upload.\n","  out_image_asset = user_folder + '/' + out_image_base\n","  !earthengine upload image --asset_id={out_image_asset} {out_image_file} {jsonFile}"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"LZqlymOehnQO","colab_type":"text"},"source":["Now there's all the code needed to run the prediction pipeline, all that remains is to specify the output region in which to do the prediction, the names of the output files, where to put them, and the shape of the outputs.  In terms of the shape, the model is trained on 256x256 patches, but can work (in theory) on any patch that's big enough with even dimensions ([reference](https://www.cv-foundation.org/openaccess/content_cvpr_2015/papers/Long_Fully_Convolutional_Networks_2015_CVPR_paper.pdf)).  Because of tile boundary artifacts, give the model slightly larger patches for prediction, then clip out the middle 256x256 patch.  This is controlled with a kernel buffer, half the size of which will extend beyond the kernel buffer.  For example, specifying a 128x128 kernel will append 64 pixels on each side of the patch, to ensure that the pixels in the output are taken from inputs completely covered by the kernel.  "]},{"cell_type":"code","metadata":{"id":"FPANwc7B1-TS","colab_type":"code","colab":{}},"source":["# This has a read-only asset in it:\n","user_folder = 'users/defendersofwildlifeGIS'\n","\n","# Base file name to use for TFRecord files and assets.\n","li_image_base = 'li_parking_deeplab512Pred'\n","# Half this will extend on the sides of each patch.\n","li_kernel_buffer = [256, 256]\n","# Huntington\n","li_region = ee.Feature(towns.filterMetadata(\"TOWN\", 'equals', 'Huntington').first()).geometry()\n","print(li_region.area().getInfo())"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"lLNEOLkXWvSi","colab_type":"code","cellView":"both","colab":{}},"source":["#@title Don't run\n","# Run the export.\n","doExport(featureStack, li_image_base, li_kernel_buffer, li_region)"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"KxACnxKFrQ_J","colab_type":"code","cellView":"both","colab":{}},"source":["#@title Don't run\n","# Run the prediction.\n","doPrediction(nc_image_base, user_folder, nc_kernel_buffer, nc_region)"],"execution_count":0,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"uj_G9OZ1xH6K","colab_type":"text"},"source":["# Display the output\n","\n","One the data has been exported, the model has made predictions and the predictions have been written to a file, and the image imported to Earth Engine, it's possible to display the resultant Earth Engine asset.  Here, display the impervious area predictions over Beijing, China."]},{"cell_type":"code","metadata":{"id":"Jgco6HJ4R5p2","colab_type":"code","colab":{}},"source":["out_image = ee.Image(user_folder + '/' + bj_image_base)\n","mapid = out_image.getMapId({'min': 0, 'max': 1})\n","map = folium.Map(location=[39.898, 116.5097])\n","folium.TileLayer(\n","    tiles=EE_TILES.format(**mapid),\n","    attr='Google Earth Engine',\n","    overlay=True,\n","    name='predicted impervious',\n","  ).add_to(map)\n","map.add_child(folium.LayerControl())\n","map"],"execution_count":0,"outputs":[]}]}