R/SENTINEL-2Processing.R

In dwalin93/SENTINEL2Processing:

test = function(link,image){
  shape = rgdal::readOGR(dsn = link, layer = basename(strsplit(link, "\\.")[[1]])[1])
  raster = raster::raster(image)
  raster::plot(image)
  raster::plot(shape,add=T)
  }

#'Returns the image with differences as Polygons on top
#'@param NDVI1 NDVI Result Image 1
#'@param NDVI2 NDVI Result Image 2
showDifferencesOnImage = function(NDVI1,NDVI2,TCI1,TCI2){
  greyscale = gray.colors(100,
                          start = 0.0,
                          end = 1.0,
                          gamma=2.2,
                          alpha = NULL)
  mycol = heat.colors(5,alpha=0.3)
  polygons = NDVIDifferences(NDVI1,NDVI2)
  TCI1 = raster::raster(TCI1)
  TCI2 = raster::raster(TCI2)
  #polygons[polyRaster<=0] = NA
  graphics::par(bg=NA,mar=c(0,0,0,0),oma=c(0,0,0,0))
  raster::image(TCI1,axes=FALSE, frame=FALSE, col=greyscale)
  raster::image(polygons,axes=FALSE, frame=FALSE,col=mycol,add=T)
  raster::image(TCI2,axes=FALSE, frame=FALSE,col=greyscale)
  raster::image(polygons,axes=FALSE, frame=FALSE,col=mycol,add=T)
}


#' Returns the Polygons used to show differences
#' @param NDVI1 NDVi Result Image 1
#' @param NDVI2 NDVI Result Image 2
NDVIDifferences = function(NDVI1,NDVI2){
  normalizeFirst = rescale(raster::raster(NDVI1))
  normalizeSecond = rescale(raster::raster(NDVI2))
  result = normalizeFirst-normalizeSecond
  values = result[result <= -0.01:0.01] = NA
  return (result)
}

#' Rescales the NDVI image from -1 to 1
#' @param x The NDVI image
#' @param x.min min value of image
#' @param x.max max value of image
#' @param new.min -1
#' @param new.max 1
rescale = function(x, x.min = 0, x.max = 255, new.min = -1, new.max = 1) {
  if(is.null(x.min)) x.min = min(x)
  if(is.null(x.max)) x.max = max(x)
  new.min + (x - x.min) * ((new.max - new.min) / (x.max - x.min))
}

#' Checks if a shapefile overlaps the image
#' @param image image to proof
#' @param shape shape to be proofed
checkExtent = function(image,shape){
  raster = raster::raster(image)
  shape = readShape(shape)
  ref = sp::proj4string(shape)
  raster::crs(raster)=ref
  print(raster)
  return (is.null(raster::intersect(raster::extent(raster),shape)))
}


#' Calculates the NDVI only in the Shapefile areas
#' @param NIR NIR Band of Image
#' @param Red Red Band of image
#' @param shapeLink Link to shapefile
calcNDVI = function(NIR,Red,shapeLink){
  if (checkExtent(NIR,shapeLink) == FALSE){
    print('HEHO')
    NIR = cropRaster(NIR,shapeLink)
    Red = cropRaster(Red,shapeLink)
    NDVI = NDVI_Result(NIR,Red)
    graphics::par(bg=NA,mar=c(0,0,0,0),oma=c(0,0,0,0))
    raster::image(NDVI,axes=FALSE, frame=FALSE, col=base::rev(terrain.colors(3)))
  } else {
    NIR = raster::raster(NIR)
    Red = raster::raster(Red)
    NDVI = NDVI_Result(NIR,Red)
    graphics::par(bg=NA,mar=c(0,0,0,0),oma=c(0,0,0,0))
    raster::image(NDVI,axes=FALSE, frame=FALSE, col=base::rev(terrain.colors(3)))
  }
}


#'Returns the raster based on the shapefiles
#'@param raster Raster to be croped
#'@param shapeLink link to shapefile
cropRaster = function(raster,shapeLink){
  link = shapeLink
  shape = readShape(link)
  plot(shape)
  band = raster::raster(raster)
  ref = sp::proj4string(shape)
  raster::crs(band)=ref
  crop = crop(band,raster::extent(shape), snap="out")
  frame = raster::rasterize(shape,crop)
  layer = raster::mask(x=crop,mask=frame)
  raster::plot(shape)
  graphics::par(bg=NA,mar=c(0,0,0,0),oma=c(0,0,0,0))
  return(layer)
}


#'Reads in the Shapefile
#'@param link Link to the shapefile
readShape = function(link){
  shape = rgdal::readOGR(dsn = link, layer = basename(strsplit(link, "\\.")[[1]])[1])

}

#' Returns the NDVI (NIR-Red/NIR+Red)
#' @param x NIR Band
#' @param y Red Band
NDVI = function(x,y){
  result = (x-y)/(x+y)
  return (result)
}


#' Returns the NDVI image
#' @param x NIR Band
#' @param y Red band
#' @return return NDVI image
NDVI_Result = function(x,y){
  output = raster::overlay(x,y, fun = NDVI)
  return (output)

}

#'Plots a false color composite image
#'@param R Red channel
#'@param G Green Channel
#'@param B Blue Channel
FCC = function (R,G,B){
  rasterR = raster::raster(R)
  rasterG = raster::raster(G)
  rasterB = raster::raster(B)
  RGBBrick = raster::brick(rasterR,rasterG,rasterB)
  graphics::par(bg=NA,mar=c(0,0,0,0),oma=c(0,0,0,0))
  raster::plotRGB(RGBBrick,r=1,g=2,b=3,axes=FALSE,xaxs="i", yaxs="i")

}