BioVeL ESW STACK - ENM Statistical Workflow with raster stack computation

Script

library(rgdal)
library(raster)

size <- as.integer(png_size)
#lnames <- as.character(namelist)

#layernames names for the layerstack for covariance first
layernames <- c()
nlistOK <- FALSE
checkNames <- FALSE
if(!isTRUE(!is.na(pmatch('none',namelist)))){	
	if(!is.na(pmatch(',',namelist))){
		stop("The layer namelist is not comma separated or has only one entry")
	}else {
		layernames <- gsub(" ","",namelist,fixed=TRUE)
		namesplit <- strsplit(layernames, ",")[[1]];
		layernames <- c(namesplit)
		checkNames <- TRUE
	}
}

llist <- gsub(" ","",layerlist,fixed=TRUE)
dest <- as.numeric(Sys.time());

if(!is.na(pmatch(',',llist))){
	stop("The input filelist is not comma separated or has only one entry")
}

spl <- strsplit(llist, ",")[[1]];

#check if layerlist has the same length as the layerNameList
if(checkNames){
	if (length(spl) != length(namesplit)){
		stop("The input filelist has a different number of entrys to the list of raster layers")	
	} else {
		nlistOK <- TRUE	
	}
}

layer <- spl[1]
lastst <- strsplit(spl[1], "/")[[1]];
destname <- paste(dest,"_",lastst[length(lastst)], sep="")



if(is.na(pmatch('http',layer))){
	stop("the Layer URL has to start with 'http'")
}

if(isTRUE(toString(Sys.info()['sysname'])=="Windows")){	
	# windows download
	download.file(url=layer, destfile=destname, method = "internal", mode="wb", quiet=FALSE)
} else {
	#linux-mac download
	download.file(url=layer, destfile=destname, method = "auto", quiet=FALSE)
}

x <- readGDAL(destname)
rx <- raster(x)

#layernames names for the layerstack for covariance first
#layernames <- c()
if(!nlistOK){
	layernames <- append(layernames, lastst[length(lastst)])
}

#read rastervalues values from GDALinfo
#GDALinfo
gi1 <- GDALinfo(destname)
#get the datatype from GDALinfo
dataType <- toString(attributes(gi1)$df[,1])
#STOP if the data type is not Byte
if(!isTRUE(dataType == "Byte")){
	stop("one of the input files is not a 'Byte' datatype. Currently are only 'Byte' datatypes allowed")	
}
#get GDAL raster driver
driver <- attributes(gi1)$driver# <- HFA or GTiff
#get noData value from GDALinfo
novalue <- attributes(gi1)$df[,3]
NAvalue(rx) <- novalue

layerstack <-stack(rx)

for(i in 2:length(spl)){
	layer <- spl[i]
	lastst <- strsplit(spl[i], "/")[[1]];
	if(!nlistOK){
		layernames <- append(layernames, lastst[length(lastst)])
	}
	destname <- paste(dest,"_",lastst[length(lastst)], sep="")
	#download.file(url=layer, destfile=destname, method = "auto", quiet=FALSE)

	if(is.na(pmatch('http',layer))){
		next
	}

	if(isTRUE(toString(Sys.info()['sysname'])=="Windows")){	
		# windows download
		download.file(url=layer, destfile=destname, method = "internal", mode="wb", quiet=FALSE)
	} else {
		#linux-mac download
		download.file(url=layer, destfile=destname, method = "auto", quiet=FALSE)
	}

	x <- readGDAL(destname)
	rx <- raster(x)

	#read rastervalues values from GDALinfo
	#GDALinfo
	gi1 <- GDALinfo(destname)
	#get the datatype from GDALinfo
	dataType <- toString(attributes(gi1)$df[,1])
	#STOP if the data type is not Byte
	if(!isTRUE(dataType == "Byte")){
		stop("one of the input files is not a 'Byte' datatype. Currently are only 'Byte' datatypes allowed")	
	}
	#get GDAL raster driver
	driver1 <- attributes(gi1)$driver# <- HFA or GTiff
	if(driver != driver1){
		stop("different file formats. Please use the same format HDF or GTiff for all files")
	}
	#get noData value from GDALinfo
	novalue <- attributes(gi1)$df[,3]
	NAvalue(rx) <- novalue

	if(!isTRUE((origin(rx)==origin(layerstack))&&(extent(rx)==extent(layerstack)) )){
		cropped <- TRUE
		is <- intersect(rx, layerstack)
		rx <-crop (rx, is)
		ry <-crop (layerstack, is)
		#resample rx to ry to get the same origin
		rx <- resample(rx, layerstack, method="ngb")
	}

	layerstack <-stack(layerstack, rx)
}
#layernames names for the layerstack
#layernames <- c()

#download and unzip the shapefile
tmpDir <- tempdir()
shapeDir <- paste(tmpDir,"tdwg_level_4.shp",sep='/')
if(!isTRUE(file.exists(shapeDir))){
	shapeDownload <- "file download"
	url <- "http://biovel.iais.fraunhofer.de/biostif/data/tdwg_level_4.zip"
	file <- basename(url)
	download.file(url, file)
	unzip(file, exdir = tmpDir )	
}
shapeFile <- readOGR(tmpDir,"tdwg_level_4")

#covariance
names(layerstack) <- layernames
covariance_table <- layerStats(layerstack,"cov",na.rm=TRUE)
write.table(covariance_table, file=covariance, sep = ",", quote = FALSE)

#pearson
pearson_table <- layerStats(layerstack,"pearson",na.rm=TRUE)
write.table(pearson_table, file=pearson, sep = ",", quote = FALSE)

#summary file with raster statistics
stacksummary <- paste("Layername","ncell","mean","median","cv","sd","min","max","countNA", sep=",")
for(i in 1:nlayers(layerstack)){
cuRast <- raster(layerstack, layer=i)
fileName <- strsplit(spl[i], "/")[[1]];
##cat(paste(fileName[length(fileName)],",",ncell(cuRast),",",cellStats(cuRast, 'mean'),",",cellStats(cuRast, median),",",cellStats(cuRast,cv),",",cellStats(cuRast, 'sd'),",",cellStats(cuRast,'min'),",",cellStats(cuRast, 'max'),",",cellStats(cuRast, 'countNA'),"\n",sep=""))
ilayer <- paste(fileName[length(fileName)],ncell(cuRast),cellStats(cuRast, 'mean'),cellStats(cuRast, median),cellStats(cuRast,cv),cellStats(cuRast, 'sd'),cellStats(cuRast,'min'),cellStats(cuRast, 'max'),cellStats(cuRast, 'countNA'),sep=",")
stacksummary <- paste(stacksummary,ilayer,sep="\n")
#layernames <- append(layernames, fileName[length(fileName)])
}

#--------------
#compute overall statistics
overall_stat <- paste("Layername","Overall_coverage","Overall_intensity", sep = ",")

for(i in 1:nlayers(layerstack)){
cuRast <- raster(layerstack, layer=i)
fileName <- strsplit(spl[i], "/")[[1]];

count_zero <- 0
if(isTRUE(toString(version["major"][1]) == 2)){
	count_zero <- count(cuRast,0)
} else {
	count_zero <- as.numeric(toString(freq(cuRast)[1,"count"]))
}

overall_cover <- ((ncell(cuRast)-cellStats(cuRast, 'countNA')-count_zero)/ncell(cuRast))*100
overall_intensity <- (cellStats(cuRast, sum)/(ncell(cuRast)-cellStats(cuRast, 'countNA')))
ilayer <- paste(fileName[length(fileName)],overall_cover,overall_intensity,sep=",")
overall_stat <- paste(overall_stat,ilayer,sep="\n")
}

sos <- sum(layerstack, na.rm=FALSE)
maxValue <- tempfile()
if(driver == "GTiff"){
	maxValue <- 254
} else {
	maxValue <- 100
}
sos <- (sos/(nlayers(layerstack)*maxValue))*maxValue
count_zero <- 0
if(isTRUE(toString(version["major"][1]) == 2)){
	count_zero <- count(sos,0)
} else {
	count_zero <- as.numeric(toString(freq(sos)[1,"count"]))
}

overall_cover <- ((ncell(sos)-cellStats(sos, 'countNA')-count_zero)/ncell(sos))*100
overall_intensity <- (cellStats(sos, sum)/(ncell(sos)-cellStats(sos, 'countNA')))
sl <- paste("average_sumlayer",overall_cover,overall_intensity,sep=",")
overall_stat <- paste(overall_stat,sl,sep="\n")
ilayer <- paste("average_sumlayer",ncell(sos),cellStats(sos, 'mean'),cellStats(sos, median),cellStats(sos,cv),cellStats(sos, 'sd'),cellStats(sos,'min'),cellStats(sos, 'max'),cellStats(sos, 'countNA'),sep=",")
stacksummary <- paste(stacksummary,ilayer,sep="\n")
layerstack <-stack(layerstack, sos)
layernames <- append(layernames, "average_sumlayer")

#create pngs
png(filename=PNGoverview,height=size/2, width=size, bg="white")
cols <- c(rgb(1.0,1.0,1.0),rgb(0.996,0.898,0.851),rgb(0.988,0.682,0.569),rgb(0.984,0.416,0.290),rgb(0.871,0.176,0.149),rgb(0.647,0.059,0.082))

names(layerstack) <- layernames
fun <- function(){plot(shapeFile, add=TRUE)}
plot(layerstack,col=cols,addfun=fun)
dev.off()


png(filename=stackPNG,height=0.94*size*nlayers(layerstack), width=size,bg="white")
par(mfrow=c(nlayers(layerstack),1))
cols <- c(rgb(1.0,1.0,1.0),rgb(0.996,0.898,0.851),rgb(0.988,0.682,0.569),rgb(0.984,0.416,0.290),rgb(0.871,0.176,0.149),rgb(0.647,0.059,0.082))
breakpoints <- tempfile()
for(i in 1:nlayers(layerstack)){
	if(driver == "GTiff"){
		breakpoints <- c(0,2,63,127,190,254)
	}else{
		breakpoints <- c(0,1,25,50,75,100)
	}
	plot(layerstack,i,breaks=breakpoints,col=cols,alpha=0.7)
	plot(shapeFile, add=TRUE)
}
dev.off()

#store stack
if(driver == "HFA"){
	writeRaster(layerstack, filename='stackImg.img', format="HFA", overwrite=TRUE)
	file.rename('stackImg.img','stackImg.png')
} else {
	writeRaster(layerstack, filename='stackImg.tif', format="GTiff", overwrite=TRUE)
	file.rename('stackImg.tif','stackImg.png')
}

R Server