4 - Earth Engine and TensorFlow.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Earth Engine and TensorFlow in Cloud Datalab\n",
    "\n",
    "This notebook walks you through a simple example of using Earth Engine and TensorFlow together in Cloud Datalab.\n",
    "\n",
    "Specifically, we will train a neural network to recognize cloudy pixels in a Landsat scene. For this simple example we will use the output of the Fmask cloud detection algorithm as training data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Configure the Environment\n",
    "We begin by importing a number of useful libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import ee\n",
    "from IPython import display\n",
    "import math\n",
    "from matplotlib import pyplot\n",
    "import numpy\n",
    "from osgeo import gdal\n",
    "import tempfile\n",
    "import tensorflow as tf\n",
    "import urllib\n",
    "import zipfile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the Earth Engine client. This assumes that you have already configured Earth Engine credentials in this Datalab instance. If not, see the \"`Earth Engine Datalab Initialization.ipynb`\" notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspect the Input Data\n",
    "\n",
    "Load a Landsat image with corresponding Fmask label data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "input_image = ee.Image('LANDSAT/LT5_L1T_TOA_FMASK/LT50100551998003CPE00')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's define a helper function to make it easier to print thumbnails of Earth Engine images. (We'll be adding a library with utility functions like this one to the Earth Engine Python SDK, but for now we can do it by hand.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def print_image(image):\n",
    "  display.display(display.Image(ee.data.getThumbnail({\n",
    "      'image': image.serialize(),\n",
    "      'dimensions': '360',\n",
    "  })))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use our helper function to quickly visualize the image and label data. The Fmask values are:\n",
    "\n",
    "0 | 1 | 2 | 3 | 4\n",
    ":---:|:---:|:---:|:---:|:---:\n",
    "Clear | Water | Shadow | Snow | Cloud"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print_image(input_image.visualize(\n",
    "    bands=['B3', 'B2', 'B1'],\n",
    "    min=0,\n",
    "    max=0.3,\n",
    "))\n",
    "print_image(input_image.visualize(\n",
    "    bands=['fmask'],\n",
    "    min=0,\n",
    "    max=4,\n",
    "    palette=['808080', '0000C0', '404040', '00FFFF', 'FFFFFF'],\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fetch the Input Data\n",
    "\n",
    "First we define some helper functions to download raw data from Earth Engine as `numpy` arrays.\n",
    "\n",
    "We use the `getDownloadId()` function, which only works for modestly sized datasets. For larger datasets, a better approach would be to initiate a batch Export from Earth Engine to Cloud Storage, which you could easily manage right here in Datalab too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def download_tif(image, scale):\n",
    "  url = ee.data.makeDownloadUrl(ee.data.getDownloadId({\n",
    "        'image': image.serialize(),\n",
    "        'scale': '%d' % scale,\n",
    "        'filePerBand': 'false',\n",
    "        'name': 'data',\n",
    "  }))\n",
    "  local_zip, headers = urllib.urlretrieve(url)\n",
    "  with zipfile.ZipFile(local_zip) as local_zipfile:\n",
    "    return local_zipfile.extract('data.tif', tempfile.mkdtemp())\n",
    "\n",
    "def load_image(image, scale):\n",
    "  local_tif_filename = download_tif(image, scale)\n",
    "  dataset = gdal.Open(local_tif_filename, gdal.GA_ReadOnly)\n",
    "  bands = [dataset.GetRasterBand(i + 1).ReadAsArray() for i in range(dataset.RasterCount)]\n",
    "  return numpy.stack(bands, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use that function to load the data from Earth Engine, including a valid data band, as a `numpy` array. This may take a few seconds. We also convert the Fmask band  to a binary cloud label (i.e. `fmask`=4)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mask = input_image.mask().reduce('min')\n",
    "data = load_image(input_image.addBands(mask), scale=240)\n",
    "data[:,:,7] = numpy.equal(data[:,:,7], 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the local data. This time, for variety, we display it as an NRG false-color image. We can use `pyplot` to display local `numpy` arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pyplot.imshow(numpy.clip(data[:,:,[3,2,1]] * 3, 0, 1))\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocess the Input Data\n",
    "\n",
    "Select the valid pixels and hold out a fraction for use as validation data.  Compute per-band means and standard deviations of the training data for normalization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "HOLDOUT_FRACTION = 0.1\n",
    "\n",
    "# Reshape into a single vector of pixels.\n",
    "data_vector = data.reshape([data.shape[0] * data.shape[1], data.shape[2]])\n",
    "\n",
    "# Select only the valid data and shuffle it.\n",
    "valid_data = data_vector[numpy.equal(data_vector[:,8], 1)]\n",
    "numpy.random.shuffle(valid_data)\n",
    "\n",
    "# Hold out a fraction of the labeled data for validation.\n",
    "training_size = int(valid_data.shape[0] * (1 - HOLDOUT_FRACTION))\n",
    "training_data = valid_data[0:training_size,:]\n",
    "validation_data = valid_data[training_size:-1,:]\n",
    "\n",
    "# Compute per-band means and standard deviations of the input bands.\n",
    "data_mean = training_data[:,0:7].mean(0)\n",
    "data_std = training_data[:,0:7].std(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "valid_data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build the TensorFlow Model\n",
    "\n",
    "We start with a helper function to build a simple TensorFlow neural network layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_nn_layer(input, output_size):\n",
    "  input_size = input.get_shape().as_list()[1]\n",
    "  weights = tf.Variable(tf.truncated_normal(\n",
    "      [input_size, output_size],\n",
    "      stddev=1.0 / math.sqrt(float(input_size))))\n",
    "  biases = tf.Variable(tf.zeros([output_size]))\n",
    "  return tf.matmul(input, weights) + biases"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we define our TensorFlow model, a neural network with two hidden layers with tanh() nonlinearities. The main network has two outputs, continuous-valued “logits” representing non-cloud and cloud, respectively. The binary output is intepreted as the argmax of these outputs.\n",
    "\n",
    "We define a training step, which uses Kingma and Ba's Adam algorithm to minimize the cross-entropy between the logits and the training data. Finally, we define a simple overall percentage accuracy measure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "NUM_INPUT_BANDS = 7\n",
    "NUM_HIDDEN_1 = 20\n",
    "NUM_HIDDEN_2 = 20\n",
    "NUM_CLASSES = 2\n",
    "\n",
    "input = tf.placeholder(tf.float32, shape=[None, NUM_INPUT_BANDS])\n",
    "labels = tf.placeholder(tf.float32, shape=[None])\n",
    "\n",
    "normalized = (input - data_mean) / data_std\n",
    "hidden1 = tf.nn.tanh(make_nn_layer(normalized, NUM_HIDDEN_1))\n",
    "hidden2 = tf.nn.tanh(make_nn_layer(hidden1, NUM_HIDDEN_2))\n",
    "logits = make_nn_layer(hidden2, NUM_CLASSES)\n",
    "outputs = tf.argmax(logits, 1)\n",
    "\n",
    "int_labels = tf.to_int64(labels)\n",
    "cross_entropy = tf.nn.sparse_softmax_cross_entropy_with_logits(logits, int_labels, name='xentropy')\n",
    "train_step = tf.train.AdamOptimizer().minimize(cross_entropy)\n",
    "\n",
    "correct_prediction = tf.equal(outputs, int_labels)\n",
    "accuracy = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train the Neural Network\n",
    "\n",
    "Now train the neural network, using batches of training data drawn randomly from the training data pool. We periodically compute the accuracy against the validation data. When we're done training, we apply the model to the complete input data set.\n",
    "\n",
    "This simple notebook performs all TensorFlow operations locally. However, for larger analyses you could bring up a cluster of TensorFlow workers to parallelize the computation, all controlled from within Datalab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "BATCH_SIZE = 1000\n",
    "NUM_BATCHES = 1000\n",
    "\n",
    "with tf.Session() as sess:\n",
    "  sess.run(tf.initialize_all_variables())\n",
    "\n",
    "  validation_dict = {\n",
    "    input: validation_data[:,0:7],\n",
    "    labels: validation_data[:,7],\n",
    "  }\n",
    "\n",
    "  for i in range(NUM_BATCHES):\n",
    "    batch = training_data[numpy.random.choice(training_size, BATCH_SIZE, False),:]\n",
    "    train_step.run({input: batch[:,0:7], labels: batch[:,7]})\n",
    "\n",
    "    if i % 100 == 0 or i == NUM_BATCHES - 1:\n",
    "      print('Accuracy %.2f%% at step %d' % (accuracy.eval(validation_dict) * 100, i))\n",
    "\n",
    "  output_data = outputs.eval({input: data_vector[:,0:7]})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspect the Results\n",
    "\n",
    "Here we dislay the results. The red band corresponds to the TensorFlow output and the blue band corresponds to the labeled training data, so pixels that are red and blue correspond to disagreements between the model and the training data. (There aren't many: look carefully around the fringes of the clouds.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "output_image = output_data.reshape([data.shape[0], data.shape[1]])\n",
    "red = numpy.where(data[:,:,8], output_image, 0.5)\n",
    "blue = numpy.where(data[:,:,8], data[:,:,7], 0.5)\n",
    "green = numpy.minimum(red, blue)\n",
    "\n",
    "comparison_image = numpy.dstack((red, green, blue))\n",
    "pyplot.figure(figsize = (12,12))\n",
    "pyplot.imshow(comparison_image)\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can zoom in on a particular region over on the right side of the image to see some of the disagreements. Red pixels represent comission errors and blue pixels represent omission errors relative to the labeled input data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pyplot.figure(figsize = (12,12))\n",
    "pyplot.imshow(comparison_image[300:500,600:,:], interpolation='nearest')\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}