{
"cells": [
{
"cell_type": "markdown",
"id": "4f26dc3c",
"metadata": {},
"source": [
"Lecture 09: Machine Learning 3/3 (exercises)"
]
},
{
"cell_type": "markdown",
"id": "293ce49e",
"metadata": {},
"source": [
"# EX: Semantic classification of satellite image (part 2/2)"
]
},
{
"cell_type": "markdown",
"id": "ffa45a51",
"metadata": {},
"source": [
"\n",
"This exercice has been developped by Andreas Ley and Ronny Hänsch (TU-Berlin) in the framework of the GEO.X 2017 Autumn school.
\n",
"
\n",
"The exercice is designed to implement the classification of a satellite image (Sentinel-2) into semantic labels defining land use: forest, fields, urban, water.
\n",
"
\n",
"The exercice is split into 2 parts: this first part (lecture 08) will implement a PCA analisis on the image crops in order to reduce dimensionality. The second part (lecture 09) will use the outputs of the PCA and implement a SVM to classify the pixels into different land use classes.\n",
" \n",
"
"
]
},
{
"cell_type": "markdown",
"id": "4be3257c",
"metadata": {},
"source": [
"\n",
"
Description of the dataset:
\n",
"
\n",
"The data stored in the file \"images/s2_training_data.npz\" contains image crops of a Sentinel-2 multispectral image, and some semantic labels of land use.
\n",
"
\n",
" The file can be loaded using
numpy.load(), and contains 2 parts:
\n",
"
\n",
" - data: contains 20,000 patches (crops), each of 15x15 pixels and 4 channels (R + G + B + shortwave-infrared). The data is organized as one numpy array of shape (20000, 4, 15, 15), where the first dimension are the instances, the next dimension are the channels, and the latter two are width and height.
The first 5000 crops contain forest in the central pixel, the next 5000 contain fields/lower vegetation, the next 5000 are urban areas, and the last 5000 are water. \n",
" - labels: contains semantic labels of land use (forest, fields, urban, water), which will be later used for training the classifier.
\n",
"
\n",
"
\n",
"In this second part of the exercice, in addition to the training data you will be using two other files:
\n",
"
\n",
" - s2_testing_data.npz: the file is structured like the training file \"s2_training_data.npz\", and will be used to evaluate our SVM classifier.
\n",
"
\n",
" - s2_application_data.npz: the file only contains a \"data\" field, which stores a single Sentinel-2 image of 512x512 pixels, with 4 channels (R+G+B+SWIR). This image will be used to apply the trained SVM classifier in order to classify the pixels into the different land use classes (forest, fields, urban, water).
\n",
"
\n",
"
\n",
" The datasets contain 5k crops per class. Each crop is 15x15 pixels in size with 4 channels per pixel. The label of the central pixel is given in a separate array. Since the input dimensionality of each crop is quite high (4x15x15 = 900), we want to compress it using PCA.
\n",
" The following function have been developped in the 1st part of the exercice.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "b8891fd1",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 45,
"id": "7b8b8773",
"metadata": {},
"outputs": [],
"source": [
"def show_raw_image(img):\n",
"\n",
" img2 = np.log(img[[2,1,0],:,:])\n",
" \n",
" img2[0,:,:] = img2[0,:,:].copy() * 1.05303 + -6.32792\n",
" img2[1,:,:] = img2[1,:,:].copy() * 1.74001 + -10.8407\n",
" img2[2,:,:] = img2[2,:,:].copy() * 1.20697 + -6.73016\n",
"\n",
" img2 = np.clip(img2 / 6 + 0.5, 0.0, 1.0)\n",
"\n",
" plt.imshow(np.transpose(img2, (1, 2, 0)))\n",
" #plt.show()\n",
"\n",
"\n",
"def compute_mean_PCs(X):\n",
" mean = np.mean(X, 0) # Second parameter is the axis along which we want to average (0 == across instances)\n",
" mean_free = X-mean\n",
"\n",
" vars_per_img = X.shape[1] * X.shape[2] * X.shape[3]\n",
" num_imgs = mean_free.shape[0]\n",
"\n",
" # Flatten data in to vectors\n",
" mean_free_vectorized = np.reshape(mean_free, (num_imgs, vars_per_img))\n",
"\n",
" # Increase this to speed up debugging\n",
" covar_subsampling = 2\n",
"\n",
" # Accumulate covar matrix\n",
" covar = np.zeros((vars_per_img, vars_per_img))\n",
" print(\"Image: 0\")\n",
" for i in range(0, num_imgs, covar_subsampling):\n",
" print(\"\\rImage: {}\".format(i))\n",
" covar += np.outer(mean_free_vectorized[i,:], mean_free_vectorized[i,:])\n",
"\n",
" covar /= num_imgs/covar_subsampling\n",
"\n",
" eig_val, eig_vec = np.linalg.eig(covar)\n",
"\n",
" # Sort by importance\n",
" idx = np.argsort(eig_val)[::-1]\n",
" eig_vec = eig_vec[:,idx]\n",
" eig_val = eig_val[idx]\n",
"\n",
" # Reshape data back into images. Note that eig_vec is the transpose of what you might expect it to be.\n",
" principal_components = np.transpose(eig_vec, (1,0)).reshape((vars_per_img, X.shape[1], X.shape[2], X.shape[3]))\n",
"\n",
" return mean, principal_components, eig_val\n",
"\n",
"\n",
"\n",
"def show_first_principal_components(pcs):\n",
"\n",
" f, axarr = plt.subplots(8,8)\n",
" for i in range(0,8):\n",
" for j in range(0,8):\n",
" img2 = pcs[i*8+j,[2,1,0],:,:]\n",
" img2 = np.clip(img2 * 10 + 0.5, 0.0, 1.0)\n",
" axarr[i,j].imshow(np.transpose(img2, (1, 2, 0)))\n",
"\n",
" plt.show()\n",
"\n",
"def compute_features(X, mean, principal_components, count):\n",
" X_mean_free = X - mean\n",
" features = np.zeros((X.shape[0], count))\n",
" for i in range(0, X.shape[0]):\n",
" for j in range(0, count):\n",
" # Note: The [i,:,:,:] is being very explicit here. [i] would also work.\n",
" features[i,j] = X_mean_free[i,:,:,:].flatten().dot(principal_components[j,:,:,:].flatten())\n",
" return features\n",
" \n",
"def reconstruct_image(feature, mean, principal_components):\n",
" reconstruction = np.copy(mean)\n",
" for i in range(0, feature.shape[0]):\n",
" reconstruction += feature[i] * principal_components[i,:,:,:]\n",
" return reconstruction "
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "97b5cbee",
"metadata": {},
"outputs": [],
"source": [
"def testAccuracy(classifier, test_X, test_Y):\n",
" inferred_Y = classifier.predict(test_X)\n",
" return np.mean(test_Y == inferred_Y)"
]
},
{
"cell_type": "markdown",
"id": "0f07f28b",
"metadata": {},
"source": [
"## run PCA on training_data"
]
},
{
"cell_type": "markdown",
"id": "23deb622",
"metadata": {},
"source": [
"\n",
" Use the functions implemented for the PCA exercise to compute the mean and principal components on the training data (s2_training_data.npz).
\n",
" Compute a 16-dimensional feature vector for each training sample by using the \"compute_features\" function.\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "e1915902",
"metadata": {},
"source": [
"## train a SVM on PCA features"
]
},
{
"cell_type": "markdown",
"id": "d3fa7895",
"metadata": {},
"source": [
"\n",
" Train a support vector machine on the computed features and the corresponding labels. Use
sklearn.svm.SVC().\n",
"
\n",
" Compute the accuracy on the training dataset using the function testAccuracy() defined above.\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "3eb63d67",
"metadata": {},
"source": [
"## load test dataset and evaluate classifier"
]
},
{
"cell_type": "markdown",
"id": "3451b503",
"metadata": {},
"source": [
"\n",
" 1. Load the test dataset from the file \"s2_testing_data.npz\".
\n",
" 2. Compute features for it as well (remember to compute the logarithm of the data first) and evaluate the classifier on the test dataset by testing it's accuracy.
\n",
" 3. Optional: try out the training accuracy for 16, 32, and 256 principal components. What is the test accuracy for those? How long does training/inference take (roughly)?\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "51f52097",
"metadata": {},
"source": [
"## classify pixels in application dataset"
]
},
{
"cell_type": "markdown",
"id": "1004d454",
"metadata": {},
"source": [
"\n",
" 1. Load the image \"data/s2_application_data.npz\" and plot using show_raw_image().
\n",
" Have a look at it: can you spot the forest, urban, field, and water regions?
\n",
"
\n",
" 2. Compute features on the training dataset with 64 components, and train SVM classifier with these features.
\n",
"
\n",
" 3. Classify the image pixels into the four semantic classes (Forest, Field, Urban, Water):
\n",
"
\n",
" Use a \"sliding window\" approach where you crop a window of 15x15 pixels from the source image, that you will use to infer the label of the central pixel. To do so, compute the crop features, and use these to \"predict\" the pixel class with your SVM classifier. Store the predicted pixel value in an array \"application_labels\" (shape = (application_X.shape[1]-15, application_X.shape[2]-15)). Then shift the crop by one pixel and repeat these steps in order to classify the next pixel. Continue sliding the crop across the entire image in order to classify the entire image.
\n",
"Note: compute_features expects an array of crops instead of a single crop which means that you need to turn your crop into an array containing just one crop.
\n",
"
\n",
" 4. Render the final label map with: plt.imshow(application_labels).
\n",
" Is the classifier able to correctly label the different regions (forest, urban, field, water) you spotted in the image?
\n",
" How could you improve it?\n",
"
"
]
}
],
"metadata": {
"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.8.8"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {
"height": "calc(100% - 180px)",
"left": "10px",
"top": "150px",
"width": "282.883px"
},
"toc_section_display": true,
"toc_window_display": true
}
},
"nbformat": 4,
"nbformat_minor": 5
}