{"id":388,"date":"2018-05-03T14:36:26","date_gmt":"2018-05-03T12:36:26","guid":{"rendered":"https:\/\/blogs.fu-berlin.de\/reseda\/?page_id=388"},"modified":"2018-09-25T19:38:29","modified_gmt":"2018-09-25T17:38:29","slug":"area-adjusted-accuracies","status":"publish","type":"page","link":"https:\/\/blogs.fu-berlin.de\/reseda\/area-adjusted-accuracies\/","title":{"rendered":"Area Adjusted Accuracies"},"content":{"rendered":"<p>The area of land cover change obtained from a map may differ greatly from the true area of change because of misclassifications in your map. In regard of this, Olofsson et al. published two papers in <a href=\"https:\/\/ac.els-cdn.com\/S0034425712004191\/1-s2.0-S0034425712004191-main.pdf?_tid=c4857f5c-edcf-48ca-b372-d5a40a6e6005&amp;acdnat=1537182066_0c994600d448d5ef70ccbcf8d7465a6f\" rel=\"noopener\" target=\"_blank\">2013<\/a> and <a href=\"http:\/\/reddcr.go.cr\/sites\/default\/files\/centro-de-documentacion\/olofsson_et_al._2014_-_good_practices_for_estimating_area_and_assessing_accuracy_of_land_change.pdf\" rel=\"noopener\" target=\"_blank\">2014<\/a>, introducing an error-adjusted estimator of area, that can be easily produced once an accuracy assessment has been performed and an error matrix was constructed. Furthermore, they provide a confidence interval for the area of land change to quantify the uncertainty of the estimations. The following script reproduces the findings of <a href=\"http:\/\/reddcr.go.cr\/sites\/default\/files\/centro-de-documentacion\/olofsson_et_al._2014_-_good_practices_for_estimating_area_and_assessing_accuracy_of_land_change.pdf\" rel=\"noopener\" target=\"_blank\">Olofsson et al. 2014<\/a>. <\/p>\n<p><strong>What do you need?<\/strong><\/p>\n<ul>\n<li>your classification map<\/li>\n<li>your training shapefile<\/li>\n<li>your validation shapefile<\/li>\n<\/ul>\n<p><strong>What do you get?<\/strong><\/p>\n<ul>\n<li>error-adjusted area estimates of your classes<\/li>\n<li>confidence intervals for all accuracies<\/li>\n<\/ul>\n<p>The script is commented comprehensively and the equations from the paper are marked accordingly in the comments of this script.<\/p>\n<pre class=\"theme:amityreseda\">\r\n# area adjusted accuracy assessment \r\n# calculated as in Olofsson et al. 2014\r\n\r\nlibrary(raster)\r\n\r\n# import classification image and train and validation shapefiles\r\nsetwd(\"\/media\/sf_exchange\/landsatdata\/\")\r\nimg.class &lt;- raster(&quot;classification_RF.tif&quot;)\r\nshp.train &lt;- shapefile(&quot;training_data.shp&quot;)\r\nshp.valid &lt;- shapefile(&quot;validation_RF.shp&quot;)\r\n\r\n# create regular accuracy matrix \r\nconfmat &lt;- table(as.factor(extract(img.class, shp.valid)), as.factor(shp.valid$validclass))\r\n# get number of pixels per class and convert in km\u00b2\r\nimgVal &lt;- as.factor(getValues(img.class))\r\nnclass &lt;- length(unique(shp.train$classes))\r\nmaparea &lt;- sapply(1:nclass, function(x) sum(imgVal == x))\r\nmaparea &lt;- maparea * res(img.class)[1] ^ 2 \/ 1000000\r\n\r\n# set confidence interval\r\nconf &lt;- 1.96\r\n\r\n# total  map area\r\nA &lt;- sum(maparea)\r\n# proportion of area mapped as class i\r\nW_i &lt;- maparea \/ A\r\n# number of reference points per class\r\nn_i &lt;- rowSums(confmat) \r\n# population error matrix (Eq.4)\r\np &lt;- W_i * confmat \/ n_i\r\np[is.na(p)] &lt;- 0\r\n\r\n# area estimation\r\np_area &lt;- colSums(p) * A\r\n# area estimation confidence interval (Eq.10)\r\np_area_CI &lt;- conf * A * sqrt(colSums((W_i * p - p ^ 2) \/ (n_i - 1))) \r\n\r\n# overall accuracy (Eq.1)\r\nOA &lt;- sum(diag(p))\r\n# producers accuracy (Eq.2)\r\nPA &lt;- diag(p) \/ colSums(p)\r\n# users accuracy (Eq.3)\r\nUA &lt;- diag(p) \/ rowSums(p)\r\n\r\n# overall accuracy confidence interval (Eq.5)\r\nOA_CI &lt;- conf * sqrt(sum(W_i ^ 2 * UA * (1 - UA) \/ (n_i - 1)))\r\n# user accuracy confidence interval (Eq.6)\r\nUA_CI &lt;- conf * sqrt(UA * (1 - UA) \/ (n_i - 1)) \r\n# producer accuracy confidence interval (Eq.7)\r\nN_j &lt;- sapply(1:nclass, function(x) sum(maparea \/ n_i * confmat[ , x]) )\r\ntmp &lt;- sapply(1:nclass, function(x) sum(maparea[-x] ^ 2 * confmat[-x, x] \/ n_i[-x] * ( 1 - confmat[-x, x] \/ n_i[-x]) \/ (n_i[-x] - 1)) )\r\nPA_CI &lt;- conf * sqrt(1 \/ N_j ^ 2 * (maparea ^ 2 * ( 1 - PA ) ^ 2 * UA * (1 - UA) \/ (n_i - 1) + PA ^ 2 * tmp))\r\n\r\n# gather results\r\nresult &lt;- matrix(c(p_area, p_area_CI, PA * 100, PA_CI * 100, UA * 100, UA_CI * 100, c(OA * 100, rep(NA, nclass-1)), c(OA_CI * 100, rep(NA, nclass-1))), nrow = nclass)\r\nresult &lt;- round(result, digits = 2) \r\nrownames(result) &lt;- levels(as.factor(shp.train$classes))\r\ncolnames(result) &lt;- c(&quot;km\u00b2&quot;, &quot;km\u00b2\u00b1&quot;, &quot;PA&quot;, &quot;PA\u00b1&quot;, &quot;UA&quot;, &quot;UA\u00b1&quot;, &quot;OA&quot;, &quot;OA\u00b1&quot;)\r\nclass(result) &lt;- &quot;table&quot;\r\nresult\r\n<\/pre>\n<pre class=\"theme:amityreseda minimize:true\" title=\"Click here to show the script with intermediate results for comparison\">\r\nlibrary(raster)\r\n\r\nsetwd(\"\/media\/sf_exchange\/landsatdata\/\")\r\nimg.class &lt;- raster(&quot;classification_RF.tif&quot;)\r\nshp.train &lt;- shapefile(&quot;training_data.shp&quot;)\r\nshp.valid &lt;- shapefile(&quot;validation_RF.shp&quot;)\r\n\r\nconfmat &lt;- table(as.factor(extract(img.class, shp.valid)), as.factor(shp.valid$validclass))\r\nconfmat\r\n##    \r\n##      1  2  3  4  5  6\r\n##   1 31  0  7  2  5  0\r\n##   2  0 47  2  0  0  1\r\n##   3  0  7 39  0  4  0\r\n##   4  0  0  0 40  9  0\r\n##   5  0  0  1  2 47  0\r\n##   6  0  0  0  0  0 49\r\n\r\nimgVal &lt;- as.factor(getValues(img.class))\r\nhead(imgVal)\r\n## [1] 5 5 5 5 5 5\r\n## Levels: 1 2 3 4 5 6\r\n\r\nnclass &lt;- length(unique(shp.train$classes))\r\nnclass\r\n## [1] 6\r\n\r\nmaparea &lt;- sapply(1:nclass, function(x) sum(imgVal == x))\r\nmaparea\r\n## [1]  26741 121740  36345 129130 185222  19942\r\n\r\nmaparea &lt;- maparea * res(img.class)[1] ^ 2 \/ 1000000\r\nmaparea\r\n## [1]  24.0669 109.5660  32.7105 116.2170 166.6998  17.9478\r\n\r\nconf &lt;- 1.96\r\n\r\nA &lt;- sum(maparea)\r\nA\r\n## [1] 467.208\r\n\r\nW_i &lt;- maparea \/ A\r\nW_i\r\n## [1] 0.05151217 0.23451225 0.07001271 0.24874788 0.35679997 0.03841501\r\n\r\nn_i &lt;- rowSums(confmat) \r\nn_i\r\n##  1  2  3  4  5  6 \r\n## 45 50 50 49 50 49\r\n\r\np &lt;- W_i * confmat \/ n_i\r\np[is.na(p)] &lt;- 0\r\nround(p, digits = 4)\r\n##    \r\n##          1      2      3      4      5      6\r\n##   1 0.0355 0.0000 0.0080 0.0023 0.0057 0.0000\r\n##   2 0.0000 0.2204 0.0094 0.0000 0.0000 0.0047\r\n##   3 0.0000 0.0098 0.0546 0.0000 0.0056 0.0000\r\n##   4 0.0000 0.0000 0.0000 0.2031 0.0457 0.0000\r\n##   5 0.0000 0.0000 0.0071 0.0143 0.3354 0.0000\r\n##   6 0.0000 0.0000 0.0000 0.0000 0.0000 0.0384\r\n\r\np_area &lt;- colSums(p) * A\r\np_area\r\n##         1         2         3         4         5         6 \r\n##  16.57942 107.57151  36.97457 102.60865 183.33473  20.13912\r\n\r\np_area_CI &lt;- conf * A * sqrt(colSums((W_i * p - p ^ 2) \/ (n_i - 1))) \r\np_area_CI\r\n##         1         2         3         4         5         6 \r\n##  3.292170  7.948700  9.994001 15.744342 17.208163  4.294987\r\n\r\nOA &lt;- sum(diag(p))\r\nOA\r\n## [1] 0.8874041\r\n\r\nPA &lt;- diag(p) \/ colSums(p)\r\nPA\r\n##         1         2         3         4         5         6 \r\n## 1.0000000 0.9574286 0.6900470 0.9245908 0.8547088 0.8911909\r\n\r\nUA &lt;- diag(p) \/ rowSums(p)\r\nUA\r\n##         1         2         3         4         5         6 \r\n## 0.6888889 0.9400000 0.7800000 0.8163265 0.9400000 1.0000000\r\n\r\nOA_CI &lt;- conf * sqrt(sum(W_i ^ 2 * UA * (1 - UA) \/ (n_i - 1)))\r\nOA_CI\r\n## [1] 0.04079462\r\n\r\nUA_CI &lt;- conf * sqrt(UA * (1 - UA) \/ (n_i - 1)) \r\nUA_CI\r\n##          1          2          3          4          5          6 \r\n## 0.13679244 0.06649632 0.11598896 0.10954451 0.06649632 0.00000000\r\n\r\nN_j &lt;- sapply(1:nclass, function(x) sum(maparea \/ n_i * confmat[ , x]) )\r\ntmp &lt;- sapply(1:nclass, function(x) sum(maparea[-x] ^ 2 * confmat[-x, x] \/ n_i[-x] * ( 1 - confmat[-x, x] \/ n_i[-x]) \/ (n_i[-x] - 1)) )\r\nPA_CI &lt;- conf * sqrt(1 \/ N_j ^ 2 * (maparea ^ 2 * ( 1 - PA ) ^ 2 * UA * (1 - UA) \/ (n_i - 1) + PA ^ 2 * tmp))\r\nPA_CI\r\n##          1          2          3          4          5          6 \r\n## 0.00000000 0.02843232 0.17545912 0.08399238 0.06198829 0.19006061\r\n\r\nresult &lt;- matrix(c(p_area, p_area_CI, PA * 100, PA_CI * 100, UA * 100, UA_CI * 100, c(OA * 100, rep(NA, nclass-1)), c(OA_CI * 100, rep(NA, nclass-1))), nrow = nclass)\r\nresult &lt;- round(result, digits = 2) \r\nrownames(result) &lt;- levels(as.factor(shp.train$classes))\r\ncolnames(result) &lt;- c(&quot;km\u00b2&quot;, &quot;km\u00b2\u00b1&quot;, &quot;PA&quot;, &quot;PA\u00b1&quot;, &quot;UA&quot;, &quot;UA\u00b1&quot;, &quot;OA&quot;, &quot;OA\u00b1&quot;)\r\nclass(result) &lt;- &quot;table&quot;\r\nresult\r\n##              km\u00b2   km\u00b2\u00b1     PA    PA\u00b1     UA    UA\u00b1     OA    OA\u00b1\r\n## baresoil   16.58   3.29 100.00   0.00  68.89  13.68  88.74   4.08\r\n## forest    107.57   7.95  95.74   2.84  94.00   6.65              \r\n## grassland  36.97   9.99  69.00  17.55  78.00  11.60              \r\n## urban_hd  102.61  15.74  92.46   8.40  81.63  10.95              \r\n## urban_ld  183.33  17.21  85.47   6.20  94.00   6.65              \r\n## water      20.14   4.29  89.12  19.01 100.00   0.00\r\n<\/pre>\n","protected":false},"excerpt":{"rendered":"<p>The area of land cover change obtained from a map may differ greatly from the true area of change because of misclassifications in your map. In regard of this, Olofsson et al. published two papers in 2013 and 2014, introducing an error-adjusted estimator of area, that can be easily produced once an accuracy assessment has &hellip; <a href=\"https:\/\/blogs.fu-berlin.de\/reseda\/area-adjusted-accuracies\/\" class=\"more-link\">Continue reading<span class=\"screen-reader-text\"> &#8220;Area Adjusted Accuracies&#8221;<\/span><\/a><\/p>\n","protected":false},"author":3237,"featured_media":0,"parent":0,"menu_order":0,"comment_status":"closed","ping_status":"closed","template":"","meta":{"footnotes":""},"class_list":["post-388","page","type-page","status-publish","hentry"],"_links":{"self":[{"href":"https:\/\/blogs.fu-berlin.de\/reseda\/wp-json\/wp\/v2\/pages\/388","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/blogs.fu-berlin.de\/reseda\/wp-json\/wp\/v2\/pages"}],"about":[{"href":"https:\/\/blogs.fu-berlin.de\/reseda\/wp-json\/wp\/v2\/types\/page"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.fu-berlin.de\/reseda\/wp-json\/wp\/v2\/users\/3237"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.fu-berlin.de\/reseda\/wp-json\/wp\/v2\/comments?post=388"}],"version-history":[{"count":11,"href":"https:\/\/blogs.fu-berlin.de\/reseda\/wp-json\/wp\/v2\/pages\/388\/revisions"}],"predecessor-version":[{"id":2535,"href":"https:\/\/blogs.fu-berlin.de\/reseda\/wp-json\/wp\/v2\/pages\/388\/revisions\/2535"}],"wp:attachment":[{"href":"https:\/\/blogs.fu-berlin.de\/reseda\/wp-json\/wp\/v2\/media?parent=388"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}