BioVeL ESW STACK - ENM Statistical Workflow with raster stack computation

Script

##############################################
# FUNCTIONS	        	 	         #
##############################################

## function for getting the correct aspect ratio of the raster in order to be plotted correctly - 15/05/2014
getAsp <- function(rasterObj) {
  if (raster:::.couldBeLonLat(rasterObj, warnings=FALSE)){ 
    ym_rasterObj <- mean(c(ymax(rasterObj), ymin(rasterObj)))
    asp <- 1/cos((ym_rasterObj * pi)/180) 
  }
  else { 
    asp <- 1 
  }
  asp
}

## function for printing the layer stack object received as parameter with the layer names, 
## shape file and occurrences points also passed as parameters
plotLayerStack <- function(layerstackObj,layernamesObj,shapeFileObj,driverObj,pointsSeriesObj) {
  if(driverObj == "GTiff"){
    breakpoints <- c(0,2,63,127,190,254)
  }else{
    breakpoints <- c(0,1,25,50,75,100)
  }
  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))
  
  for(i in 1:nlayers(layerstackObj)){
    currentLayer <- subset(layerstackObj, i, drop=TRUE)
    plot(currentLayer,main=layernamesObj[i],breaks=breakpoints,col=cols,alpha=0.7,asp=getAsp(currentLayer)) 
    plot(shapeFile, add=TRUE)	
    
    #Print points	  
    if (length(pointsSeriesObj) > 0 ){	
      if(i <= length(pointsSeriesObj)){	
        occurencePoints<-pointsSeriesObj[[i]]
        coordinates(occurencePoints) <- c("decimalLongitude", "decimalLatitude")
        points(coordinates(occurencePoints),col = "black", bg = "yellow",pch = 21)
      }else{
        #summatory layer
        for(j in 1:length(pointsSeriesObj)){
          occurencePoints<-pointsSeriesObj[[j]]
          coordinates(occurencePoints) <- c("decimalLongitude", "decimalLatitude")		
          points(coordinates(occurencePoints),col = "black", bg = "yellow",pch = 21)
        }
      }
	}
  }
}

##############################################
# MAIN BODY OF THE 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 and check if it has the same length as the layerNameList
#--------------------------------------------------------------
spl <- strsplit(llist, ",")[[1]];

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
#--------------------------------------------------------------
pointsSeries <- list()
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 {
      for(i in 1:length(pointFiles)){
        download.file(url=pointFiles[i], destfile=csv_name, method = "auto", quiet=FALSE)
        pointsSeries[[length(pointsSeries)+1]] <- read.csv(csv_name, sep=",", header = T)
      }
    }
  } else {
    download.file(url=csv_string, destfile=csv_name, method = "auto", quiet=FALSE)
    pointsSeries[[length(pointsSeries)+1]] <- read.csv(csv_name, sep=",", header = T)
  }
} else if(!is.na(pmatch('none',csv_string))){
  pointsSeries <- list()
} else {
  con <- textConnection(csv_string)
  pointsSeries[[length(pointsSeries)+1]]<-read.csv(con, sep=",", header = T)
  close(con)
}


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

#--------------------------------------------------------------
# Download the rest of the layers checking that have the same 
# format as the first layer. Also cropping layers when necessary
#--------------------------------------------------------------
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="")
  
  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)
}

#--------------------------------------------------------------
#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"
  url <- "http://biovel.googlecode.com/svn/trunk/esw/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, col.names=NA)

#--------------------------------------------------------------
#pearson
#--------------------------------------------------------------
#pearson_table <- layerStats(layerstack,"pearson",na.rm=TRUE)
#write.table(pearson_table, file=pearson, sep = ",", quote = FALSE, col.names=NA)
# 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)  
  ilayer <- paste(layernames[i],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")
}

#--------------------------------------------------------------
#compute overall statistics
#--------------------------------------------------------------
overall_stat <- paste("Layername","Overall_coverage","Overall_intensity", sep = ",")
for(i in 1:nlayers(layerstack)){
  cuRast <- raster(layerstack, layer=i)
  
  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)-cellStats(cuRast, 'countNA')))*100
  overall_intensity <- (cellStats(cuRast, sum)/(ncell(cuRast)-cellStats(cuRast, 'countNA')))
  ilayer <- paste(layernames[i],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)-cellStats(sos, 'countNA')))*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
#--------------------------------------------------------------
numColumnsOverView <- ceiling(sqrt(nlayers(layerstack)))
numRowsOverView <- ceiling(nlayers(layerstack)/numColumnsOverView)
weightOverview <- size
heightOverview <- ceiling((size/numColumnsOverView)*numRowsOverView)

png(filename=PNGoverview,width=weightOverview,height=heightOverview,bg="white")
par(mfrow=c(numRowsOverView,numColumnsOverView),mar=c(5.1,2.1,2.1,5.1))
plotLayerStack(layerstack,layernames,shapeFile,driver,pointsSeries)
dev.off()

png(filename=stackPNG,height=size*nlayers(layerstack), width=size,bg="white")
par(mfrow=c(nlayers(layerstack),1))
plotLayerStack(layerstack,layernames,shapeFile,driver,pointsSeries)
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);