changes

Author: unknown

.Rbuildignore

@@ -1,2 +1,3 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
+^\.travis\.yml$

.travis.yml

@@ -0,0 +1,14 @@
+# Sample .travis.yml for R projects
+
+language: r
+warnings_are_errors: true
+sudo: required
+
+env:
+ global:
+   - CRAN: http://cran.rstudio.com
+
+notifications:
+  email:
+    on_success: change
+    on_failure: change

DESCRIPTION

@@ -1,17 +1,20 @@
+Type: Package
 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.
+Description: An automated cloud and cloud shadow masking system for Landsat MSS imagery. It provides a means of more easily incorporating MSS imagery in large-area and time series analysis by providing an efficient 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
+SystemRequirements: GDAL binaries
 License: GPL-2
+URL: http://www.msscvm.jdbcode.com/
 LazyData: true

R/MSScvm.r

@@ -16,7 +16,9 @@
 
 
 MSScvm = function(imgDir, demFile, classify=F){
-
+  
+  
+  
   #get the image file and check
   files = list.files(imgDir, full.names=T)
   refl_match = grep("toa_reflectance.tif$", files)
@@ -36,12 +38,22 @@ MSScvm = function(imgDir, demFile, classify=F){
   }
 
   imgfile = files[match]
+  print(paste("Making cloud & shadow mask for",basename(imgfile)))
 
   #get some info about the image
   ref = raster::raster(imgfile) #read in the MSS file for information on image extent and as a template for holding values later
   dem = raster::raster(demFile) #read in the DEM file - raster package function
   demproj = raster::projection(dem)
   imgproj = raster::projection(ref)
+  demres = raster::xres(dem)
+  imgres = raster::xres(ref)
+  
+  #check to make the DEM and image resolutions are the same
+  if(demres != imgres){
+    print(paste("The DEM file:",demFile,"does not have the same pixel resolution as the image file."))
+    print("Please make sure the DEM file has the same pixel resolution as the image file. Use the function 'reprojectDEM' to assist in getting it in the same resolution")
+    stop("Stopping MSScvm")
+  }
 
   #check to make the DEM and image projections are the same
   if(demproj != imgproj){
@@ -99,13 +111,15 @@ MSScvm = function(imgDir, demFile, classify=F){
     b4 = refl(imgfile, 4 ,info$b4gain, info$b4bias, sunzen, d, esun[4]) #get TOA reflectance for MSS band 4
   } 
 
+
   if(imgtype == "refl") {
     #load in the image bands
     b1 = raster::as.matrix(raster::raster(imgfile, 1)) #band 1
     b2 = raster::as.matrix(raster::raster(imgfile, 2)) #band 2
     b4 = raster::as.matrix(raster::raster(imgfile, 4)) #band 4
   }
 
+  
   #crop the hillshade layer
   dem_ex  = raster::alignExtent(dem, ref, snap="near") #aligned the extent of the DEM to the image - raster package function
   raster::extent(dem) = dem_ex #set the DEM extend to aligned extent - raster package function

R/MSSdn2rad.r

@@ -15,9 +15,10 @@
 
 MSSdn2rad = function(imgFile){
 
+  print(paste("Converting",basename(imgFile),"to TOA radiance"))
   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 = raster::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

R/MSSdn2refl.r

@@ -17,13 +17,14 @@
 
 MSSdn2refl = function(imgFile){
 
+  print(paste("Converting",basename(imgFile),"to TOA reflectance"))
   #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 = raster::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

R/MSSunpack.r

@@ -26,15 +26,15 @@
 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
+  print(paste("Unpacking and preparing",basename(imgFile)))
   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"){
+  mtlfile = list.files(tempdir, pattern = "MTL.txt", full.names = T, recursive = T) #find the mtl metadata file
+  info = getMetadata(mtlfile)
+  if(useL1G == F & info$datatype == "L1G"){
+      unlink(tempdir, recursive=T, force=T) #delete the temp folder and its contents
       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
@@ -54,8 +54,8 @@ MSSunpack = function(imgFile, toaRad=FALSE, toaRefl=FALSE, useL1G=FALSE){
   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
+  s = raster::stack(tiffiles[1],tiffiles[2],tiffiles[3],tiffiles[4]) #stack the band files
+  img = raster::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

R/mosaicDEMs.r

@@ -3,10 +3,11 @@
 #' A helper function to create the large-extent DEM file required by the 'MSScvm' function.
 #' @param dir directory name (character). Full path to a directory containing digital elevation model (DEM) files to be mosaiced.
 #' @param projRef filename (character). Full path to an image file produced by the \code{\link{MSSunpack}} function.
-#' @param NAvalue numeric. What is the background value of the DEM files in the directory. 
+#' @param srcNodata numeric. What is the background value of the DEM files in the directory. If there is no background value, use NA (default) 
+#' @param dstNodata numeric. Specify the value to represent background pixels in the mosaic DEM. -32768 is the default
 #' @details The provided directory path should only contain decompressed digital elevation files from the same source (SRTM, NED, GTOPO, etc).
 #' The function will search the directory and include all files found in the mosaic. It is important that each file have the same background value
-#' and that it is correctly assigned to the 'NAvalue' parameter, if not, intersection between DEMs could have unexpected results. 
+#' and that it is correctly assigned to the 'srcNodata' parameter, if not, intersection between DEMs could have unexpected results. 
 #' Each individual DEM file will be adjusted to match the projection and pixel resolution of the 'proRef' image. 
 #' Then they will be merged using the mean value of intersecting pixels.  
 #' @return A GeoTIFF raster file representing the union of all individual DEM files found in the provided directory path. The mosaic file will be written to the 
@@ -16,12 +17,12 @@
 #' 
 #' mosaicDEMs(dir = "C:/mss/dems", 
 #'            projRef = "C:/mss/LM10360321973191AAA04/LM10360321973191AAA04_dn.tif",
-#'            NAvalue = -32768)
+#'            srcNodata = -9999, dstNodata= -32768)
 #' }
 #' @export
 
 
-mosaicDEMs = function(dir, projRef, NAvalue){
+mosaicDEMs = function(dir, projRef, srcNodata=NA, dstNodata=-32768){
 
   align = function(img, refimg){
     img = raster::raster(img)
@@ -35,16 +36,26 @@ mosaicDEMs = function(dir, projRef, NAvalue){
   proj = raster::projection(template)
   demfiles = normalizePath(list.files(dir, full.names=T))
 
+  #check files to see if they are compressed
+  zip = grep("zip", demfiles)
+  tar = grep("tar", demfiles)
+  gz =  grep("gz", demfiles)
+  if(sum(length(zip),length(tar),length(gz)) != 0){
+    print(paste("There are compressed files in this directory:",dir))
+    print("make sure the directory only contains decompressed files")
+    stop("Stopping mosaicDEMs")
+  }
+  
   for(i in 1:length(demfiles)){
     demfile = demfiles[i]
-    print(paste("Reprojecting file:",demfile))
+    print(paste("Reprojecting",basename(demfile)))
     bname = basename(demfile)
     extension = substr(bname,(nchar(bname)-3),nchar(bname))
     dstfile = sub(extension, "_reprojected.tif", demfile) 
 
     gdalUtils::gdalwarp(srcfile=demfile, dstfile=dstfile, 
-      t_srs=proj, of="GTiff",
-      r="near", dstnodata=NAvalue, multi=T,
+      t_srs=proj, of="GTiff", overwrite=TRUE,
+      r="near", srcnodata=srcNodata, dstnodata=dstNodata, multi=T,
       tr=c(reso,reso), co="INTERLEAVE=BAND")
   }
 
@@ -58,21 +69,17 @@ mosaicDEMs = function(dir, projRef, NAvalue){
 
   refimg = raster::raster(demfiles[1])
 
+  outfile = file.path(dir,"dem_mosaic.tif")
   big = NA #create a dummy variable to be overwritten on-the-fly with eval
   for(i in 1:len){
     if(i == 1){mergeit = "r1"} else {mergeit = paste(mergeit,",r",i, sep="")}
     dothis = paste("r",i,"=align(demfiles[",i,"], refimg)", sep="")
     eval(parse(text=dothis))
-    if(i == len){mergeit = paste("big = raster::mosaic(",mergeit,",fun=mean,na.rm=T,tolerance=0.5)", sep="")}
+    if(i == len){mergeit = paste("big = raster::mosaic(",mergeit,",fun=mean,na.rm=T,tolerance=0.5, filename=outfile, format='GTiff', datatype = 'INT2S',overwrite=T,options=c('COMPRESS=NONE'))", sep="")}
   }
 
-
   print("Creating DEM mosaic")
   eval(parse(text=mergeit))
-  
-  print("Writing DEM mosaic")
-  outfile = file.path(dir,"dem_mosaic.tif")
-  raster::writeRaster(big, outfile, format="GTiff", datatype = "INT2S",overwrite=T,options=c("COMPRESS=NONE"))
   unlink(demfiles)
 
 }

R/reprojectDEM.r

@@ -4,30 +4,33 @@
 #' conform to the properties of an image file prior to using it as an input to the 'MSScvm' masking function
 #' @param demFile filename (character). Full path to DEM file.
 #' @param projRef filename (character). Full path to an image file produced by the \code{\link{MSSunpack}} function.
+#' @param srcNodata numeric. Specify the background value in the input DEM. If there is no background value, use NA (default)
+#' @param dstNodata numeric. Specify the value to represent background pixels in the reprojected DEM. -32768 is the default
 #' @details The DEM file will be adjusted to match the projection and pixel resolution of the 'proRef' image.
 #' @return A GeoTIFF raster file with '_reprojected.tif' replacing the last 4 characters of the input DEM filename
 #' @seealso \code{\link{mosaicDEMs}}
 #' @examples \dontrun{
 #' 
 #' reprojectDEM(demFile = "C:/mss/dem/wrs1_p036r032_dem.tif",
-#'              projRef = "C:/mss/LM10360321973191AAA04/LM10360321973191AAA04_dn.tif")
+#'              projRef = "C:/mss/LM10360321973191AAA04/LM10360321973191AAA04_dn.tif",
+#'              srcNodata= -9999, dstNodata= -32768)
 #' }
 #' @export
 
-reprojectDEM = function(demFile, projRef){
+reprojectDEM = function(demFile, projRef, srcNodata=NA, dstNodata=-32768){
 
   template = raster::raster(projRef) #read in the projection reference image
   reso = raster::xres(template) #get the pixel resolution
   proj = raster::projection(template) #get the projection
 
-  print(paste("Reprojecting file:",demFile))
+  print(paste("Reprojecting",demFile))
   bname = basename(demFile) #get the basename of the dem file
   extension = substr(bname,(nchar(bname)-3),nchar(bname)) #extract the extension of the dem file
   dstfile = sub(extension, "_reprojected.tif", demFile) #define the filename
 
   gdalUtils::gdalwarp(srcfile=demFile, dstfile=dstfile, 
            t_srs=proj, of="GTiff",
-           r="near", multi=T,
+           r="near", multi=T, srcnodata=srcNodata, dstnodata=dstNodata, overwrite=TRUE,
            tr=c(reso,reso), co="INTERLEAVE=BAND") #reproject the image using gdalwarp through gdalUtils
 
 }

man/mosaicDEMs.Rd

@@ -4,14 +4,16 @@
 \alias{mosaicDEMs}
 \title{Create a DEM mosaic from a direcory of DEM's}
 \usage{
-mosaicDEMs(dir, projRef, NAvalue)
+mosaicDEMs(dir, projRef, srcNodata = NA, dstNodata = -32768)
 }
 \arguments{
 \item{dir}{directory name (character). Full path to a directory containing digital elevation model (DEM) files to be mosaiced.}
 
 \item{projRef}{filename (character). Full path to an image file produced by the \code{\link{MSSunpack}} function.}
 
-\item{NAvalue}{numeric. What is the background value of the DEM files in the directory.}
+\item{srcNodata}{numeric. What is the background value of the DEM files in the directory. If there is no background value, use NA (default)}
+
+\item{dstNodata}{numeric. Specify the value to represent background pixels in the mosaic DEM. -32768 is the default}
 }
 \value{
 A GeoTIFF raster file representing the union of all individual DEM files found in the provided directory path. The mosaic file will be written to the
@@ -23,7 +25,7 @@ A helper function to create the large-extent DEM file required by the 'MSScvm' f
 \details{
 The provided directory path should only contain decompressed digital elevation files from the same source (SRTM, NED, GTOPO, etc).
 The function will search the directory and include all files found in the mosaic. It is important that each file have the same background value
-and that it is correctly assigned to the 'NAvalue' parameter, if not, intersection between DEMs could have unexpected results.
+and that it is correctly assigned to the 'srcNodata' parameter, if not, intersection between DEMs could have unexpected results.
 Each individual DEM file will be adjusted to match the projection and pixel resolution of the 'proRef' image.
 Then they will be merged using the mean value of intersecting pixels.
 }
@@ -32,7 +34,7 @@ Then they will be merged using the mean value of intersecting pixels.
 
 mosaicDEMs(dir = "C:/mss/dems",
            projRef = "C:/mss/LM10360321973191AAA04/LM10360321973191AAA04_dn.tif",
-           NAvalue = -32768)
+           srcNodata = -9999, dstNodata= -32768)
 }
 }
 \seealso{