{"id":1964,"date":"2018-07-09T15:23:21","date_gmt":"2018-07-09T13:23:21","guid":{"rendered":"https:\/\/blogs.fu-berlin.de\/reseda\/?page_id=1964"},"modified":"2018-09-30T17:04:48","modified_gmt":"2018-09-30T15:04:48","slug":"svm-classification","status":"publish","type":"page","link":"https:\/\/blogs.fu-berlin.de\/reseda\/svm-classification\/","title":{"rendered":"SVM Classification"},"content":{"rendered":"

Just as with the Random Forest, there are quite a few R packages that provide SVM, e.g., caret<\/a>, e1071<\/a>, or kernLab<\/a>. We will use the e1071 package, as it offers an interface to the well-known libsvm <\/a>implementation. <\/p>\n

Below you can see a complete code implementation. While this is already executable with your input data, you should read the comprehensive in-depth guide in the following to understand the code in detail.<\/p>\n

\r\n# import packages\r\nlibrary(raster)\r\nlibrary(e1071)\r\n \r\n# import image (img) and shapefile (shp)\r\nsetwd(\"\/media\/sf_exchange\/landsatdata\/\")\r\nimg <- brick("LC081930232017060201T1-SC20180613160412_subset.tif")\r\nshp <- shapefile("training_data.shp")\r\n \r\n# extract samples with class labels and put them all together in a dataframe\r\nnames(img) <- c("b1", "b2", "b3", "b4", "b5", "b6", "b7")\r\nsmp <- extract(img, shp, df = TRUE)\r\nsmp$cl <- as.factor( shp$classes[ match(smp$ID, seq(nrow(shp)) ) ] )\r\nsmp <- smp[-1]\r\n \r\n# optional: shuffle samples and undersample if necessary\r\nsmp <- smp[sample(nrow(smp)),]\r\nsmp <- smp[ave(1:(nrow(smp)), smp$cl, FUN = seq) <= min(summary(smp$cl)), ]\r\n \r\n# tune and train model via gridsearch (gamma & cost)\r\ngammas = 2^(-8:5)\r\ncosts = 2^(-5:8)\r\nsvmgs <- tune(svm,\r\n              train.x = smp[-ncol(smp)],\r\n              train.y = smp$cl,\r\n              scale = TRUE,\r\n              kernel = "radial", \r\n              type = "C-classification",\r\n              ranges = list(gamma = gammas,\r\n                            cost = costs),\r\n              tunecontrol = tune.control(cross = 5)\r\n              )\r\n\r\n# save svmmodel\r\nsvmmodel <- svmgs$best.model\r\nsave(svmmodel, file = "svmmodel.RData")\r\n \r\n# predict image data with svm model\r\nresult <- predict(img,\r\n                  svmmodel,\r\n                  filename = "classification_svm.tif",\r\n                  overwrite = TRUE\r\n                  )\r\n<\/pre>\n
<\/br><\/a><\/p>\n
In-depth Guide<\/font><\/h1>\n
In order to be able to use the functions of the e1071 package, we must additionally load the library into the current session via library()<\/span>. If you do not use our VM, you must first download and install the packages with install.packages()<\/span>:<\/p>\n
\r\n#install.packages(\"raster\")\r\n#install.packages(\"e1071\")\r\nlibrary(raster)\r\nlibrary(e1071)\r\n<\/pre>\n
First, it is necessary to process the training samples in the form of a data frame. The necessary steps are described in detail in the previous section<\/a>.<\/p>\n
\r\nnames(img) <- c("b1", "b2", "b3", "b4", "b5", "b6", "b7")\r\nsmp <- extract(img, shp, df = TRUE)\r\nsmp$cl <- as.factor( shp$classes[ match(smp$ID, seq(nrow(shp)) ) ] )\r\nsmp <- smp[-1]\r\n<\/pre>\n
After that, you can identify the number of available training samples per class with summary()<\/span>. There is often an imbalance in the number of those training pixels, i.e., one class is represented by a large number of pixels, while another class has very few samples:<\/p>\n
\r\nsummary(smp$cl)\r\n##  baresoil    forest grassland  urban_hd  urban_ld     water \r\n##       719      2074      1226      1284       969       763\r\n<\/pre>\n
This often leads to the problem that classifiers favor and overclass strongly-represented classes in the classification. Unfortunately, a SVM can not handle this problem as elegantly as Random Forest<\/a> methods, why you might need to prepare your training data:
\nOnly if your data is imbalanced<\/strong>, you should consider to do the following undersampling: We shuffle our entire training dataset and pick out the same number of samples for each class. The purpose of shuffle is to prevent sampling of all pixels that are next to each other (spatial autocorrelation). The training data is mixed using the sample()<\/span> method. We sample all rows of our dataset nrow(smp)<\/span> in a random manner in line 10:<\/p>\n
\r\nhead(smp)\r\n##    b1  b2  b3  b4  b5  b6 b7    cl\r\n## 1 192 229 321 204 161 100 72 water\r\n## 2 179 203 233 130 156  93 68 water\r\n## 3 189 221 272 164 173 109 81 water\r\n## 4 159 179 188  97 125  63 40 water\r\n## 5 171 194 203 116 146  82 56 water\r\n## 6 164 188 196 108 138  71 51 water\r\n\r\nsmp <- smp[sample(nrow(smp)), ]\r\n\r\nhead(smp)\r\n##       b1  b2  b3  b4   b5   b6   b7       cl\r\n## 2024 176 195 324 234 2146 1022  479   forest\r\n## 5545 409 517 844 988 3038 3005 1818 baresoil\r\n## 2877 195 212 363 293 2500 1435  690   forest\r\n## 6210 303 322 482 411 2465 1410  872 urban_hd\r\n## 6613 367 418 672 600 2855 1876 1210 urban_ld\r\n## 321  163 178 201 136  124   76   56    water\r\n<\/pre>\n
Compare the first six entries of the smp<\/span> dataset before and after the command in line 10 using the head()<\/span> function. You will notice that the order of the rows is shuffled!
\nWith the min()<\/span> function, we can select the minimum value from the class column’s summary and use the ave()<\/span> function to make a selection of the same size for each class:<\/p>\n
\r\nsummary(smp$cl)\r\n##  baresoil    forest grassland  urban_hd  urban_ld     water \r\n##       719      2074      1226      1284       969       763\r\n\r\nsmp.maxsamplesize <- min(summary(smp$cl))\r\nsmp.maxsamplesize\r\n## [1] 719\r\n\r\nsmp <- smp[ave(1:(nrow(smp)), smp$cl, FUN = seq) <= smp.maxsamplesize, ]\r\n\r\nsummary(smp$cl)\r\n##  baresoil    forest grassland  urban_hd  urban_ld     water \r\n##       719       719       719       719       719       719\r\n<\/pre>\n
Let the training begin!
\nClassifications using an SVM require two parameters: one gamma \\(\\gamma\\) and one cost \\(C\\) value (for more details refer to the SVM section<\/a>). These hyperparameters significantly determine the performance of the model. Finding the best hyparameters is not trivial and the best combination can not be determined in advance. Thus, we try to find the best combination by trial and error. Therefore, we create two vectors comprising all values that should be tried out:<\/p>\n
\r\ngammas = 2^(-8:5)\r\ngammas\r\n##  [1]  0.00390625  0.00781250  0.01562500  0.03125000  0.06250000\r\n##  [6]  0.12500000  0.25000000  0.50000000  1.00000000  2.00000000\r\n## [11]  4.00000000  8.00000000 16.00000000 32.00000000\r\n\r\ncosts = 2^(-5:8)\r\ncosts\r\n##  [1]   0.03125   0.06250   0.12500   0.25000   0.50000   1.00000   2.00000\r\n##  [8]   4.00000   8.00000  16.00000  32.00000  64.00000 128.00000 256.00000\r\n<\/pre>\n
So we have 14 different values for \\(\\gamma\\) and 14 different values for \\(C\\). Thus, the whole training process is done for 196 (14 * 14) models, often referred to as gridsearch<\/strong>. Conversely, this means that the more parameters we check, the longer the training process takes.
\nWe start the training with the tune()<\/span> function. We need to specify the training samples as train.x<\/span>, i.e., all columns of our smp<\/span> dataframe except the last one, and the corresponding class labels as train.y<\/span>, i.e. the last column of our smp<\/span> dataframe:<\/p>\n
\r\nsvmgs <- tune(svm,\r\n              train.x = smp[-ncol(smp)],\r\n              train.y = smp[ ncol(smp)],\r\n              type = "C-classification",\r\n              kernel = "radial", \r\n              scale = TRUE,\r\n              ranges = list(gamma = gammas, cost = costs),\r\n              tunecontrol = tune.control(cross = 5)\r\n              )\r\n<\/pre>\n
We have to set the type<\/span> of the SVM to “C-classification”<\/span> in line 4 in order to perform a classification task. Furthermore, we set the kernel used in training and predicting to a RBF kernel<\/a> via “radial”<\/span>. We set the argument scale<\/span> to TRUE<\/span> in order to initiate the z-transformation of our data<\/a>. The argument ranges<\/span> in line 7 takes a named list of parameter vectors spanning the sampling range. We put our gammas<\/span> and costs<\/span> vectors in this list. By using the tunecontrol<\/span> argument in line 8, you can set k for the k-fold cross validation on the training data<\/a>, which is necessary to assess the model performance. <\/p>\n
Depending on the complexity of the data, this step may take some time. Once completed, you can check the output by calling the resultant object name:<\/p>\n
\r\nsvmgs\r\n## \r\n## Parameter tuning of 'svm':\r\n## \r\n## - sampling method: 5-fold cross validation \r\n## \r\n## - best parameters:\r\n##  gamma cost\r\n##      2    8\r\n## \r\n## - best performance: 0.02132796\r\n<\/pre>\n
In the course of the cross-validation, the overall accuracies were compared and the best parameters were determined: In our example, those are 2 and 8 for \\(\\gamma\\) and \\(C\\), respectively. Furthermore, the error of the best model is displayed: 2.1% error rate (which is quite good!).
\nWe can plot the errors of all 196 different hyperparameter combinations in the so called gridsearch using the standard plot()<\/span> function:<\/p>\n
\r\nplot(svmgs)\r\n<\/pre>\n
\"\"<\/a><\/p>\n
The plot appears very homogeneous, as the performance varies only very slightly and is already at a very high level. It can happen that the highest accuracies are found on the edge of the grid search. Then the training with other ranges for \\(\\gamma\\) and \\(C\\) should be done again.
\nWe can extract the best model out of our svmgs<\/span> to use for image prediction:<\/p>\n
\r\nsvmmodel <- svmgs$best.model\r\nsvmmodel\r\n## \r\n## Call:\r\n## best.tune(method = svm, train.x = smp[-ncol(smp)], train.y = smp$cl, \r\n##     ranges = list(gamma = gammas, cost = costs), tunecontrol = tune.control(cross = 5), \r\n##     scale = TRUE, kernel = "radial", type = "C-classification")\r\n## \r\n## \r\n## Parameters:\r\n##    SVM-Type:  C-classification \r\n##  SVM-Kernel:  radial \r\n##        cost:  8 \r\n##       gamma:  2 \r\n## \r\n## Number of Support Vectors:  602\r\n<\/pre>\n
Save the best model by using the save()<\/span> function. This function saves the model object svmmodel<\/span> to your working directory, so that you have it permanently stored on your hard drive. If needed, you can load it any time with load()<\/span>.<\/p>\n
\r\nsave(svmmodel, file = \"svmmodel.RData\")\r\n#load(\"svmmodel.RData\")\r\n<\/pre>\n
Since your SVM model is now completely trained, you can use it to predict all the pixels in your image. The command method predict()<\/span> takes a lot of work from you: It is recognized that there is an image which will be processed pixel by pixel. As with the training pixels, each image pixel is now individually classified and finally reassembled into your final classification image. Use the argument filename = <\/span> to specify the name of your output map:<\/p>\n
\r\nresult <- predict(img,\r\n                  svmmodel,\r\n                  filename = "classification_svm.tif",\r\n                  overwrite = TRUE\r\n                  )\r\n\r\nplot(result, \r\n     axes = FALSE, \r\n     box = FALSE,\r\n     col = c(\"#fbf793\", # baresoil\r\n             \"#006601\", # forest\r\n             \"#bfe578\", # grassland\r\n             \"#d00000\", # urban_hd\r\n             \"#fa6700\", # urban_ld\r\n             \"#6569ff\"  # water\r\n             )\r\n     )\r\n<\/pre>\n
Support Vector Machine classification completed!
\nIf you store the result of the predict()<\/span> function in an object, e.g., result<\/span>, you can plot the map using the standard plot command and passing this object:<\/p>\n
\"\"<\/a><\/p>\n
<\/br><\/br><\/p>\n
\n
\n