{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"UNET_regression_demo.ipynb","provenance":[{"file_id":"https://github.com/google/earthengine-api/blob/master/python/examples/ipynb/UNET_regression_demo.ipynb","timestamp":1586992475463}],"private_outputs":true,"collapsed_sections":[],"toc_visible":true,"machine_shape":"hm"},"kernelspec":{"name":"python3","display_name":"Python 3"},"accelerator":"GPU"},"cells":[{"cell_type":"code","metadata":{"id":"esIMGVxhDI0f","colab_type":"code","colab":{}},"source":["#@title Copyright 2020 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":"aV1xZ1CPi3Nw","colab_type":"text"},"source":["
\n","\n"," Run in Google Colab\n","\n"," View source on GitHub
"]},{"cell_type":"markdown","metadata":{"id":"_SHAc5qbiR8l","colab_type":"text"},"source":["# Introduction\n","\n","This is an Earth Engine <> TensorFlow demonstration notebook. Suppose you want to predict a continuous output (regression) from a stack of continuous inputs. In this example, the output is impervious surface area from [NLCD](https://www.mrlc.gov/data) and the input is a Landsat 8 composite. The model is a [fully convolutional neural network (FCNN)](https://www.cv-foundation.org/openaccess/content_cvpr_2015/papers/Long_Fully_Convolutional_Networks_2015_CVPR_paper.pdf), specifically [U-net](https://arxiv.org/abs/1505.04597). This notebook shows:\n","\n","1. Exporting training/testing patches from Earth Engine, suitable for training an FCNN model.\n","2. Preprocessing.\n","3. Training and validating an FCNN model.\n","4. Making predictions with the trained model and importing them to Earth Engine."]},{"cell_type":"markdown","metadata":{"id":"_MJ4kW1pEhwP","colab_type":"text"},"source":["# Setup software libraries\n","\n","Authenticate and import 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":"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","print(tf.__version__)"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"n1hFdpBQfyhN","colab_type":"code","colab":{}},"source":["# Folium setup.\n","import folium\n","print(folium.__version__)"],"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":["# INSERT YOUR BUCKET HERE:\n","BUCKET = 'your-bucket-name'"],"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 = 'fcnn-demo'\n","TRAINING_BASE = 'training_patches'\n","EVAL_BASE = 'eval_patches'\n","\n","# Specify inputs (Landsat bands) to the model and the response variable.\n","opticalBands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7']\n","thermalBands = ['B10', 'B11']\n","BANDS = opticalBands + thermalBands\n","RESPONSE = 'impervious'\n","FEATURES = BANDS + [RESPONSE]\n","\n","# Specify the size and shape of patches expected by the model.\n","KERNEL_SIZE = 256\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 = 16000\n","EVAL_SIZE = 8000\n","\n","# Specify model training parameters.\n","BATCH_SIZE = 16\n","EPOCHS = 10\n","BUFFER_SIZE = 2000\n","OPTIMIZER = 'SGD'\n","LOSS = 'MeanSquaredError'\n","METRICS = ['RootMeanSquaredError']"],"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","l8sr = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')\n","\n","# Cloud masking function.\n","def maskL8sr(image):\n"," cloudShadowBitMask = ee.Number(2).pow(3).int()\n"," cloudsBitMask = ee.Number(2).pow(5).int()\n"," qa = image.select('pixel_qa')\n"," mask1 = qa.bitwiseAnd(cloudShadowBitMask).eq(0).And(\n"," qa.bitwiseAnd(cloudsBitMask).eq(0))\n"," mask2 = image.mask().reduce('min')\n"," mask3 = image.select(opticalBands).gt(0).And(\n"," image.select(opticalBands).lt(10000)).reduce('min')\n"," mask = mask1.And(mask2).And(mask3)\n"," return image.select(opticalBands).divide(10000).addBands(\n"," image.select(thermalBands).divide(10).clamp(273.15, 373.15)\n"," .subtract(273.15).divide(100)).updateMask(mask)\n","\n","# The image input data is a cloud-masked median composite.\n","image = l8sr.filterDate('2015-01-01', '2017-12-31').map(maskL8sr).median()\n","\n","# Use folium to visualize the imagery.\n","mapid = image.getMapId({'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.3})\n","map = folium.Map(location=[38., -122.5])\n","folium.TileLayer(\n"," tiles=mapid['tile_fetcher'].url_format,\n"," attr='Map Data © Google Earth Engine',\n"," overlay=True,\n"," name='median composite',\n"," ).add_to(map)\n","\n","mapid = image.getMapId({'bands': ['B10'], 'min': 0, 'max': 0.5})\n","folium.TileLayer(\n"," tiles=mapid['tile_fetcher'].url_format,\n"," attr='Map Data © Google Earth Engine',\n"," overlay=True,\n"," name='thermal',\n"," ).add_to(map)\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":"e0wHDyxVirec","colab_type":"code","colab":{}},"source":["nlcd = ee.Image('USGS/NLCD/NLCD2016').select('impervious')\n","nlcd = nlcd.divide(100).float()\n","\n","mapid = nlcd.getMapId({'min': 0, 'max': 1})\n","map = folium.Map(location=[38., -122.5])\n","folium.TileLayer(\n"," tiles=mapid['tile_fetcher'].url_format,\n"," attr='Map Data © Google Earth Engine',\n"," overlay=True,\n"," name='nlcd impervious',\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"," nlcd.select(RESPONSE)\n","]).float()\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":["trainingPolys = ee.FeatureCollection('projects/google/DemoTrainingGeometries')\n","evalPolys = ee.FeatureCollection('projects/google/DemoEvalGeometries')\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=[38., -100.], zoom_start=5)\n","folium.TileLayer(\n"," tiles=mapid['tile_fetcher'].url_format,\n"," attr='Map Data © 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":"both","colab":{}},"source":["# Convert the feature collections to lists for iteration.\n","trainingPolysList = trainingPolys.toList(trainingPolys.size())\n","evalPolysList = evalPolys.toList(evalPolys.size())\n","\n","# These numbers determined experimentally.\n","n = 200 # Number of shards in each polygon.\n","N = 2000 # Total sample size in each polygon.\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":"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 tuple 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"," 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.io.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","\tdataset = get_dataset(glob)\n","\tdataset = dataset.shuffle(BUFFER_SIZE).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":"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","\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":"markdown","metadata":{"id":"9JIE7Yl87lgU","colab_type":"text"},"source":["# Model\n","\n","Here we use the Keras implementation of the U-Net model. 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":"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=optimizers.get(OPTIMIZER), \n","\t\tloss=losses.get(LOSS),\n","\t\tmetrics=[metrics.get(metric) for metric in METRICS])\n","\n","\treturn model"],"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","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)"],"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":"code","metadata":{"id":"-RJpNfEUS1qp","colab_type":"code","colab":{}},"source":["# Load a trained model. 50 epochs. 25 hours. Final RMSE ~0.08.\n","MODEL_DIR = 'gs://ee-docs-demos/fcnn-demo/trainer/model'\n","m = tf.keras.models.load_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":"M3WDAa-RUpXP","colab_type":"code","colab":{}},"source":["def doExport(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),\n"," description = out_image_base,\n"," bucket = BUCKET,\n"," fileNamePrefix = FOLDER + '/' + out_image_base,\n"," region = region.getInfo()['coordinates'],\n"," scale = 30,\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}\n","\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.io.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.io.parse_single_example(example_proto, imageFeaturesDict)\n","\n"," def toTupleImage(inputs):\n"," inputsList = [inputs.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 + '/' + out_image_base + '.TFRecord'\n"," writer = tf.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"," 'impervious': 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":["# Output assets folder: YOUR FOLDER\n","user_folder = 'users/username' # INSERT YOUR FOLDER HERE.\n","\n","# Base file name to use for TFRecord files and assets.\n","bj_image_base = 'FCNN_demo_beijing_384_'\n","# Half this will extend on the sides of each patch.\n","bj_kernel_buffer = [128, 128]\n","# Beijing\n","bj_region = ee.Geometry.Polygon(\n"," [[[115.9662455210937, 40.121362012835235],\n"," [115.9662455210937, 39.64293313749715],\n"," [117.01818643906245, 39.64293313749715],\n"," [117.01818643906245, 40.121362012835235]]], None, False)"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"lLNEOLkXWvSi","colab_type":"code","colab":{}},"source":["# Run the export.\n","doExport(bj_image_base, bj_kernel_buffer, bj_region)"],"execution_count":0,"outputs":[]},{"cell_type":"code","metadata":{"id":"KxACnxKFrQ_J","colab_type":"code","colab":{}},"source":["# Run the prediction.\n","doPrediction(bj_image_base, user_folder, bj_kernel_buffer, bj_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=mapid['tile_fetcher'].url_format,\n"," attr='Map Data © 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":[]}]}