BioVeL ESW DIFF Basic

None

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

Name	Description
layerlist	comma separated list of image URLs e.g. http://biovel.googlecode.com/svn/trunk/esw/sample_data/0ubESU.tif,http://biovel.googlecode.com/svn/trunk/esw/sample_data/mXGm7g.tif,http://biovel.googlecode.com/svn/trunk/esw/sample_data/NtNUdB.tif,http://biovel.googlecode.com/svn/trunk/esw/sample_data/wVM8bF.tif
png_size	size of the PNG overviews and thumbnails in pixel
namelist	comma separated list of user layernames for using in in the PNG overlays and the matrices, e.g.: layername_1,layername_2,layername_3, layername_4 the list must have the same order as the input raster list and the same number of elements when own user layer names are not needed, please type "none" as input value (this is the standard value), in this case the raster names will be used as layernames
occurrence_csv	URL to the CSV File with occurrence points coming from ENM. If you don't want to display any points, type "none" at the begin, else paste a CSV URL or the content from a CSV or a CSV file as input It is possible to display different pointsets to the corresponding PNGs. In that case, the point URLs should be entered comma separated in the same order as the input files. When only one URL is entered, the pointset of that URL will be added to all PNGs. example: http://biovel.googlecode.com/svn/trunk/esw/sample_data/ED_4.csv,http://biovel.googlecode.com/svn/trunk/esw/sample_data/AT_1.csv,http://biovel.googlecode.com/svn/trunk/esw/sample_data/CG_3.csv,http://biovel.googlecode.com/svn/trunk/esw/sample_data/AM_2.csv

Name	Description
stackPNG	PNG overview of the input current layer with a background country borders based on the TDWG level4 layer Used color ramp: ERDASImagine file: breakpoints = (0,1,25,50,75,100) GeoTIFF file breakpoints = c(0,2,63,127,190,254) colors = 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))
covariance	computed covariance matrix from the input raster layers.
rastersummary	This file contains the summary for the current, prediction and the difference layer: - Number of cells (ncell), - Mean and median value (mean,median) - Coefficient of variation (cv) - Standard deviation (sd) - current minimal and maximal values (min,max) - Number of noData values in file (countNA)
overall_statistics	This file contains the changes computation in the submitted files: Overall coverage in the input raster files and the average sum layer overall_coverage = ((ncell(x)-cellStats(x, 'countNA')-count(x,0))/(ncell(x)-cellStats(x, 'countNA')))*100 Overall intensity absolute and in percent in the raster files and the average sum layer overall_intensity = (cellStats(x, sum)/(ncell(x)-cellStats(x, 'countNA')))
stackImg	stack multilayer file as GeoTIFF or ERDASImagine file contains all input raster layers and the computed average sum layer
PNGoverview	PNG overview of the input current layer with a background country borders based on the TDWG level4 layer Used color ramp: ERDASImagine file: breakpoints = (0,1,25,50,75,100) GeoTIFF file breakpoints = c(0,2,63,127,190,254) colors = 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))

Source	Sink
layerlist	Rshell:layerlist
png_size	Rshell:png_size
namelist	Rshell:namelist
parse_csv:csv_string	Rshell:csv_string
occurrence_csv	parse_csv:csv_content
Rshell:stackPNG	stackPNG
Rshell:covariance	covariance
Rshell:stacksummary	rastersummary
Rshell:overall_stat	overall_statistics
Rshell:stackImg	stackImg
Rshell:PNGoverview	PNGoverview

BioVeL ESW DIFF Basic

Preview

Run

Run this Workflow in the Taverna Workbench...

Workflow Components

Script

R Server

Script

Reviews (0)

Comments (0)

Other workflows that use similar services (0)