BioVeL ESW STACK - ENM Statistical Workflow with raster stack computation

Script

library(rgdal)
library(raster)

size <- as.integer(png_size)

#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());
csv_name <- paste(dest, 'csv', sep="")

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

#create Layerlist
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	
	}
}

#download and read the CSV occurence points
displayPoints <- FALSE
pointFileSeries <-FALSE
if(!is.na(pmatch('http',csv_string))){

#check if one or more files are referenced
pointFiles <- strsplit(csv_string, ",")[[1]];
if (length(pointFiles)> 1 ){
#check if layerlist has the same length as the layerNameList
	if (length(spl) != length(pointFiles)){
		stop("The input csv_list has a different number of entrys to the list of raster layers")	
	} else {
		pointFileSeries <- TRUE
		displayPoints <- TRUE	
	}
} else {
	download.file(url=csv_string, destfile=csv_name, method = "auto", quiet=FALSE)
	pts<-data.frame()
	pts<-read.csv(csv_name, sep=",", header = T)
	coordinates(pts) <- c("decimalLongitude", "decimalLatitude")
	displayPoints <- TRUE	
	}
} else if(!is.na(pmatch('none',csv_string))){
	displayPoints <- FALSE
} else {
	pts<-data.frame()
	con <- textConnection(csv_string)
	pts<-read.csv(con, sep=",", header = T)
	close(con)
	coordinates(pts) <- c("decimalLongitude", "decimalLatitude")
	displayPoints <- TRUE
}


#download the first layer
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)
# output port when needed again = pearson as text-file

#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)

	#if(isTRUE(displayPoints)){
	#	points(coordinates(pts),col = "black", bg = "yellow",pch = 21)
	#}

}
dev.off()

png(filename=stackPNGwithPoints,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)

if(isTRUE(displayPoints)){
	if(isTRUE(pointFileSeries)){
	if(i <= (length(pointFiles))){	
		pointfile <- pointFiles[i]
		download.file(url=pointfile, destfile=csv_name, method = "auto", quiet=FALSE)
		pts<-data.frame()
		pts<-read.csv(csv_name, sep=",", header = T)
		coordinates(pts) <- c("decimalLongitude", "decimalLatitude")
		points(coordinates(pts),col = "black", bg = "yellow",pch = 21)
	}}else{
		points(coordinates(pts),col = "black", bg = "yellow",pch = 21)
	}
	}

}
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

Script

import java.io.StringReader;
import java.io.BufferedReader;

BufferedReader reader = new BufferedReader(new StringReader(csv_content));

    String temp = null;
    StringBuffer text = new StringBuffer();
    while ((temp = reader.readLine()) != null) {
            text.append(temp + "\n");

    };

String csv_string = text.toString().substring(0,text.toString().length()-1);