#' 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))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.