R/percent_cover_parallel.R

In Tomhigg/fracCover: Generate and Validate Fractional Coverage products using Landsat/MODIS

#' Percentage Cover Calculations 

#'@description
#' This script reads training data from the CSV file created using the "create_csv" 
#' script.  The script then uses the X and Y coordinates from the training data file to select
#' the pixel values (predictor values) for each sample point in the input image. The predictor 
#' values and the percent cover data from the training data file (response variable) are 
#' combined and used as input to the random forests model. After the model is created percent 
#' cover predictions are made on the input image to create an output image with percent cover 
#' values ranging from 0 to 1. 

#'@param no_cores number of cores to implement on 
#'@param pointdata Data frame of training data from the rf_csv function, does not need to be projected
#'@param LS.stack A stack of Landsat data, ideally generated by the team LUCC package 
#'@param LS.no.data No data value for the Landsat image, normally 0 
#'@param outImage. Name and path for the output classified image, use NA if not needed 
#'@param point_CRS a CRS object for the training data 
#'@return a raster image classified into percentage coverage



percent_cover_parallel<-function(no_cores,inImage, pointdata,outImage,LS.no.data,point_CRS){
# Start processing
print("Set variables and start processing")
startTime <- Sys.time()
cat("Start time", format(startTime),"\n")


cl <- makeCluster(spec = no_cores)
# Register the cluster with foreach:
registerDoParallel(cl)

# Load the moderate resolution image
for (b in 1:nlayers(inImage)) { NAvalue(inImage@layers[[b]]) <- LS.no.data}

#load point data
#point.csv <- read.csv(file = pointData,header = TRUE)
xy <- SpatialPoints(pointdata[,1:2],proj4string = point_CRS)
xy <-spTransform(x = xy,(crs(inImage)))
response <- as.numeric(pointdata[,3])


# Get pixel DNs from the input image for each sample point
print("Getting the pixel values under each point")
trainvals <- cbind(response, extract(inImage, xy)) 

# Remove NA values from trainvals
trainvals_no_na <- na.omit(trainvals)

# Run Random Forest
print("Starting to calculate random forest object")

randfor <- randomForest(response ~. , data=trainvals_no_na)


print("Enter the Matrix")
randomForest_raster_predict <- function(inraster,rfModel,...)
{
  # We need to load randomForest (note this is only
  # loaded once per worker on most cluster setups):
  require(randomForest)
  
  # First, preserve the names:
  band_names <- dimnames(inraster)[3][[1]]
  
  # This will “flatten” the array to a matrix (we lose the names here):
  inraster_matrix <- inraster
  dim(inraster_matrix) <- c(dim(inraster)[1]*dim(inraster)[2],dim(inraster)[3])
  
  # Now, let’s coerce this to a data.frame:
  inraster.df <- as.data.frame(inraster_matrix)
  
  # We need to re-set the names because randomForest requires it:
  names(inraster.df) <- band_names
  
  # Now we can use this in a predict statement:
  out_predictions <- predict(rfModel,inraster.df)
  
  # We must return this as a numeric array (a raster cannot
  # use character classes), and set the dims to match the input.
  # Also, right now rasterEngine requires ALL outputs to be double
  # precision floating point, so we cannot use “as.numeric” on the
  # factor, because that will return an integer.
  
  out_predictions_array <- array(as.double(out_predictions),dim=c(dim(inraster)[1:2],1))
  
  return(out_predictions_array)
}

# Now, rasterEngine. Notice we pass the randomForest model along
# via the args= setting.
system.time(
  
  cover_output_array <-
    rasterEngine(
      # Match the variable name in the function to the raster:
      inraster=inImage,
      # Assign the function:
      fun=randomForest_raster_predict,
      args=list(rfModel= randfor)
    ))
print("Unplug")
cover_output <- setMinMax(cover_output_array)


stopCluster(cl) # Stops the cluster

if (outImage !=NA) {
  writeRaster(x = cover_output,filename = outImage, format="GTiff", overwrite=TRUE)
} 


# Calculate processing time
timeDiff <- Sys.time() - startTime
cat("Processing time", format(timeDiff), "\n")
return(cover_output)

}