Classifying High Resolution Aerial Imagery - Part 2

I have been attempting to use random forests to classify high resolution aerial imagery. Part one of this post series was my first attempt. The aerial imagery dataset that I am working on is made up of many ortho tiles that I need to classify into vegetation categories. The first attempt was to classify vegetation on one tile. This note documents classifying vegetation across tiles.

See the r script below.

Notes

Set Up

The training polygon shapefile consists of 8 classes (it should be 7 but I accidentally used two names for bare ground: "BG_Soil" and "BG_Rock").

The Results

The out of bag (OOB) estimate of error rate was 1.85% which is pretty good. However, the training dataset has very imbalanced classes.

Classes were very imbalanced in the training dataset. Both "BG_Soil" and "BG_Rock" should be combined making soil account for over 50% of the area. {: .caption}

You can see the class imbalance in the confusion matrix as well.

The most difficult divisions for the random forests was between "Black_Sage" and "Sage". Which makes a lot of sense. About 6.3% of the time Sage was classified as Black_Sage, out of a total error of 10%, and about 8.7% of the time "Black_Sage" was classified as Sage, out of a total error of 16%.

Making Improvements (Iterate or Die)

Class Size: Obviously, I need to improve the class balance. I think the best way to do this is to both draw more training polygons in those classes with less (black_sage, sage and other_shrub) coverage and two add a balancing step in the script. Which brings me to sampsize().

Sample Size I also need to take advantage of randomForest's ability to manage sample size with sampsize(). I have had a hard time automating the sample size to be used in the function so that I don't have to find the smallest class and then manually put in that number as many times as there are classes. For instance if I had five classes and the minimum class size was 500 pixels I would need to set sampsize to be sampsize(500,500,500,500,500). Or I could sample the combined training set by randomly sampling pixels before I run random forests.

Segmentation - I still think segmentation has a role to play in this classification process, but my first attempt at segmentation was a complete failure. I need to find a tutorial on classification before I can make any improvement.

Adding Data - Using only five variables, red, blue, green, alpha, and ndvi, seems like very little data to classify imagery. There are a few datasets that I am considering adding:

• Soils data
• Elevation type data (slope, flow direction, terrain roughness, etc. )
• Landfire data, which is essentially landsat data that has been classified.

But if this process is going to scale to the rest of the state, the datasources that I use need to be available to everyone.

Scripts

## Load libraries
require(dplyr)
require(stringr)
require(raster)
require(rgdal)
require(rgeos)
require(randomForest)



## Function to work subset training dataset of one tile
# - inputs (all are automatically created by build_model):
# -- tile_name: is the name of the tile from a single iteration in the list
# -- folder_of_tiles: is the directory that all of the tiles are stored.
# -- file_path_to_shape: the path to the shapefile that will classify.
# - steps:
# -- RASTER STEPS
# -- read raster brick
# -- change names to b1 to b4 for each band
# -- create an NDVI layer
# -- combine NDVI to four Bands
# -- renmae each band
# -- SHAPEFILE STEPS
# -- split file path name
# -- read shapefile with split filepath
# -- Don't know what the layer paste is all about, I think I can take it out.
# -- convert Class, the supervised classes for training, to numeric
# -- clip the shapefile to the extent of the loaded tile.
# -- rasterized the clipped shapefile with values being class.
# -- rename heading to class
# -- remove temporary variables
# -- MORE RASTER STEPS
# -- clip the raster to the shapefile extent
# -- add training raster (converted shapefile) to the ortho raster data.
# -- convert raster to values (table)
# -- convert values to data.frame, the above step may be redundant.
# -- filter na's, because when you clip a raster to a geometry, the pixels outside of the geometry are still there but they have values of NA.  This step is the most memory intensive step in the whole process.
# -- save the RDS to be modeled in random forests. There is definitely a better way to do this.....
# -- rm data_to_be_modeled otherwize each iteration will fail on the remove NA step.


subset_training_data_from_tile<- function( tile_name, folder_of_tiles, file_path_to_shape){


  ## read raster make ndvi layer and add it to the brick.
  raster<-brick(paste0(paste(folder_of_tiles, tile_name ,sep="/"), ".tif" ))
  names(raster)<-c("b1", "b2", "b3", "b4") ## Rename Bands for ease of use
  ndvi<-overlay(raster$b4, raster$b3, fun=function(x,y){(x-y)/(x+y)})
  cov<-addLayer(raster,ndvi) ## Merge NDVI and the tile imagery
  names(cov)<-c("b1", "b2", "b3", "b4", "ndvi")

  ## parse and read shapefile
  temp<-stringr::str_split(file_path_to_shape, "/")
  temp2<-readOGR(dsn=paste(head(temp[[1]], -1) , collapse = "/"), layer=paste(tail(temp[[1]], 1))) ## Don't remember what this does?
  temp2@data$code<-as.numeric(temp2@data$Class)
  shape<-crop(temp2, extent(cov))%>%
    rasterize(ndvi, field="code")  
  names(shape)<-"class"

  rm(ndvi, temp, temp2, raster)

  data_to_be_modeled<-crop(cov, extent(shape))%>%
    mask(shape)%>%
    addLayer(shape)%>%
    getValues()%>%
    as.data.frame()%>%
    filter(!is.na(class))

  saveRDS(data_to_be_modeled, paste0("training_datasets/", tile_name, "training_data.rds"))

  rm(data_to_be_modeled)
}


## Run random forest on combined training dataset
# - reads the folder with tiles
# -- takes all file names and stores them in a list to be transfered to the `subset_training_data_from_tile` function which
# - reads the shapefile
# -- adds the shapefile to `subset_training_data_from_tile`
# - tile names are looped over and sent to `subset_training_data_from_tile` individually with the shapefile.
# - all individually returned tile training datasets are combined together with rbind.fill
# - combined dataset, is run through random forests.
build_model<-function(folder_of_tiles, file_path_to_shape){
  t<-base::list.files(folder_of_tiles)%>%
    tools::file_path_sans_ext()%>%
    tools::file_path_sans_ext()%>%
    base::unique() ## May want to change this to just a vector and use this for just apply

  t2<-folder_of_tiles
  t3<-file_path_to_shape

  lapply(t, subset_training_data_from_tile, t2, t3)

  t4<-base::list.files("training_datasets")

  t5<-lapply(t4, function(x)readRDS(paste0("training_datasets/", x)))
  t6 <- do.call(plyr::rbind.fill, t5)

  head(t6)

  rm(t5)

  fit<-randomForest(as.factor(class)~.,
                    data=t6,
                    importance=TRUE,
                    ntree=500, norm.votes = FALSE
  )
  rm(t6)
  return(fit)
}


## Run the function `build_model`
set.seed(420)
start_time <- Sys.time()

test_model<-build_model("test_tiles", "test_shape/training_polygons12N_12102018")

end_time <- Sys.time()

end_time - start_time