---
title: "RasterMaster"
output:
html_document:
toc: TRUE
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.align = 'center', warning = FALSE)
library(rgdal)
library(rgeos)
library(raster)
library(rasterVis)
library(ggplot2)
library(dplyr)
library(readr)
library(RStoolbox)
MAXPIXELS <- 250000
```
*tl;dr: If you use the `raster` package in R the image alignment with the polygons should be easy.*
### Rasters
In this kernel, we show to plot the rasters with the `raster` package and with `ggplot` plus `rasterVis`.
We also show how to align the rasters with the polygons in R if you are using the `raster` classes.
We'll be using the image 6060_2_3, which is in the training set, so we have polygons for it.
```{r data.paths, echo=FALSE}
img1.path <- '../input/three_band/6060_2_3.tif'
bandM.path <- '../input/sixteen_band/6060_2_3_M.tif'
bandP.path <- '../input/sixteen_band/6060_2_3_P.tif'
bandA.path <- '../input/sixteen_band/6060_2_3_A.tif'
zip.path <- '../input/three_band.zip'
gold.path <- '../input/train_wkt_v4.csv'
```
### Basic Load and Plot
We load the data as a `raster::brick`, which is a data-aligned stack of rasters.
```{r load.img}
img1 <- brick(img1.path)
names(img1) <- c('R','G','B')
img1
cellStats(img1, summary)
```
The default plot method plots the bands separately.
```{r base.plot}
plot(img1, maxpixels=MAXPIXELS)
```
There is also a `plotRGB` function in package `raster` for visualization.
The images come out best with the parameter `stretch='lin'`
```{r RBGplot.img}
plotRGB(img1, stretch='lin', maxpixels=MAXPIXELS)
```
### Using ggplot + rasterVis
The `rasterVis` package allows us to plot rasters in `ggplot` using `geom_tile`.
It uses a wrapper function called `gplot` to start things off.
*NOTE: There is only one "g" in gplot*.
```{r ggplot}
g <- gplot(img1) + geom_tile(aes(fill=value)) + facet_wrap(~variable)
g <- g + scale_fill_gradient(low = 'white', high = 'blue') + coord_equal()
g
```
### ggplot + RStoolbox
*Hat Tip to Chippy for the pointer to the `RStoolbox` package used here*
It looks like this is how we get a 3-band image as a ggplot, using the `RStoolbox` package.
```{r ggRGB}
g <- ggRGB(img1, r=1, g=2, b=3, maxpixels=MAXPIXELS, stretch='lin')
g
```
### Histograms
The `rasterVis` package can take quick histograms by channel of a raster object,
including stack/brick.
```{r histogram}
histogram(img1)
```
On smaller regions, this might aid in classification.
### Labeled Polygons
We have ground truth polygons for 25 images.
The labels are coded 1-10 as follows:
1. Building
2. Misc. Structures
3. Roads
4. Paths
5. Trees
6. Crops
7. Waterway
8. Standing Water
9. Large Vehicle
10. Small Vehicle
```{r gold}
# readr::read_csv tries to read the ID's as date w/o the column spec.
gold.data <- read_csv(gold.path, col_types='cic')
gold <- filter(gold.data, MultipolygonWKT != 'MULTIPOLYGON EMPTY',
ImageId=='6060_2_3')
glimpse(gold)
```
### Plotting Polygons
Only 6 classes are present in this image.
We'll look at 3 classes that take up some area, using the default plots.
```{r plot.polygons}
cropsWKT <- gold$MultipolygonWKT[gold$ClassType==6]
crops.shape <- readWKT(cropsWKT)
plot(crops.shape, col='lightblue')
treesWKT <- gold$MultipolygonWKT[gold$ClassType==5]
trees.shape <- readWKT(treesWKT)
plot(trees.shape, col='green')
bldgWKT <- gold$MultipolygonWKT[gold$ClassType==1]
bldg.shape <- readWKT(bldgWKT)
plot(bldg.shape, col='red')
plot(crops.shape, col='lightblue')
plot(trees.shape, col='green', add=T)
plot(bldg.shape, col='red', add=T)
```
In this image crops (class 6) take up most of the area.
### Polygons in ggplot
This recipe plots the polygons in `ggplot`.
```{r ggpolygons}
g <- ggplot(crops.shape, aes(x=long, y=lat, group=group))
g + geom_polygon(fill='lightblue') + coord_equal()
g <- ggplot(trees.shape, aes(x=long, y=lat, group=group))
g + geom_polygon(fill='green') + coord_equal()
g <- ggplot(bldg.shape, aes(x=long, y=lat, group=group))
g + geom_polygon(fill='red', color='black') + coord_equal()
g <- ggplot(crops.shape, aes(x=long, y=lat, group=group))
g <- g + geom_polygon(fill='lightblue')
g <- g + geom_polygon(data=trees.shape, fill='darkgreen')
g <- g + geom_polygon(data=bldg.shape, fill='red', color='black')
g + coord_equal()
```
This might be one of those rare cases where the default plots are better than
those from `ggplot`.
### Aligning the Rasters and Polygons
In R, using the `raster` package, there is no real alignment problem.
We just set the extent of the raster to the values from `grid_sizes.csv`.
This changes the scale of the rasters to the scale the polygons are on.
That is also the scale of the submission.
For this image, we have $xmax=0.009188$ and $ymin=-0.00904$ (from grid_sizes.csv).
```{r alignment}
xmin(img1) <- 0
xmax(img1) <- 0.009188
ymin(img1) <- -0.00904
ymax(img1) <- 0
plotRGB(img1, stretch='lin', maxpixels=MAXPIXELS)
plot(trees.shape, add=TRUE, col='lightgreen')
```
The paths (class 4) will really test this.
```{r re-alignment}
pathsWKT <- gold$MultipolygonWKT[gold$ClassType==4]
paths.shape <- readWKT(pathsWKT)
plotRGB(img1, stretch='lin', maxpixels=MAXPIXELS)
plot(paths.shape, add=TRUE, col='red')
```
Paths are hard to distinguish, but as far as I can tell, that looks right.
### Aligning Polygons and Rasters using ggplot
Aligning the polygons over the rasters using `ggplot` just involves combining earlier examples.
```{r ggalign}
g <- ggRGB(img1, r=1, g=2, b=3, maxpixels=MAXPIXELS, stretch='lin')
g <- g + geom_polygon(data=trees.shape, aes(x=long, y=lat, group=group), color='black', fill='green', alpha=0.5, inherit.aes=FALSE)
g + coord_equal()
```
### Sixteen Band Data
Finally, we'll take a quick look at the 16-band data.
```{r load16band, results='hold'}
bandM <- brick(bandM.path)
bandP <- brick(bandP.path)
bandA <- brick(bandA.path)
slotNames(bandM) # The others are the same
c(bandM=nlayers(bandM), bandP=nlayers(bandP), bandA=nlayers(bandA))
```
We'll plot them separately since they don't make any obvious composite image.
```{r plot16band}
plot(bandM, maxpixels=MAXPIXELS)
plot(bandP, maxpixels=MAXPIXELS)
plot(bandA, maxpixels=MAXPIXELS)
```