{"id":2614,"date":"2018-09-27T15:03:58","date_gmt":"2018-09-27T13:03:58","guid":{"rendered":"https:\/\/blogs.fu-berlin.de\/reseda\/?page_id=2614"},"modified":"2018-09-30T16:16:37","modified_gmt":"2018-09-30T14:16:37","slug":"prepare-samples-for-regression","status":"publish","type":"page","link":"https:\/\/blogs.fu-berlin.de\/reseda\/prepare-samples-for-regression\/","title":{"rendered":"Prepare Samples for Regression"},"content":{"rendered":"

We want to create training samples in the form of a data frame \u2013 which most regressors in R can handle easily. For this we will have to use a few GIS operations. Fortunately, these can also be done automatically in R, so that the entire workflow do not have to be performed manually again and again for other datasets. Finally, we will utilize the raster <\/a>and the rgeos <\/a>packages for this. <\/p>\n

The complete code is shown below. The shapefile reg_train_data.shp<\/em> and the imagery LC081930232017060201T1-SC20180613160412_subset.tif<\/em> can be downloaded here<\/a>.
\nHowever, we strongly advise you to have a look at the detailed description in the following.<\/p>\n

\r\nlibrary(raster)\r\nlibrary(rgeos)\r\n\r\n# import low resolution image and digitized polygon training areas\r\nsetwd(\"\/media\/sf_exchange\/regression\/\")\r\nimg <- brick("LC081930232017060201T1-SC20180613160412_subset.tif")\r\nshp <- shapefile("reg_train_data.shp")\r\n\r\n# rename bands if needed\r\nnames(img) <- c("b1", "b2", "b3", "b4", "b5", "b6", "b7")\r\n# create buffer with no width to avoid topology problems\r\nshp <- gBuffer(shp, byid=TRUE, width=0)\r\n\r\n# crop image to extent of training polygons\r\nimg.subset <- crop(img, shp)\r\n# extract pixel completely covered by training polygons\r\nimg.mask <- rasterize(shp, img.subset, getCover = TRUE)\r\nimg.mask[img.mask < 100] <- NA\r\nimg.subset <- mask(img.subset, img.mask)\r\n# transform those pixels to polygons\r\ngrid <- rasterToPolygons(img.subset)\r\n# create a new ID attribute for addressing\r\ngrid$ID <- seq.int(nrow(grid))\r\n\r\n# create new dataframe for final samples\r\nsmp <- data.frame()\r\n# do for each cell in grid:\r\nfor (i in 1:length(grid)) {\r\n  # intersect grid and shapefile for cell i\r\n  cell <- intersect(grid[i, ], shp)\r\n  # extract all target class polygons\r\n  cell <- cell[cell$Class_name == "building" | cell$Class_name == "impervious", ]\r\n  # if target class is in cell: calculate area, otherwise area is zero\r\n  if (length(cell) > 0) {\r\n    areaPercent <- sum(area(cell) \/ area(grid)[1])\r\n  } else {\r\n    areaPercent <- 0\r\n  }\r\n  # create new sample and append to training samples\r\n  newsmp <- cbind(grid@data[i, 1:nlayers(img)], areaPercent)\r\n  smp <- rbind(smp, newsmp)\r\n}\r\n<\/pre>\n
<\/br><\/a><\/p>\n
In-depth Guide<\/h1>\n
In order to use functionalities of the raster and rgeos packages, load them into your 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(\"rgeos\")\r\nlibrary(raster)\r\nlibrary(rgeos)\r\n<\/pre>\n
Next: set your working directory, in which all your image and shapefile data is stored by giving a character (do not forget the quotation marks ” “<\/span>) variable to setwd()<\/span>. Check your path with getwd()<\/span> and the stored files in it via dir()<\/span>:<\/p>\n
\r\nsetwd(\"\/media\/sf_exchange\/regression\/\")\r\n \r\ngetwd()\r\n## [1] \"\/media\/sf_exchange\/regression\"\r\n\r\ndir()  \r\n##  [1] \"LC081930232017060201T1-SC20180613160412_subset.tif\"                                                \r\n##  [2] \"reg_train_data.dbf\"     \r\n##  [3] \"reg_train_data.prj\"\r\n##  [4] \"reg_train_data.qpj\"\r\n##  [5] \"reg_train_data.shp\" \r\n##  [6] \"reg_train_data.shx\"   \r\n<\/pre>\n
If you do not get your files listed, there is an error in your working path – check again! Everything ready to go? Fine, then import your raster file as img<\/span> and your shapefile as shp<\/span> and have a look at them:<\/p>\n
\r\nimg <- brick("LC081930232017060201T1-SC20180613160412_subset.tif")\r\nimg\r\n## class       : RasterBrick \r\n## dimensions  : 504, 1030, 519120, 7  (nrow, ncol, ncell, nlayers)\r\n## resolution  : 30, 30  (x, y)\r\n## extent      : 369795, 400695, 5812395, 5827515  (xmin, xmax, ymin, ymax)\r\n## coord. ref. : +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 \r\n## data source : \/media\/sf_exchange\/regression\/LC081930232017060201T1-SC20180613160412_subset.tif \r\n## names       : LC0819302\/\/2_subset.1, LC0819302\/\/2_subset.2, LC0819302\/\/2_subset.3, LC0819302\/\/2_subset.4, LC0819302\/\/2_subset.5, LC0819302\/\/2_subset.6, LC0819302\/\/2_subset.7 \r\n## min values  :                 -2000,                 -2000,                 -1502,                 -1384,                  -107,                   -72,                   -27 \r\n## max values  :                 13240,                 13539,                 14574,                 14549,                 14243,                 14896,                 15388\r\n\r\nshp <- shapefile("reg_train_data.shp")\r\nshp\r\n## class       : SpatialPolygonsDataFrame \r\n## features    : 286 \r\n## extent      : 390233.4, 390698.4, 5817147, 5817720  (xmin, xmax, ymin, ymax)\r\n## coord. ref. : +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 \r\n## variables   : 1\r\n## names       : Class_name \r\n## min values  :  bare_soil \r\n## max values  :      water\r\n\r\ncompareCRS(shp, img)\r\n## [1] TRUE\r\n<\/pre>\n
Both the brick()<\/span> and shapefile()<\/span> functions are provided by the raster package. As shown above, they create objects of the class RasterBrick and SpatialPolygonsDataFrame respectively. The L8 raster provides 7 bands, and our example shapefile 286 features, i.e., polygons. You can check whether the projections of the two datasets are identical or not by executing compareCRS(shp, img)<\/span>. If this is not the case (output equals FALSE<\/em>), the raster package will automatically re-project your data on the fly later on. However, we recommend to adjust the projections manually in advance to prevent any future inaccuracies.
\nPlot your data to make sure everything is imported properly (check Visualize in R<\/a> for an intro to plotting):<\/p>\n
\r\nplot(shp)\r\n<\/pre>\n
\"\"<\/a><\/p>\n
With the argument add = TRUE<\/span> in line 2 several data layers can be displayed one above the other:<\/p>\n
\r\nplot(img[[4]], col=gray.colors(100))\r\nplot(shp, add=TRUE)\r\n<\/pre>\n
\"\"<\/a><\/p>\n
In Line 1, initially only the fourth channel of the image img[[4]]<\/span> was displayed with a gray color representation using 100 shades col=gray.colors(100)<\/span>. The black spots in the upper plot represent the polygons of the shapefile. It becomes clear that they cover only a very small part of our study area.<\/p>\n
Optional<\/strong>: Let us take a look at the naming of the raster bands via the names()<\/span> function. Those names can be quite bulky and cause problems in some illustrations when used as axis labels. You can easily rename it to something more concise by overriding the names with any string vector of the same length:<\/p>\n
\r\nnames(img)\r\n## [1] \"LC081930232017060201T1.SC20180613160412_subset.1\"\r\n## [2] \"LC081930232017060201T1.SC20180613160412_subset.2\"\r\n## [3] \"LC081930232017060201T1.SC20180613160412_subset.3\"\r\n## [4] \"LC081930232017060201T1.SC20180613160412_subset.4\"\r\n## [5] \"LC081930232017060201T1.SC20180613160412_subset.5\"\r\n## [6] \"LC081930232017060201T1.SC20180613160412_subset.6\"\r\n## [7] \"LC081930232017060201T1.SC20180613160412_subset.7\"\r\n \r\nnames(img) <- c("b1", "b2", "b3", "b4", "b5", "b6", "b7")\r\n \r\nnames(img)\r\n## [1] "b1" "b2" "b3" "b4" "b5" "b6" "b7"\r\n<\/pre>\n
In order to accelerate all computation processes in the following, we first trim the image to the extent of the shapefile using the crop()<\/span> function. It is often recommended to use a buffer of zero width wrapped around the polygons to avoid any topological errors in the upcoming process using gbuffer()<\/span> as shown in line 1:<\/p>\n
\r\nshp <- gBuffer(shp, byid=TRUE, width=0)\r\n\r\nimg.subset <- crop(img, shp)\r\n<\/pre>\n
Now we just want to extract the pixels that are 100% covered by the polygons. Full coverage is important to be able to extract reliable predictions about the ratio of an target class. In order to assess this coverage, we use the rasterize()<\/span> function with the argument getCover = TRUE<\/span> to generate a mask. In line 2, we set all pixels which have less than 100% coverage to NA<\/span>. In the end, we use the mask to subset our image data again using the mask()<\/span> function:<\/p>\n
\r\nimg.mask <- rasterize(shp, img.subset, getCover = TRUE)\r\nimg.mask[img.mask < 100] <- NA\r\nimg.subset <- mask(img.subset, img.mask)\r\n<\/pre>\n
Plot the image data again with the polygons in order to see our achievements:<\/p>\n
\r\nplot(shp)\r\nplot(img.subset[[4]], add=T)\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Only the raster cells within our polygons are left – great!
\nIn a next step, we transform those pixels to polygons, forming a grid with the information of the Landsat 8 bands maintained by using the  rasterToPolygons()<\/span> function. In addition, we want to assign each of these cells in the grid with an ID in order to address them (line 2):<\/p>\n
\r\ngrid <- rasterToPolygons(img.subset)\r\ngrid$ID <- seq.int(nrow(grid))\r\n\r\ngrid \r\n## class       : SpatialPolygonsDataFrame \r\n## features    : 93 \r\n## extent      : 390255, 390675, 5817165, 5817705  (xmin, xmax, ymin, ymax)\r\n## coord. ref. : +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 \r\n## variables   : 8\r\n## names       :   b1,   b2,   b3,   b4,   b5,   b6,   b7, ID \r\n## min values  :  144,  198,  339,  279,  678,  563,  398,  1 \r\n## max values  : 1175, 1245, 1452, 1422, 3686, 2382, 2193, 93\r\n\r\nplot(grid)\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Finally, we have generated 93 cells (polygons) in this grid<\/span> for which we now want to make a statement of how much of their areas is impervious, i.e., all polygons that have either building <\/em> or impervious<\/em> as Class_name<\/span> attribute.
\nFor this, we generate a new, empty data frame in which we write the training samples one after the other in line 1. We use a for-loop to iterate over each cell of the grid<\/span> in line 2:<\/p>\n
\r\nsmp <- data.frame()\r\nfor (i in 1:length(grid)) {\r\n  cell <- intersect(grid[i, ], shp)\r\n  cell <- cell[cell$Class_name == "building" | cell$Class_name == "impervious", ]\r\n  if (length(cell) > 0) {\r\n    areaPercent <- sum(area(cell) \/ area(grid)[1])\r\n  } else {\r\n    areaPercent <- 0\r\n  }\r\n  newsmp <- cbind(grid@data[i, 1:nlayers(img)], areaPercent)\r\n  smp <- rbind(smp, newsmp)\r\n}\r\n\r\nhead(smp)\r\n##      b1   b2   b3   b4   b5   b6  b7  areaPercent\r\n## 2  1175 1245 1452 1422 1742 1321 950 9.765630e-01\r\n## 3   497  597  736  767 1335 1175 947 7.986431e-01\r\n## 5   277  330  477  462 1221  960 743 3.669043e-01\r\n## 4   265  310  444  428 1077  789 545 1.727224e-01\r\n## 41  230  267  390  343  809  613 404 9.765625e-06\r\n## 42  168  217  339  312  678  563 398 3.540961e-02\r\n<\/pre>\n
Explanation<\/strong>:
\nLine 3: we select the ith cell of the grid and intersect this cell with the shp<\/span>, which gives us the geometries and attributes of both the shapefile and the cell.
\nLine 4: only the polygons belonging to the target classes were extracted, we discard all other polygons.
\nLine 5: if a polygon of a target class exists in the cell, do everything inside the curly brackets. If there is no polygon in the cell, jump to line 7.
\nLine 6: calculate the proportion of the area of each polygon and sum the proportions (area(grid)[1]<\/span> equals 900 m\u00b2 in our example due to the geometric resolution).
\nLine 7: if there is no polygon of class building <\/em>or impervious<\/em>, set areaPercent<\/span> to zero.
\nLine 10: use the column bind (cbind()<\/span>) function to combine the feature values (reflectances of Landsat image) and the areaPercent<\/span> value to form a new sample newsmp<\/span>.
\nLine 10: use row bind (rbind()<\/span>) to add the new sample newsmp<\/span> to our training dataset smp<\/span>.<\/p>\n
We now prepared a training data set smp<\/span> with all input features (Landsat 8 bands 1 to 7) and the percentage value of our target class (imperviousness<\/em>). <\/p>\n
Time to train a Support Vector Machine in the next section<\/a>!<\/p>\n
<\/br><\/br><\/p>\n
\n
\n