changed msscvm2 to msscvm
changed msscvm2 to msscvm

jdbcode/LandsatLinkr,DougFirErickson/LandsatLinkr
#' Create a cloud and cloud shadow mask for Landsat MSS imagery
#' 
#' Takes in any numeric value and squares it.
#' @param file character. MSS reflectance image filename (full system path to MSS file) 
#' @param demfile character. DEM filename (full system path to spatially coincident DEM file)
#' @param topoprep logical. TRUE if slope and aspect are already created in the "topo" folder and FALSE if not
#' @param test logical. If TRUE clouds, cloud shadows and clear pixels have unique values, if FALSE obscured are 0 and clear are 1
#' @return A binary raster with the same dimensions as the MSS image where pixels with value 1 represent clear pixel and 0 as obsured by either cloud or cloud shadow
#' @import raster
#' @import rgdal
#' @import SDMTools
#' @import igraph
#' @export


msscvm = function(file, demfile, topoprep, test=F, overwrite=F){
  #msscvm
  #get the metadata
#   file="E:/llr_test/mixed/mss/wrs1/036032/images/1973/LM10360321973191_reflectance.tif"
#   demfile="E:/llr_test/mixed/mss/wrs1/036032/topo/wrs1_036032_60m_dem.tif"
#   topoprep = T
#   test=F
  
  check = file_check(file,"cloudmask.tif",overwrite)
  print(check)
  if(check == 0){return(0)}
  
  ref = raster(file)
  info = get_metadata(file)
  
  b1 = as.matrix(raster(file, 1))
  b2 = as.matrix(raster(file, 2))
  b4 = as.matrix(raster(file, 4))
  
  #apply topographic correction to band 4 for identifying cloud shadows
  k=0.55
  sunzen = info$sunzen*(pi/180)
  
  #crop the hillshade layer and convert to illumination
  if(topoprep == T){
    dname = dirname(file)
    scenedir = substr(dname,1,nchar(dname)-12)
    topodir = file.path(scenedir,"topo")
    slopefile = list.files(topodir,"slope.tif$",full.name=T)
    aspectfile = list.files(topodir,"aspect.tif$",full.name=T)
    slope = raster(slopefile)
    aspect = raster(aspectfile)
    
    slope_ex  = alignExtent(slope, ref, snap="near")
    extent(slope) = slope_ex
    slope = crop(slope,ref)
    
    aspect_ex  = alignExtent(aspect, ref, snap="near")
    extent(aspect) = slope_ex
    aspect = crop(aspect,ref)
    
    ill = as.matrix(hillShade(slope, aspect, angle=info$sunelev, direction=info$sunaz, normalize=F))
#     hill = hillShade(slope, aspect, angle=info$sunelev, direction=info$sunaz, normalize=F)
#     hillfile = sub("reflectance","hillshade",file)
#     writeRaster(hill, hillfile)
#     ill = as.matrix(hill)
  }
  if(topoprep == F){
    dem = raster(demfile)
    dem_ex  = alignExtent(dem, ref, snap="near")
    extent(dem) = dem_ex
    dem = crop(dem,ref)
    
    #dem = crop(raster(demfile), raster(file))
    slope = terrain(dem, opt="slope")
    aspect = terrain(dem, opt="aspect")
    ill = as.matrix(hillShade(slope, aspect, angle=info$sunelev, direction=info$sunaz, normalize=F))
  }
  
  
  #apply the correction
  c = (cos(sunzen)/ill)^k
  b4topoc = round(b4*c)
  
  #find clouds
  cloud = (b1 - b2)/(b1 + b2)
  clouds = which(cloud > 0.0 & b1 > 1750 | b1 > 3900)
  
  #find the shadows  
  b4topoc[clouds] = NA
  
  b4nocldmean = mean(b4topoc, na.rm=T)
  shadowthresh1 = round(0.40*b4nocldmean + 247.97) #rad = 536.94
  nocldorshdw = which(b4topoc > shadowthresh1) #find non-shadow pixels first pass
  b4nocldorshdw = b4topoc[nocldorshdw]
  b4nocldorshdwmean = mean(b4nocldorshdw, na.rm=T)
  shadowthresh2 = round(0.47*b4nocldorshdwmean + 73.23) #refl = 73.23 rad = 158.56
  shadows = which(b4topoc <= shadowthresh2)
  
  #find the water
  slope = as.matrix(slope)  
  ndvi = (b4-b2)/(b4+b2)
  waterindex = which(ndvi < 0.0850 & slope < (0.5*(pi/180))) #| ndvi < -0.05 & & slope < (1.5*(pi/180))
  
  #create a set of blank rasters for filtering
  b1[] = 0
  water=shadow=cloud=b1 
  b1=0
  
  water[waterindex] = 1
  clumps = .Call("ccl", water, PACKAGE = "SDMTools")  
  clumps = setValues(ref, clumps)
  fre = freq(clumps)
  these = which(fre[,2] < 7)
  values = fre[these,1]
  m = match(as.matrix(clumps), values) 
  these = which(is.na(m) == F)
  water[these] = 0
  
  water = setValues(ref,water)
  water = focal(water, w=matrix(1,5,5), fun=max, na.rm=T)  
  waterindex = which(as.matrix(water) == 1)
  
  shadow[shadows] = 1 
  shadow[waterindex] = 0
  cloud[clouds] = 1
  
  #filter the aggregated cloud and shadow 
  clumps = .Call("ccl", cloud, PACKAGE = "SDMTools")  
  clumps = setValues(ref, clumps)
  fre = freq(clumps)
  these = which(fre[,2] < 10) #10
  values = fre[these,1]
  m = match(as.matrix(clumps), values) 
  these = which(is.na(m) == F)
  cloud[these] = 0
  
  cloud = setValues(ref,cloud)
  cloud = focal(cloud, w=matrix(1,5,5), fun=max, na.rm=F, pad=T, padValue=0)
  
  reso = xres(ref)
  if(reso == 60){
    r = raster(ncol=31,nrow=31)
    ext = extent(0, 31, 0, 31)
    extent(r) = ext
    projection(r) = set_projection(file)
    r[] = NA
    r[16,16] = 1
    dist = gridDistance(r, 1)
    kernal = dist <= 16
  }
  if(reso == 30){
    r = raster(ncol=61,nrow=61)
    ext = extent(0, 61, 0, 61)
    extent(r) = ext
    projection(r) = set_projection(file)
    r[] = NA
    r[31,31] = 1
    dist = gridDistance(r, 1)
    kernal = dist <= 31
  }
  
  shadowproj = focal(cloud, w=as.matrix(kernal), fun=max, na.rm=F, pad=T, padValue=0)
  shiftstart = 1000/tan((pi/180)*info$sunelev)
  shiftend = 7000/tan((pi/180)*info$sunelev)
  
  shiftlen = seq(shiftstart,shiftend,900)
  
  shiftit = function(shadowproj,shiftlen,info,reso){
    if(info$sunaz > 90 & info$sunaz < 180){
      angle = info$sunaz-90  
      yshift = round((sin((pi/180)*angle) * shiftlen)/reso)*reso
      xshift = round((cos((pi/180)*angle) * shiftlen * -1)/reso)*reso
      shadowproj = shift(shadowproj, x=xshift, y=yshift)
    }
    if(info$sunaz > 0 & info$sunaz < 90){
      angle = 90-info$sunaz
      angle = 10
      yshift = round((sin((pi/180)*angle) * shiftlen *-1)/reso)*reso
      xshift = round((cos((pi/180)*angle) * shiftlen *-1)/reso)*reso
      shadowproj = shift(shadowproj, x=xshift, y=yshift)
    }
    return(shadowproj)
  }
  
  for(m in 1:length(shiftlen)){
    if(m == 1){mergeit = "r1"} else {mergeit = paste(mergeit,",r",m, sep="")}
    dothis = paste("r",m,"=shiftit(shadowproj,shiftlen[m],info,reso)", sep="")
    eval(parse(text=dothis))
    if(m == length(shiftlen)){mergeit = paste("shadowproj = mosaic(",mergeit,",fun=max)", sep="")}
  }
  
  #run the merge function
  if(length(shiftlen) == 1){shadowproj = r1} else {eval(parse(text=mergeit))}
  
  #ake sure that all values are finite, the mosaicing function with max can cause some problems where there are no pixels
  shadowproj[!is.finite(values(shadowproj))] = 0
  
  #extend the layer so it has a full union with the cloud layer
  shadowproj = extend(shadowproj, cloud, value=0)
  
  #crop the cloud projection layer by the cloud layer
  shadowproj = crop(shadowproj, cloud)
  
  #convert the shadow layer to a matrix
  shadow = setValues(ref,shadow)
  
  #get the intersection of shadow and cloud projection
  shadow = shadow*shadowproj
  
  #convert the shadow layer to a matrix for spatial sieve
  shadow = as.matrix(shadow)
  
  #filter the aggregated cloud and shadow 
  clumps = .Call("ccl", shadow, PACKAGE = "SDMTools")  
  clumps = setValues(ref, clumps)
  fre = freq(clumps)
  these = which(fre[,2] < 10) #10
  values = fre[these,1]
  m = match(as.matrix(clumps), values) 
  these = which(is.na(m) == F)
  shadow[these] = 0
  shadow = setValues(ref, shadow)
  
  shadow = focal(shadow, w=matrix(1,5,5), fun=max, na.rm=F, pad=T, padValue=0)
  if(test == T){cloud = cloud*2
                cloudshadow = mosaic(cloud,shadow,fun=max, na.rm=T)
  } else {cloudshadow = sum(cloud, shadow, na.rm=T)}
  
  
  if(test == F){cloudshadow = setValues(ref,as.numeric(values(cloudshadow) == 0))}
  cloudshadow[is.na(ref)] = NA
  projection(cloudshadow) = set_projection(file)
  cloudshadow = as(cloudshadow, "SpatialGridDataFrame")
  if(test == F){outfile = sub("reflectance", "cloudmask", file)} else {outfile = sub("reflectance", "cloudmask_test", file)}
  writeGDAL(cloudshadow, outfile, drivername = "GTiff", type = "Byte", mvFlag=255, options="INTERLEAVE=BAND")
  
}
#' Create a cloud and cloud shadow mask for Landsat MSS imagery
#' 
#' Takes in any numeric value and squares it.
#' @param file MSS image filename (full system path to MSS file) 
#' @param demfile DEM filename (full system path to spatially coincident DEM file)
#' @return A binary raster with the same dimensions as the MSS image where pixels with value 1 represent clear pixel and 0 as obsured by either cloud or cloud shadow
#' @import raster
#' @import SDMTools
#' @export


msscvm = function(file, demfile){
  
#    library(SDMTools)
#    library(raster)
#    library(MSScvm)
# #   
#    file = "E:/mss/wrs1/036032/images/1975/LM20360321975226_radiance.tif"
#    demfile = "K:/gis_data/dems/wrs_dem/wrs1_036032_60m_dem.tif"
  
  #get the metadata
  info = get_metadata(file)
  
  b1 = as.matrix(raster(file, 1))
  b2 = as.matrix(raster(file, 2))
  b4 = as.matrix(raster(file, 4))
  
  #apply topographic correction to band 4 for identifying cloud shadows
  k=0.55
  sunzen = info$sunzen*(pi/180)
  
  #crop the hillshade layer and convert to illumination
  dem = raster(demfile)
  ref = raster(file)
  dem_ex  = alignExtent(dem, ref, snap="near")
  extent(dem) = dem_ex
  dem = crop(dem,ref)

#dem = crop(raster(demfile), raster(file))
  slope = terrain(dem, opt="slope")
  aspect = terrain(dem, opt="aspect")
  ill = as.matrix(hillShade(slope, aspect, angle=info$sunelev, direction=info$sunaz, normalize=F))
  
  #apply the correction
  c = (cos(sunzen)/ill)^k
  b4topoc = round(b4*c)
  
  #find clouds
  cloud = (b1 - b2)/(b1 + b2)
  cloudtight = which(cloud > 0.0450 & b1 > 7500 | b1 > 14000)
  
  #find the shadows  
  b4topoc[cloudtight] = NA
  
  b4nocldmean = mean(b4topoc, na.rm=T)
  shadowthresh1 = round(0.40*b4nocldmean + 536.94)
  nocldorshdw = which(b4topoc > shadowthresh1) #find non-shadow pixels firest pass
  b4nocldorshdw = b4topoc[nocldorshdw]
  b4nocldorshdwmean = mean(b4nocldorshdw, na.rm=T)
  shadowthresh2 = round(0.47*b4nocldorshdwmean + 158.56)
  shadows = which(b4topoc <= shadowthresh2)
  
  #find the water
  slope = as.matrix(slope)  
  ndvi = (b4-b2)/(b4+b2)
  waterindex = which(ndvi < -0.2 & slope < (0*(pi/180)))    #slope = 0.0 degrees
  shoreindex = which(ndvi < -0.3 & slope < (0.5*(pi/180)))  #slope = 0.5 degrees
  
  b2[] = 0
  b2[c(waterindex,shoreindex)] = 1
  clumps = .Call("ccl", b2, PACKAGE = "SDMTools")  
  clumps = setValues(raster(file), clumps)
  #clumps = ConnCompLabel(b2)
  
  fre = freq(clumps)
  these = which(fre[,2] < 4)
  values = fre[these,1]
  m = match(as.matrix(clumps), values) 
  these = which(is.na(m) == F)
  #b2 = as.matrix(b2)
  b2[these] = 0
  b2 = setValues(raster(file),b2)
  
  b2 = focal(b2, w=matrix(1,5,5), fun=max) 
  b2 = as.matrix(b2)
  
  #create a blank raster and fill it in
  b1[] = 0
  b1[cloudtight] = 1  
  b1[shadows] = 1
  b1[b2 == 1] = 0 #set water back to clear
  #b1 = setValues(raster(file), b1)
  clumps = .Call("ccl", b1, PACKAGE = "SDMTools")  
  clumps = setValues(raster(file), clumps)
  #clumps = ConnCompLabel(b1)  
  
  fre = freq(clumps)
  these = which(fre[,2] < 10)
  values = fre[these,1]
  m = match(as.matrix(clumps), values) 
  these = which(is.na(m) == F)
  #b1 = as.matrix(b1)
  b1[these] = 0
  b1 = setValues(raster(file),b1)
  
  b1 = focal(b1, w=matrix(1,5,5), fun=max)
  b1[] = as.numeric(values(b1) == 0) 
  b1[is.na(raster(file))] = NA
  projection(b1) = set_projection(file)
  b1 = as(b1, "SpatialGridDataFrame")
  outfile = sub("radiance", "cloudmask", file)
  writeGDAL(b1, outfile, drivername = "GTiff", type = "Byte", mvFlag=255, options="INTERLEAVE=BAND")
}





