Processing LANDSAT 8 data in R (using packages Raster, RSToolbox, etc.) in order to calculate NDVI and land cover classification? [closed]

Question

I am analyzing several LANDSAT scenes from the LANDSAT 8 OLI mission (18 in total, all from Jan 2017 - cloud cover < 20%).

I now want to calculate the NDVI and do some land cover classification on these scenes using R and its RS packages, but first need to apply some corrections to these scenes.

If I'm correct, I need to (with examples of codes provided from RStooolbox and raster packages):

1. Do the DN calibration

p198r053_rad<-radCor(p198r053, metaData = metap198r053, method = "rad")

2. Apply an atmospheric correction

p198r053_haze<-estimateHaze(p198r053, darkProp = 0.01, hazeBand= c("B1_dn", "B2_dn", "B3_dn", "B4_dn"))
p198r053_sdos<-radCor(p198r053, metaData = metap198r053, hazeValues= p198r053_haze, hazeBand= c("B1_dn", "B2_dn", "B3_dn", "B4_dn"), method="sdos")

3. Apply a topographic correction (C-correction)

using a DEM and RStoolbox::topCor

4. Mask the clouds using the QA layer

using raster::mask()

5. Check the geolocation of each scene

I don't know how to do so yet (see questions below).

6. Calculate the NDVI for each scene

using

ndvi_formula<-function(nir,red) {
  (nir - red)/(nir + red)
}  

(i.e. bands 5 and 4 of L8 data - but from which layer?? (see question below))

7. Classify images according to land cover type

(perhaps this should be done with the Semi-Automatic Classification tool in QGIS?)

8. Check if all scenes have the same projection and reproject if necessary

Using raster::projection() and raster::projectRaster()

9. Create a mosaic of all scenes for NDVI and land cover layers

Using raster::merge()

10. Crop mosaïcs to match the country's borders

Using a shapefile of the country and raster::mask()

(I have based my approach on the one given by Wegmann et al. (2016) in "Remote Sensing and GIS for Ecologists using Open Source Software")

Are these steps correct and in the right order?

Also, more specifically:

• Steps 1-4 seem to create different rasters for each scene, how do I merge these rasters together to have a scene with all corrections from which to calculate NDVI and do the classification?
• How do I correct the geolocation of my scenes?

Answer 1

If you download surface reflectance data from espa you can skip steps 1, 2, 3, 5 and 6, and you will benefit from a more advanced atmospheric correction.

Mask the clouds using the QA layer

maskFun <- function(x, y){
    x[y %in% c(322, 386)] <- NA
    return(x)
}

# r is the raster* object you want to mask and mask a RasterLayer with
# mask values (e.g. pixel_qa layer in espa products)
overlay(r, mask)

Adjust the masking values to your need (see Landsat 8 SR product guide)

Classify images according to land cover type

Do you have training data? If yes it's a classification and you can train the classifier of your choice (randomForests, glm, randomGLM, ...), if not it's a clustering, and you also have many clustering options (kmeans, dbscan, meanshift). Here's a quick tutorial on how to perform a classification using R.

Check if all scenes have the same projection and reproject if necessary and Create a mosaic of all scenes for NDVI and land cover layers

There's no easy way to do that with raster* classes (merge and mosaic only work on rasters with same projection); your best shot if you want to do the two in one step is to use gdalwarp, which can be called directly from R using the gdalUtils package.

A good reference for the modelling part is the dismo package vignette