\n",
"This exercise has been developed by Andreas Ley and Ronny Hänsch (TU-Berlin) in the framework of the GEO.X 2017 Autumn school. \n",
" \n",
"The exercise 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 exercise is split into 2 parts: the first part implements a PCA analysis on the image crops in order to reduce dimensionality, and the second part uses the outputs of the PCA and implement a SVM to classify the pixels into different land use classes.\n",
"\n",
" \n",
" \n",
" \n",
"
\n",
"Description of the dataset: \n",
" \n",
"Part 1: In the first part of the exercise, you will use the file s2_training_data.npz. It 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",
"Part 2: In the second part of the exercise, 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",
" Load the data using numpy.load(). Store the \"data\" in a variable named \"train_X\", anf the \"labels\" in a variable named \"train_Y\". Explore the shape of these arrays.\n",
"
\n",
" Use the function below to look at the (RGB-part) of individual crops. Take a look at the first crop of each class (remember: crops 1-5k contain forest, crops 5-10k contain fields, crops 10-15k are urban areas, and crops 15-20k are water).\n",
"
\n",
" The data train_X contains the raw values which span a large range. To compress them into a more “gaussian” shape, compute the logarithm of them.\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "2e300814",
"metadata": {},
"source": [
"### build function to compute principal components of image crops\n",
"The following tasks will implement a function \"compute_mean_PCs(X)\" which takes image crops as inputs, and returns the mean image, Eigen images (principal components reshaped into images), and Eigen values."
]
},
{
"cell_type": "markdown",
"id": "24f855a9",
"metadata": {},
"source": [
"#### compute mean crop, and mean-free crops"
]
},
{
"cell_type": "markdown",
"id": "aba7264f",
"metadata": {},
"source": [
"
\n",
" Compute the mean image crop, and name the resulting array \"mean\". \n",
" Subtract this \"mean\" from all crops, and name the resulting array \"mean_free\". \n",
" \n",
" Hint: The mean image should be an array of shape (4, 15, 15), and the mean free crops should be an array of shape (20000, 4, 15, 15)).\n",
"
\n",
" Compute the covariance matrix. This should be a matrix of size 900x900 (900 = 4 * 15 * 15). \n",
" \n",
" Start by defining a matrix of zeros with shape (900,900). Loop over all 20000 crops, and for each crop compute the outer product np.outer(u, v) of the mean-free vector with itself (i.e., u = v = vectorized crop). Accumulate (sum) the outer products together and finally divide by the total number of used crops. (To speed up development, consider using only every 10th crop). \n",
" \n",
" NB: the outer product $u \\otimes v$ is equivalent to a matrix multiplication $u \\cdot v^T$\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "83ffc370",
"metadata": {},
"source": [
"#### compute Eigen values & Eigen vectors of the covariance matrix"
]
},
{
"cell_type": "markdown",
"id": "07014a41",
"metadata": {},
"source": [
"
\n",
" Compute the Eigen values \"eig_val\" and Eigen vectors \"eig_vec\" of the covariance matrix using np.linalg.eig(). \n",
" \n",
" Sort them by importance (decreasing Eigen values) using the following code: \n",
"
\n",
" Reshape the 900 eigen vectors back into the images, i.e. with shape (900, 4, 15, 15), and name this variable \"principal_components\". \n",
" \n",
" Note that the matrix eig_vec contains the Eigen vectors in its rows (not columns as you might expect), so you'll need to transpose the matrix before reshaping it.\n",
"
\n",
" Merge steps 1.4.1 to 1.4.4 into a single function \"compute_mean_PCs\". \n",
" \n",
" The function will take the image crops as inputs (i.e., variable \"train_X\"), and will return the mean image (i.e., variable \"mean\"), the Eigen images (principal components reshaped into images, i.e., variable \"principal_components\"), and Eigen values (i.e., variable \"eig_val\").\n",
"
\n",
" Run \"compute_mean_PCs\" on the satellite data, and plot the returned Eigen values. \n",
" \n",
" Look at the first 64 principal components (or the RGB part thereof) using the following function:\n",
"
\n",
" Write a function \"compute_features(X, mean, principal_components, count)\" that takes crops \"X\", subtracts the mean, and projects them onto the first \"count\" principal components. \n",
" \n",
" The returned array should have a shape (X.shape[0], count), containing the coefficients (which we will use as features in next week's exercise). \n",
" \n",
" Hint: the projection is done using the dot product.\n",
"
\n",
" Write a function \"reconstruct_image(feature, mean, principal_components)\" that restores a crop given a feature/coefficient vector. \n",
" \n",
" Use the following code to compare, side by side, original image crops and reconstructions:\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",
"
\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",
"
\n",
" Train a support vector machine on the computed features and the corresponding labels. Use sklearn.svm.SVC(). \n",
" \n",
" Use default kernel (RBF), and set the following parameters: gamma=0.001, C=100.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 41,
"id": "270675e9",
"metadata": {},
"outputs": [],
"source": [
"# Options:\n",
"#\n",
"# kernel : {'linear', 'poly', 'rbf', 'sigmoid', 'precomputed'}, default='rbf'\n",
"#\n",
"# C : float, default=1.0\n",
"# Regularization parameter. The strength of the regularization is\n",
"# inversely proportional to C. Must be strictly positive. The penalty\n",
"# is a squared l2 penalty.\n",
"#\n",
"# gamma : {'scale', 'auto'} or float, default='scale'\n",
"# Kernel coefficient for 'rbf', 'poly' and 'sigmoid'.\n",
"# - if ``gamma='scale'`` (default) is passed then it uses\n",
"# 1 / (n_features * X.var()) as value of gamma,\n",
"# - if 'auto', uses 1 / n_features."
]
},
{
"cell_type": "markdown",
"id": "aa34039e",
"metadata": {},
"source": [
"### compute the training accuracy of the classifier"
]
},
{
"cell_type": "markdown",
"id": "c5f462d2",
"metadata": {},
"source": [
"
\n",
" Compute the accuracy on the training dataset using the function testAccuracy() defined above.\n",
"
\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",
"
\n",
" 1. Load the image 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",
" \n",
"
![](\"images/application_input.png\")