-
-
Save fdbesanto2/a05fdd13a7b6a953851810dece81dd8c to your computer and use it in GitHub Desktop.
s2cnnclass.ipynb
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| {"nbformat":4,"nbformat_minor":0,"metadata":{"accelerator":"GPU","colab":{"name":"s2cnnclass.ipynb","provenance":[],"collapsed_sections":[]},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.8"}},"cells":[{"cell_type":"markdown","metadata":{"id":"8QZMxH4E2Qj7"},"source":["# Sentinel-2 image classification with Convolutional Neural Networks\n","## Mort Canty\n","## ZFL Bonn, March, 2021"]},{"cell_type":"markdown","metadata":{"id":"5rJN97BoaEGK"},"source":["## Context\n","Image recognition with deep learning techniques such as the convolutional neural network (CNN) relies heavily on the availablity of vast amounts of labeled training data as well as a lot of computing resources. For typical remote sensing image classification these are usually not available. We'll explore an alternative approach which makes use of pre-trained networks, freely available Sentinel-2 imagery, the EuroSAT Sentinel-2 training database and TensorFlow.\n"," \n","### References\n","[Géron (2019) Hands-On Machine Learning ...](https://b-ok.cc/book/5341845/f49201)\n"," \n","[Helber et al. (2017) EuroSAT: A novel dataset and deep learning benchmark ...](https://www.researchgate.net/publication/319463676_EuroSAT_A_Novel_Dataset_and_Deep_Learning_Benchmark_for_Land_Use_and_Land_Cover_Classification)\n","\n","[Leitloff (2018) github.com/jensleitloff/CNN-Sentinel](https://github.com/jensleitloff/CNN-Sentinel)\n","\n","[Sumbul et al. (2020) BigEarthNet: Dataset with a new class-nomenclature for remote sensing image understanding](https://arxiv.org/abs/2001.06372)\n","\n","### Software\n","\n","[Python scripts](https://drive.google.com/drive/folders/1D2RAEu14Zcl9NzqWkSg4J7MilHGHi0BH?usp=sharing)"]},{"cell_type":"markdown","metadata":{"id":"UMYv3ICRdMwS"},"source":["## Preliminaries"]},{"cell_type":"code","metadata":{"id":"EkILT8AAwiOq"},"source":["import numpy as np\n","import tensorflow as tf\n","import tensorflow.keras as keras\n","import tensorflow_datasets as tfds\n","import matplotlib.pyplot as plt\n","from google.colab import drive, files\n","import IPython.display as disp\n"," \n","%matplotlib inline"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"85Xux7Nu8QWi"},"source":["import ee\n","ee.Authenticate()\n","\n","ee.Initialize()"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"F7UKmTCe-Pjf"},"source":["drive.mount('/content/drive')"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"QKc0QFSc-clq"},"source":["!ls /content/drive/MyDrive/scripts"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"hpyBGnXD3FMf"},"source":["## 1. The fashion/mnist dataset\n","\n","We begin with a simple \"hello world\" type example to get familiar with convolutional networks:\n","\n","MNIST: Modified National Institute of Standards and Technology dataset\n","\n","https://medium.com/tensorflow/hello-deep-learning-fashion-mnist-with-keras-50fcff8cd74a"]},{"cell_type":"markdown","metadata":{"id":"PxZH3Zix5E5y"},"source":["Load the fashion-mnist dataset with the keras datasets model.\n","\n","\"*The tf.keras.datasets module provides a few toy datasets (already-vectorized, in Numpy format) that can be used for debugging a model or creating simple code examples.*\""]},{"cell_type":"code","metadata":{"id":"C0MEo3MJ5PK5"},"source":["(x_train, y_train), (x_test, y_test) = keras.datasets.fashion_mnist.load_data()\n","print(\"x_train shape:\", x_train.shape, \"y_train shape:\", y_train.shape)\n","\n","# normalize the train/test data\n","x_train = x_train.astype('float32') / 255\n","x_test = x_test.astype('float32') / 255\n","\n","plt.imshow(x_test[55],cmap='gray',vmin=0,vmax=1)\n","plt.show()"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"ObOf-4iz5p54"},"source":["# Further break training data into train / validation sets \n","# put 5000 into validation set and keep remaining 55,000 for train\n","(x_train, x_valid) = x_train[5000:], x_train[:5000] \n","(y_train, y_valid) = y_train[5000:], y_train[:5000]\n","\n","# Reshape input data from (28, 28) to (28, 28, 1)\n","w, h = 28, 28\n","x_train = x_train.reshape(x_train.shape[0], w, h, 1)\n","x_valid = x_valid.reshape(x_valid.shape[0], w, h, 1)\n","x_test = x_test.reshape(x_test.shape[0], w, h, 1)\n","\n","# One-hot encode the labels\n","y_train = tf.keras.utils.to_categorical(y_train, 10)\n","y_valid = tf.keras.utils.to_categorical(y_valid, 10)\n","y_test = tf.keras.utils.to_categorical(y_test, 10)\n","\n","# Print training set shape\n","print(\"x_train shape:\", x_train.shape, \"y_train shape:\", y_train.shape)\n","\n","# Print the number of training, validation, and test datasets\n","print(x_train.shape[0], 'train set')\n","print(x_valid.shape[0], 'validation set')\n","print(x_test.shape[0], 'test set')"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"AIme3k1v9FjP"},"source":["#### Build the model \n","https://towardsdatascience.com/mnist-handwritten-digits-classification-using-a-convolutional-neural-network-cnn-af5fafbc35e9\n","\n",""]},{"cell_type":"markdown","metadata":{"id":"kAVcPA6Silpv"},"source":["### TensorFlow high level sequential API"]},{"cell_type":"code","metadata":{"id":"v52uJPcb-kmz"},"source":["# Use the sequential API\n","model = tf.keras.Sequential()\n","# Must define the input shape in the first layer of the neural network\n","model.add(tf.keras.layers.Conv2D(filters=32, kernel_size=3, padding='same', activation='relu', \n"," input_shape=(28,28,1))) \n","model.add(tf.keras.layers.MaxPooling2D(pool_size=2))\n","model.add(tf.keras.layers.Conv2D(filters=64, kernel_size=3, padding='same', activation='relu'))\n","model.add(tf.keras.layers.MaxPooling2D(pool_size=2))\n","model.add(tf.keras.layers.Flatten())\n","model.add(tf.keras.layers.Dense(128, activation='relu'))\n","model.add(tf.keras.layers.Dense(10, activation='softmax'))\n","# Take a look at the model summary\n","model.summary()"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"O259IRqgfj8A"},"source":["keras.utils.plot_model(model, \"mnist_model.png\",show_shapes=True)"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"FwztPZhZ-vWC"},"source":["# conv2d: (3 times 3 pixel inputs + 1 bias) times 32 filters\n","conv2d1 = (3*3 + 1)*32\n","print('conv2d1 %i'%conv2d1)\n","\n","# conv2d: (3 times 3 pixel inputs times 32 filters + 1 bias) times 64 filters\n","conv2d2 = (3*3*32 + 1)*64\n","print('conv2d2 %i'%conv2d2)\n","\n","# flatten: 7*7*64\n","flatten = 7*7*64\n","print('flatten %i'%flatten)\n","\n","# dense (flatten inputs + 1 bias) times 128 neurons \n","dense1 = (flatten + 1)*128\n","print('dense1 %i'%dense1)\n","\n","# dense (128 inputs + 1 bias) times 10\n","dense2 = (128+1)*10\n","print('dense2 %i'%dense2)\n","\n","print('trainable parameters %i'%(conv2d1+conv2d2+dense1+dense2))"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"ju8mJS_j-0qy"},"source":["# summarize feature map shapes\n","for i in range(len(model.layers)):\n"," layer = model.layers[i]\n"," # check for convolutional layer\n"," if 'conv' not in layer.name:\n"," continue\n"," # summarize output shape\n"," print(i, layer.name, layer.output.shape)\n"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"cHshcTRJSAzS"},"source":[""]},{"cell_type":"markdown","metadata":{"id":"OWP_tW84-7wS"},"source":["#### Compile the model"]},{"cell_type":"code","metadata":{"id":"BnMPjQcj-6PE"},"source":["model.compile(loss='categorical_crossentropy',\n"," optimizer='adam',\n"," metrics=['accuracy'])"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"toBRQkXW_U_V"},"source":["#### Train the model"]},{"cell_type":"code","metadata":{"id":"l6BxEk2z_TY-"},"source":["%%time\n","from tensorflow.keras.callbacks import ModelCheckpoint,TensorBoard\n","from datetime import datetime\n","now = datetime.utcnow().strftime(\"%Y%m%d%H%M%S\")\n","logdir = 'tf_cnn/run-'+str(now)\n","# default: save checkpoint after each epoch\n","checkpointer = ModelCheckpoint(filepath='model.weights.best.hdf5', verbose = 1, save_best_only=True)\n","tensorboardlog = TensorBoard(logdir)\n","\n","model.fit(x_train,\n"," y_train,\n"," batch_size=64,\n"," epochs=10,\n"," validation_data=(x_valid, y_valid),\n"," callbacks=[checkpointer,tensorboardlog])"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"OE6GsKVNA4J1"},"source":["%load_ext tensorboard\n","%tensorboard --logdir tf_cnn"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"Dxmz1zIjFMEX"},"source":["#### Evaluate the model"]},{"cell_type":"code","metadata":{"id":"5V5jF21cFDow"},"source":["# Load the weights with the best validation accuracy\n","model.load_weights('model.weights.best.hdf5')\n","\n","# Evaluate the model on test set\n","score = model.evaluate(x_test, y_test, verbose=0)\n","\n","# Print test accuracy\n","print( '\\nTest accuracy:', score )"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"Ygp2mCkqFedJ"},"source":["#### Dig into the model\n","https://machinelearningmastery.com/how-to-visualize-filters-and-feature-maps-in-convolutional-neural-networks/"]},{"cell_type":"code","metadata":{"id":"6DU64AejFr75"},"source":["# redefine model to output right after the first pooling layer\n","model1 = tf.keras.Model(inputs=model.inputs, outputs=model.layers[1].output)\n","\n","x = x_train[0].reshape(1,28,28,1)\n","\n","featuremap = model1.predict(x)\n","featuremap.shape"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"eldDW_vxF75D"},"source":["# plot all 32 maps \n","ix = 1\n","for _ in range(4):\n"," for _ in range(8):\n"," # specify subplot and turn of axis\n"," ax = plt.subplot(4, 8, ix)\n"," ax.set_xticks([])\n"," ax.set_yticks([])\n"," # plot filter channel in grayscale\n"," plt.imshow(featuremap[0, :, :, ix-1], cmap='gray')\n"," ix += 1\n","# show the figure\n","plt.show()"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"bjOiJpWfGELA"},"source":["# redefine model to output right after the second pooling layer\n","model2 = tf.keras.Model(inputs=model.inputs, outputs=model.layers[3].output)\n","\n","featuremap = model2.predict(x)\n","featuremap.shape"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"GfNKsmy_GNpT"},"source":["# plot first 32 maps \n","ix = 1\n","for _ in range(4):\n"," for _ in range(8):\n"," # specify subplot and turn of axis\n"," ax = plt.subplot(4, 8, ix)\n"," ax.set_xticks([])\n"," ax.set_yticks([])\n"," # plot filter channel in grayscale\n"," plt.imshow(featuremap[0, :, :, ix-1], cmap='gray')\n"," ix += 1\n","# show the figure\n","plt.show()"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"H5FA1J4hwiOo"},"source":["## 2. Transfer Learning with the EuroSAT/rgb dataset\n","Transfer Learning is a deep learning technique whereby a model developed for a task is reused as the initial point for a model in another domain. Rather than training the model from scratch, the pre-trained model is modified, partially retrained and re-used. "]},{"cell_type":"markdown","metadata":{"id":"yypIvf5SsAa8"},"source":["### EUROSAT dataset (27000 Sentinel-2 RGB 64x64 images with labels)\n","\n","This dataset is part of the TFDS project. We use it to download the EuroSAT database and split it into test, train and validate subsets."]},{"cell_type":"code","metadata":{"id":"cppLiZH7wiOr"},"source":["(test_set, valid_set, train_set), info = tfds.load('eurosat/rgb', split = ['train[:10%]', 'train[10%:25]', 'train[25%:]' ], as_supervised=True, with_info=True)"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"pLS8k01QBVNr"},"source":["info.description"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"NMon_RHQBAEP"},"source":["info.features"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"buY_uc4EB7bm"},"source":["info.features['label'].names"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"LYl3MQinPLGD"},"source":["train_set"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"z5ZhVA3lTTvR"},"source":["fig=tfds.show_examples(train_set, info, rows=4, cols=4)"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"ooeWKs7GtTkD"},"source":["### Preprocess for Xception cnn\n","We will use the [Xception CNN](https://keras.io/api/applications/xception/) which has been pre-trained with the [IMAGENET](http://image-net.org/) database. The CNN requires image patches with dimension 224x224 and has its own preprocessor:"]},{"cell_type":"code","metadata":{"id":"Jh7SojXbwiOw"},"source":["def preprocess(image,label):\n"," resized_image = tf.image.resize(image, [224,224])\n"," final_image = keras.applications.xception.preprocess_input(resized_image)\n"," return final_image, label\n","\n","batch_size = 32\n","# Randomly shuffle the elements of the training dataset.\n","train_set = train_set.shuffle(1000)\n","# map the preprocessing function over the three datasets\n","train_set = train_set.map(preprocess).batch(batch_size).prefetch(1)\n","valid_set = valid_set.map(preprocess).batch(batch_size).prefetch(1)\n","test_set = test_set.map(preprocess).batch(batch_size).prefetch(1)\n","\n","test_set"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"YL-KI3eyuAWO"},"source":["For a GPU runtime, the _prefetch_ method makes the CPU prepare the next batch while the GPU is processing the current one.\n","\n","### Setup the pre-trained xception model using weights obtained for the IMAGENET dataset\n","\n","So far we have built our networks with the Keras sequential API. Here we use the more flexible _functional API_."]},{"cell_type":"code","metadata":{"id":"QJ-ck6qswiOx"},"source":["# download xception (base) model without top layers\n","base_model = keras.applications.xception.Xception(weights='imagenet', include_top=False, input_shape=(224,224,3))\n","# add a pooling layer with input from the base model (averages over all output feature maps of the base model)\n","avg = keras.layers.GlobalAveragePooling2D()(base_model.output)\n","# add an output softmax layer with 10 neurons and softmax activation\n","output = keras.layers.Dense(10, activation='softmax')(avg)\n","# define the model\n","model = keras.Model(inputs=base_model.input, outputs=output)"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"RoKRCY49gDgO"},"source":["keras.utils.plot_model(model, \"eurosat_model.png\", dpi=48, show_shapes=True)"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"BLbP2Ww_vd98"},"source":["Freeze the layers in the base model"]},{"cell_type":"code","metadata":{"id":"XZknAi1swiOx"},"source":["for layer in base_model.layers:\n"," layer.trainable = False"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"bHX0hrVPvv7_"},"source":["Define optimization method and compile the model"]},{"cell_type":"code","metadata":{"id":"AA0ncLYuwiOx"},"source":["# scaled gradient descent optimizer\n","optimizer = keras.optimizers.SGD(lr=0.2, momentum=0.9, decay=0.01)\n","# cross entropy loss (for no one-hot labels encoding)\n","loss = \"sparse_categorical_crossentropy\"\n","# only display accuracy = overall correct / total examples\n","metrics = ['accuracy']\n","\n","model.compile(loss=loss, optimizer=optimizer, metrics=metrics)"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"10XeLpSxv6uZ"},"source":["Train the model, top layer weights only (use GPU runtime!)"]},{"cell_type":"code","metadata":{"id":"ww-rh-SBwiOy","scrolled":false},"source":["history = model.fit(train_set, epochs=10, validation_data=valid_set)"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"PrqCGhcswn5l"},"source":["Fine tune the rest of the weights"]},{"cell_type":"code","metadata":{"id":"HAqxpYWKwiOy"},"source":["for layer in base_model.layers:\n"," layer.trainable = True\n","optimizer = keras.optimizers.SGD(lr=0.01, momentum=0.9, decay=0.001)\n","model.compile(loss=\"sparse_categorical_crossentropy\", optimizer=optimizer, metrics=['accuracy']) \n","history = model.fit(train_set, epochs=5, validation_data=valid_set)"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"uaLudiLZymyQ"},"source":["Test (on first 32 test examples)"]},{"cell_type":"code","metadata":{"id":"rKzuN9E-yhDl"},"source":["print( np.argmax(model.predict(test_set.take(1)),1) )\n","for image, label in test_set.take(1): \n"," print(np.array(label))"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"EDIkLl2iPXvj"},"source":["model.evaluate(test_set)"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"_Hn5O2XyECqb"},"source":[" Save the model"]},{"cell_type":"code","metadata":{"id":"tsr6P_L-gRDk"},"source":["model.save('/content/drive/MyDrive/gee/eurosat_model.h5')"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"lQxiRQ3Wf5OC"},"source":["from google.colab import files\n","files.download('/content/drive/MyDrive/gee/eurosat_model.h5')"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"9GB1m_LKHIap"},"source":["Export a Sentinel-2 image to Google Drive"]},{"cell_type":"code","metadata":{"id":"-qQjm4MVH1C2"},"source":["geoJSON = {\n"," \"type\": \"FeatureCollection\",\n"," \"features\": [\n"," {\n"," \"type\": \"Feature\",\n"," \"properties\": {},\n"," \"geometry\": {\n"," \"type\": \"Polygon\",\n"," \"coordinates\": [\n"," [\n"," [\n"," 5.9710693359375,\n"," 50.47848271564207\n"," ],\n"," [\n"," 7.2125244140625,\n"," 50.47848271564207\n"," ],\n"," [\n"," 7.2125244140625,\n"," 51.156954270151346\n"," ],\n"," [\n"," 5.9710693359375,\n"," 51.156954270151346\n"," ],\n"," [\n"," 5.9710693359375,\n"," 50.47848271564207\n"," ]\n"," ]\n"," ]\n"," }\n"," }\n"," ]\n","}\n","\n","coords = geoJSON['features'][0]['geometry']['coordinates']\n","aoi = ee.Geometry.Polygon(coords)\n","\n","s2 = ee.ImageCollection('COPERNICUS/S2') \\\n"," .filterBounds(aoi) \\\n"," .filterDate(ee.Date('2018-05-07'),ee.Date('2018-05-09')) \\\n"," .filter(ee.Filter.contains(rightValue=aoi,leftField='.geo')) \\\n"," .first() \\\n"," .select('B2','B3','B4')\n","\n","gdexport = ee.batch.Export.image.toDrive(s2.clip(aoi),\n"," description='driveExportTask_s2', \n"," folder = 'gee',\n"," fileNamePrefix = 's2_20180508',\n"," scale = 10,\n"," maxPixels = 1e11)\n","gdexport.start()"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"C4ePIXf-sXoK"},"source":["run /content/drive/MyDrive/scripts/dispms -f /content/drive/MyDrive/gee/s2_20180508.tif -e 4 -p [1,2,3] -d [500,500,4000,4000]"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"kA4KNE7qrK8l"},"source":["The classification is carried out by a Python script S2cnnclassify.py using the saved model:"]},{"cell_type":"code","metadata":{"id":"BJBU-jhQZgbj"},"source":["%run /content/drive/MyDrive/scripts/S2cnnclassify -d [500,500,4000,4000] -s 8 /content/drive/MyDrive/gee/s2_20180508.tif "],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"bnh_pnlunawU"},"source":["In order to visualize the result it is convenient to view it in the GEE code editor (download from Drive and upload to assets).\n","\n","https://code.earthengine.google.com/b00ca86e09a78b0970ef41b9582b7665"]},{"cell_type":"markdown","metadata":{"id":"FcZJP8psLXXW"},"source":["### Project\n","The 10 target classes should be seasonally invariant. Try to determine if this is the case."]},{"cell_type":"markdown","metadata":{"id":"hfxWYOFCvUM2"},"source":["### Outlook\n","In the third and final part of the course we will return to analysis of time series of SAR imagery and examine a sequential change detection method. The Colab notebooks are on the [Earthengine community tutorial website.](https://developers.google.com/earth-engine/tutorials/community/detecting-changes-in-sentinel-1-imagery-pt-1)"]}]} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment