R/BioCroOutput_to_shapefile.R

In djaiswal/BioCroRegional: Tools to perform Regional Simulations of BioCro

#'  This functions taked a data frame of BioCro Output containing three columns-Longitude, Latitude, and Yield and a Shapefile of polygon and add an additional column of given name "Variable" in the shapefile
#'  
#'  @param shapefile
#'  @param BioCroOutPut
#'  @param VariableName will be the name of new field in shapefile
#'  @param GridResolution is used in creating prediction raster Grid for interpolating BioCro Output
#'  @param InterpolationMethod- I can currently use only thin plate spline method provided in R package fields. I can extend this later providing options for kriging etc.
#'  @return A list containing mean value of zonal analysesi and raster of the BioCroOutput in ascii format
#'  @export

BioCroOutput_to_shapefile <- function(InputShapefile, BioCroOutput, VariableName, GridResolution = 0.0166, InterpolationMethod = "Tps"){
  
  #Check inputs
  
  #1
  if(class(InputShapefile)!="SpatialPolygonsDataFrame"){stop("class of Inputshapefile is not SpatialPolygonsDataFrame. It must be an output of readOGR function of rgdal package")}
  #2
  if(!("lon" %in% names(BioCroOutput))) {stop ("lon column is missing from BioCroOutput")}
  if(!("lat" %in% names(BioCroOutput))) {stop ("lat column is missing from BioCroOutput")}
  if(!("Yield" %in% names(BioCroOutput))) {stop ("Yield column is missing from BioCroOutput")}
  #3
  if(class(VariableName)!="character") {stop ("Variable name to add as new field is not string")}
     
  #find extent of shapefile to create raster
  brazilshapeextent <- extent(InputShapefile)
     
  #create a raster layer object having the same extent as Brazil
  predictedraster <- raster(brazilshapeextent)
     
  #set the resolution of raster object 
  res(predictedraster) <- GridResolution 
  
  #Interpolating Output to desired Grid
  if(InterpolationMethod=="Tps"){
  biocro.coordinates <- data.frame(lon = BioCroOutput$lon,lat = BioCroOutput$lat)
  biocro.yield <- BioCroOutput$Yield
  tps <- Tps(biocro.coordinates,biocro.yield)
  }else {stop ("Interpolation method is not implemented yet")}
  
  #predict after interpolation
  predictedraster <- interpolate(predictedraster,tps)
  
  #mask the prediction grid using Shapefile
  predictedraster <- mask(predictedraster, InputShapefile)
  
  #extracting values from interpolated raster to perform zonal statistics
  rastertoshape <- extract(predictedraster,InputShapefile,fun=mean)
  rastertoshape <- unlist(lapply(rastertoshape, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA ))
  
  #Combining predicted yield to original shape file
  InputShapefile$newfield <- rastertoshape
  
  #Change name of newly added field
  names(InputShapefile)[length(names(InputShapefile))] <- VariableName
 
  #return the InputShapefile containing additional NewVariable and raster file
  return(list(shapefile = InputShapefile,raster = predictedraster))
}