1st upload of completed package

Author: unknown

.Rbuildignore

@@ -0,0 +1,2 @@
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -0,0 +1,3 @@
+.Rproj.user
+.Rhistory
+.RData

DESCRIPTION

@@ -0,0 +1,17 @@
+Package: MSScvm
+Title: Automated Cloud Masking for Landsat MSS Images 
+Version: 1.0.0
+Date: 2015-07-23
+Authors@R: c(person("Justin", "Braaten", email = "jstnbraaten@gmail.com", role = c("aut", "cre")),
+    person("Warren", "Cohen", email = "warren.cohen@oregonstate.edu", role = "aut"),
+    person("Zhiqiang", "Yang", email = "zhiqiang.yang@oregonstate.edu", role = "aut"))
+Description: An automated cloud and cloud shadow masking system for Landsat MSS imagery. It allows MSS imagery to more easily be used in large-area and time series analysis by providing an efficent way to prevent cloud and cloud shadow pixels from contaminating mosaics, composites, and time series.
+Depends: R (>= 3.2.1)
+Imports:
+  gdalUtils,
+  igraph,
+  raster,
+  rgdal,
+  SDMTools
+License: GPL-2
+LazyData: true

LICENSE

@@ -0,0 +1,18 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source
+PackageRoxygenize: rd,collate,namespace,vignette

MSScvm.Rproj

@@ -0,0 +1,10 @@
+# Generated by roxygen2 (4.1.1): do not edit by hand
+
+export(MSScvm)
+export(MSSdn2rad)
+export(MSSdn2refl)
+export(MSSunpack)
+export(eudist)
+export(getMetadata)
+export(mosaicDEMs)
+export(reprojectDEM)

NAMESPACE

@@ -0,0 +1,34 @@
+#' Convert MSS DN values to TOA radiance
+#'
+#' Convert MSS DN values to TOA radiance
+#' @param imgFile filename (character). Full path to *dn.tif image file produced by the \code{\link{MSSunpack}} function.
+#' @details The equation used to convert DN to TOA radiance can be found \href{http://landsat.usgs.gov/how_is_radiance_calculated.php}{here}.
+#' @return A 4-band Landsat MSS GeoTIFF raster image file in units of top-of-atmosphere (TOA) radiance. The file will be placed in
+#' same directory as the 'imgFile' with the name equal to the image ID followed by 'toa_radiance'. Note that the values are scaled by 100
+#' and rounded to the nearest integer to reduce the file size.
+#' @seealso \code{\link{MSSdn2refl}}
+#' @examples \dontrun{
+#' 
+#' MSSdn2rad("C:/mss/LM10360321973191AAA04/LM10360321973191AAA04_dn.tif")
+#' }
+#' @export
+
+MSSdn2rad = function(imgFile){
+  
+  info = getMetadata(imgFile) #get image metadata
+  b = raster::brick(imgFile) #load the DN image as a brick
+  img = as.array(b)  #convert the brick to an array
+
+  img[,,1] = round(100*((info$b1gain*img[,,1])+info$b1bias)) #convert fro DN to radiance and scale by 100 and round
+  img[,,2] = round(100*((info$b2gain*img[,,2])+info$b2bias)) #convert fro DN to radiance and scale by 100 and round
+  img[,,3] = round(100*((info$b3gain*img[,,3])+info$b3bias)) #convert fro DN to radiance and scale by 100 and round
+  img[,,4] = round(100*((info$b4gain*img[,,4])+info$b4bias)) #convert fro DN to radiance and scale by 100 and round
+  
+  img[img < 0] = 0 #you can't have negative radiance - set negative values to 0
+  
+  img = raster::setValues(b,img) #convert the array to a brick 
+  img = as(img, "SpatialGridDataFrame") #convert the brick to SGHF to be written by GDAL
+  outfile = sub("dn", "toa_radiance", imgFile) #define output filename
+  rgdal::writeGDAL(img, outfile, drivername = "GTiff", type = "Int16", mvFlag = -32768, options="INTERLEAVE=BAND") #write the file
+}
+  

R/MSScvm.r

@@ -0,0 +1,58 @@
+#' Convert DN values to TOA reflectance
+#'
+#' Convert DN values to TOA reflectance
+#' @param imgFile filename (character). Full path to *dn.tif image file produced by the \code{\link{MSSunpack}} function.  
+#' @details DN values are first converted to top-of-atmosphere (TOA) radiance using the equation found \href{http://landsat.usgs.gov/how_is_radiance_calculated.php}{here}.
+#' Then TOA radiance is converted to TOA reflectance using the equation found \href{http://landsathandbook.gsfc.nasa.gov/data_prod/prog_sect11_3.html}{here}.
+#' The ESUN values used are from the publication 'Chander et al. 2009. Summary of current radiometric calibration coefficients... Remote Sensing of Environment. 113'. 
+#' @return A 4-band Landsat MSS GeoTIFF raster image file in units of top-of-atmosphere (TOA) reflectance. The file will be placed in the
+#' same directory as the 'imgFile' with the name equal to the image ID followed by 'toa_reflectance'. Note that the values are scaled by 10,000
+#' and rounded to the nearest integer to reduce the file size.
+#' @seealso \code{\link{MSSdn2rad}}
+#' @examples \dontrun{
+#' 
+#' MSSdn2refl("C:/mss/LM10360321973191AAA04/LM10360321973191AAA04_dn.tif")
+#' }
+#' @export
+
+MSSdn2refl = function(imgFile){
+  
+  #link to the equations to convert DN to TOA and TOA to SR
+  #http://landsathandbook.gsfc.nasa.gov/data_prod/prog_sect11_3.html
+  
+  #define the TOA reflectance function
+  refl = function(imgFile, band, gain, bias, sunzen, d, esun){
+    orig = raster::raster(imgFile, band) #read in the file
+    img = as.matrix(orig) #convert from raster to matrix
+    img = ((gain*img)+bias) #convert from DN to TOA radiance
+    img[img < 0] = 0 #set all values less than 0 to 0 - negative radiance does not make sense
+    img = (pi * img * (d^2))/(esun * cos(sunzen)) #convert from TOA radiance to TOA reflectance
+    img = round(img * 10000) #scale by 10,000 and round to save on file size
+    img = raster::setValues(orig,img) #convert from matrix to raster
+    return(img) #return the TOA radiance raster
+  }
+  
+  info = getMetadata(imgFile) #read in the image metadata
+  
+  #define esun values for mss (chander et al 2009 summary of current radiometric calibration coefficients... RSE 113)
+  if(info$sensor == "LANDSAT_1"){esun = c(1823,1559,1276,880.1)}
+  if(info$sensor == "LANDSAT_2"){esun = c(1829,1539,1268,886.6)}
+  if(info$sensor == "LANDSAT_3"){esun = c(1839,1555,1291,887.9)}
+  if(info$sensor == "LANDSAT_4"){esun = c(1827,1569,1260,866.4)}
+  if(info$sensor == "LANDSAT_5"){esun = c(1824,1570,1249,853.4)}
+  
+  d = eudist(info$doy) #define the earth sun distance
+  sunzen = info$sunzen*(pi/180) #get sun zenith angle in radians for cos 
+  
+  #apply the conversion to reflectance using the 'refl' function
+  b1 = refl(imgFile,1,info$b1gain, info$b1bias, sunzen, d, esun[1]) 
+  b2 = refl(imgFile,2,info$b2gain, info$b2bias, sunzen, d, esun[2])
+  b3 = refl(imgFile,3,info$b3gain, info$b3bias, sunzen, d, esun[3])
+  b4 = refl(imgFile,4,info$b4gain, info$b4bias, sunzen, d, esun[4])
+  
+  img = raster::stack(b1,b2,b3,b4) #stack the bands
+  
+  img = as(img, "SpatialGridDataFrame")   #convert the raster to SGHF so it can be written using GDAL (faster than writing it with the raster package)
+  outfile = sub("dn", "toa_reflectance", imgFile) #define the outputs file
+  rgdal::writeGDAL(img, outfile, drivername = "GTiff", type = "Int16", mvFlag = -32768, options="INTERLEAVE=BAND") #write out the TOA reflectance file
+}

R/MSSdn2rad.r

@@ -0,0 +1,84 @@
+#' Decompress and stack Landsat LPGS MSS images
+#'
+#' Decompresses and stacks Landsat LPGS MSS images provided by USGS as *.tar.gz files. Optionally 
+#' outputs top-of-atmosphere (TOA) radiance and reflectance files.
+#' @param imgFile filename (character). Full path to compressed LPGS Landsat MSS image file provided by USGS.  
+#' @param toaRad logical. If TRUE, a TOA radiance image will be created.
+#' @param toaRefl logical. If TRUE, a TOA reflectance image will be created. 
+#' @param useL1G logical. If TRUE, L1G images will be processed.
+#' @details It is important that the 'imgFile' be an unaltered tar.gz-compressed LPGS image file that you receive 
+#' from USGS through \href{http://rstudio.com}{EarthExplorer}. Note that DN values <= 1 are set to NA across all bands. 
+#' This mitigates a problem caused by bad columns on the east and west edge of images when mosaicing adjacent images together.  
+#' @return A 4-band Landsat MSS GeoTIFF raster image file in DN units. If optional 'toaRad' and/or 'toaRefl' 
+#' parameters are set to TRUE, then similar TOA radiance and reflectance image files will created. The files will be places in
+#' directory in the same location as the 'imgFile' with the name equal to the image ID. The files contain the image ID with
+#' descriptors 'dn' (digital number), 'toa_radiance' (TOA radiance), and 'toa_reflectance' (TOA reflectance).
+#' @seealso \code{\link{MSSdn2rad}}, \code{\link{MSSdn2refl}}
+#' @examples 
+#' \dontrun{
+#' 
+#' MSSunpack(imgFile = "C:/mss/LM10360321973191AAA04.tar.gz")
+#' MSSunpack(imgFile = "C:/mss/LM10360321973191AAA04.tar.gz", 
+#'           toaRad = FALSE, toaRefl = TRUE, useL1G = TRUE)
+#' }
+#' @export
+
+MSSunpack = function(imgFile, toaRad=FALSE, toaRefl=FALSE, useL1G=FALSE){
+  
+  #figure out if this is an L1G image and stop the function according to the 'useL1G' parameter
+  tempdir = file.path(dirname(imgFile),"temp")
+  untar(imgFile, exdir=tempdir)
+  if(useL1G == F){
+    mtlfile = list.files(tempdir, pattern = "MTL.txt", full.names = T, recursive = T) #find the mtl metadata file
+    info = getMetadata(mtlfile)
+    if(info$datatype == "L1G"){
+      print(paste("MSS file:",imgFile,"is L1G, if you want to use it, set the 'useL1G' parameter to TRUE"))
+      stop("Stopping MSSunpack")
+    }
+  }
+  
+  #get file and directory info
+  allfiles = list.files(tempdir, full.names=T) #find all the decompressed files 
+  tiffiles = allfiles[grep("TIF",allfiles)] #subset the tif image files
+  ancfiles = allfiles[grep("TIF",allfiles, invert=T)] #subset the other files
+  filebase = basename(tiffiles[1]) #get the basename
+  filedir = dirname(imgFile) #get the directory
+  
+  #create stack output name
+  name = paste(info$imgid, "_dn.tif", sep = "") #define the new file basename 
+  outdir = file.path(filedir, info$imgid) #define the output directory
+  finalstack = file.path(outdir, name) #define the new full filename of the output image 
+  
+  #deal with the ancillary file names
+  baseancfiles = basename(ancfiles) #get the basenames of the ancillary files
+  newancfiles =  file.path(outdir, baseancfiles) #define the new filenames for ancillary files
+  
+  #start working with the image files - stack and set bad pixels to 0
+  s = stack(tiffiles[1],tiffiles[2],tiffiles[3],tiffiles[4]) #stack the band files
+  img = as.array(s) #convert to array for fastering processing
+  b1bads = img[,,1]>1 #finds bad image edge pixels that cause problems when mosaicing 
+  b2bads = img[,,2]>1 #finds bad image edge pixels that cause problems when mosaicing
+  b3bads = img[,,3]>1 #finds bad image edge pixels that cause problems when mosaicing
+  b4bads = img[,,4]>1 #finds bad image edge pixels that cause problems when mosaicing
+  bads = b1bads*b2bads*b3bads*b4bads #combines all of the bad pixels
+  
+  img[,,1] = img[,,1]*bads #sets all the bad pixels to value 0
+  img[,,2] = img[,,2]*bads #sets all the bad pixels to value 0
+  img[,,3] = img[,,3]*bads #sets all the bad pixels to value 0
+  img[,,4] = img[,,4]*bads #sets all the bad pixels to value 0
+  
+  dir.create(outdir, recursive=T, showWarnings=F) #make a new output directory
+  
+  file.rename(ancfiles,newancfiles) #move the associated files
+  
+  #write out the dn file
+  outimg = raster::setValues(s,img) #convert from matrix to raster   
+  outimg = as(outimg, "SpatialGridDataFrame") #convert from rater to SGDF for faster writing
+  rgdal::writeGDAL(outimg, finalstack, drivername = "GTiff", options="INTERLEAVE=BAND", type = "Byte", mvFlag = 0) #write out the image.
+  
+  if(toaRad == T){MSSdn2rad(finalstack)} #calculate and write the toa radiance file if flagged
+  
+  if(toaRefl == T){MSSdn2refl(finalstack)} #calculate and write the toa reflectance file if flagged
+  
+  unlink(tempdir, recursive=T, force=T) #delete the temp folder and its contents
+}