Unsupervided Classification: LULC

# Unsupervided Classification of Landsat images
# !UTIL http://bioinf.wehi.edu.au/marray/ibc2004/lab4/lab4.html
# !UTIL http://www.digital-geography.com/unsupervised-classification-of-a-landsat-image-in-r-the-whole-story-or-part-two/
# !UTIL http://www.jstatsoft.org/v15/i09/paper
# !UTIL http://cran.r-project.org/web/packages/e1071/e1071.pdf

# Packages
kpacks <- c('raster', 'sp', 'rgdal', 'cluster')
new.packs <- kpacks[!(kpacks %in% installed.packages()[,"Package"])]
if(length(new.packs)) install.packages(new.packs)
lapply(kpacks, require, character.only=T)

# Folders
dir.tif <- '..\\Landsatfolder' # TIF file folder
dir.rst   <-'..\\Landsatfolder' # Destination folder

# Extent Kumbira study area
# (xmin, xmax, ymin, ymax):
#   416985.0000000,   425295.0000000
# -1241505.0000000, -1226985.0000000
kae <- c(416985, 425295, -1241505, -1226985)

# Capture Projection String
proj.i <- CRS("+init=epsg:32633") # UTM 33N

# Read Landsat Geotiff
# Crop image files -----------------------------------

all.files <- list.files(dir.tif, all.files = F)
# List of TIF files at dir.tif folder
list.tif <- grep(".tif$", all.files, ignore.case = TRUE, value = TRUE) 
stk <- brick()
for (i in 1:length(list.tif)) {
  # Name
  fname <- substr(list.tif[i], 1, (nchar(list.tif[i]) - 4))
  # Read Geotiff raster
  img.tmp <- raster(paste0(dir.tif, list.tif[i]), package = "raster", varname = fname)
  projection(img.tmp) <- proj.i # projection association
  # Check for a valid projection 
  if(is.null(projection(img.tmp, asText = TRUE))) {
    stop('no projection defined!')
  } 
  # apply Crop
  img.crop.m <- crop(img.tmp, kae)
  # apply Mask

} else {
  img.crop <- crop(img.tmp, ae) # apply windowing only
  stopifnot(dataType(img.crop) == 'INT2U')
}
# Write Original DN bands
if(compareRaster(img.crop, r.ref, extent=F, rowcol=F, crs=F, res=T, orig=F,
                 rotation=F, values=F, tolerance = 0, stopiffalse=F, showwarning=F)) {
  stk <- addLayer(stk, img.crop)   
}
writeRaster(img.crop, filename = paste0(dir.rst, fname,'.rst'),
            format = 'IDRISI', overwrite = TRUE,
            #datatype='INT2U',
            crs = proj.i, NAflag = -9999)
# Write IDRISI raster group rgf for DN bands
fileConn <- file(paste0(dir.rst, "ae_dn.rgf")) # raster group
writeLines(c(length(list.tif), substr(list.tif, 1, (nchar(list.tif[1]) - 4))), fileConn)
close(fileConn)
closeAllConnections()
remove(img.tmp, img.crop)
}

# Resample SRTM 90m to Landsat resolution
# CompareRaster
mask1 <- resample(dem, img.crop.m, method = "ngb")

# Extract DEM variables
# With raster package: Height, Slope, Aspect
#  Topographic Position Index
#  Terrain Ruggedness Index
# With Landsat package: Wetness Index



stk <- f.cropidrisidata(dir.tif = dir.tif, 
                        dir.rst = dir.rst,
                        crop = 'ae')