BioVeL ESW DIFF Basic

Script

#------------------------------
#load required librarys, they must be installed on the R instnce
library(rgdal)
library(raster)

#-------------------------------
#create unique names with a timestamp (Sys.time)
creation_time <- Sys.time()
dest <- as.numeric(creation_time)
destname1 <- paste(dest, '_x_raster', sep="") 
destname2 <- paste(dest, '_y_raster', sep="")
#csv_name <- paste(dest, 'csv', sep="")

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

shapeDownload <- "OK"
#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")

#download and read the CSV occurence points
displayPoints <- FALSE
#if(!is.na(pmatch('http',csv_url))){
#	download.file(url=csv_url, 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_url))){
#	displayPoints <- FALSE
#} else {
#	pts<-data.frame()
#	con <- textConnection(csv_url)
#	pts<-read.csv(con, sep=",", header = T)
#	close(con)
#	coordinates(pts) <- c("decimalLongitude", "decimalLatitude")
#	displayPoints <- TRUE
#}

#---------------------------
#variables
normalize <- FALSE
normalizeBoth <- FALSE
#I can't get the really max values, only the current values from the rasterimage
enm_tif_max <- as.numeric(254.0)
enm_img_max <- as.numeric(100.0)
novalue <- FALSE
diffxy <- tempfile()
rx_normalized <- tempfile()
ry_normalized <- tempfile()
cropped <- FALSE
output_format <- "HFA"
size <- as.integer(png_size)

#----------------------------
#read the images
x <- readGDAL(destname1)
y <- readGDAL(destname2)
# and convert to raster
rx <- raster(x)
ry <- raster(y)

#---------------------------
#read rastervalues values from GDALinfo
#GDALinfo
gi1 <- GDALinfo(destname1)
gi2 <- GDALinfo(destname2)
#check the datatype, only "Byte" with a valueRange from 0 - 255 is allowed
#TODO: check to work with other datatypes, when they are different, normalize to 100 or to 1

#get the datatype from GDALinfo
dataType1 <- toString(attributes(gi1)$df[,1])
dataType2 <- toString(attributes(gi2)$df[,1])
#STOP if the data type is not Byte
if(!isTRUE((dataType1 == "Byte")&&(dataType2 == "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
driver2 <- attributes(gi2)$driver# <- HFA or GTiff
#check wheather the input files are IMG or TIFF
if(!((isTRUE(driver1 == "HFA") || isTRUE(driver1 == "GTiff"))
	&&(isTRUE(driver2 == "HFA") || isTRUE(driver2 == "GTiff")))){
	stop("one of the input files is not a IMG or GeoTIFF file")	
}
#check if the driver is GTiff or different drivers and hat to be normalized the maxvalue to 100
if(isTRUE(driver1 == "GTiff") && isTRUE(driver2 == "GTiff")){
	normalizeBoth <- TRUE
	rx_normalized <- (rx/enm_tif_max)*100	
	ry_normalized <- (ry/enm_tif_max)*100	
} else if(!isTRUE(driver1 == driver2)){
	#in case of mix mode HFA and GTiff
	normalize <- TRUE
	if(isTRUE(driver1 == "GTiff")){
		rx_normalized <- (rx/enm_tif_max)*100	
	}else{
		rx_normalized <- rx
	}
	if(isTRUE(driver2 == "GTiff")){
		ry_normalized <- (ry/enm_tif_max)*100	
	}else{
		ry_normalized <- ry
	}
} else {
	rx_normalized <- rx
	ry_normalized <- ry	
}

#get noData value from GDALinfo
#novalue1 <- paste(attributes(gi1)$df[,3],".0",sep="")
novalue1 <- attributes(gi1)$df[,3]
novalue2 <- attributes(gi2)$df[,3]
if(isTRUE(normalize)){
	novalue <- as.numeric(101.0)
} else if(isTRUE(novalue1==novalue2)){
	novalue <- as.numeric(novalue1)
} else {
	stop("something is wrong with noData values in both files - noData values are different")	
}

#intersect both layers => create new layers with a common extent and resample the files if the origin or extent is different, else error
if(!isTRUE((origin(rx)==origin(ry))&&(extent(rx)==extent(ry)) )){
	cropped <- TRUE
	is <- intersect(rx, ry)
	rx <-crop (rx, is)
	ry <-crop (ry, is)
	#resample rx to ry to get the same origin
	rx <- resample(rx, ry, method="ngb")

	if(isTRUE(normalize) || isTRUE(normalizeBoth)){
		rx_normalized <-crop (rx_normalized, is)
		ry_normalized <-crop (ry_normalized, is)
		#resample rx to ry to get the same origin
		rx_normalized <- resample(rx_normalized, ry_normalized, method="ngb")
 	}
}
#diff from raster rx and raster ry -> 
if(isTRUE(normalize)){
	diffxy <- ry_normalized - rx_normalized
}else {
	diffxy <- (ry-rx)	
}

#for further computations and creating tables the "NA" values will be replaces with real NA values, 101 for HFA or 254 for GTiff
NAvalue(rx) <- novalue
NAvalue(ry) <- novalue
NAvalue(diffxy) <- novalue
if(isTRUE(normalize)|| isTRUE(normalizeBoth)){
	NAvalue(rx_normalized) <- novalue
	NAvalue(ry_normalized) <- novalue
}

#summary file with raster statistics
sum2 <- tempfile()
sum3 <- tempfile()
sum4 <- tempfile()

sum1 <- paste("XX","ncell","mean","median","cv","sd","min","max","countNA",sep = ",")
if(isTRUE(normalize)){
	sum2 <- paste("currentLayer",ncell(rx_normalized),cellStats(rx_normalized, stat='mean', na.rm=TRUE),cellStats(rx_normalized, median),cellStats(rx_normalized,cv),cellStats(rx_normalized, stat='sd'),cellStats(rx_normalized,stat='min'),cellStats(rx_normalized, stat='max'),cellStats(rx_normalized, stat='countNA'),sep=",")
	sum3 <- paste("predictionLayer",ncell(ry_normalized),cellStats(ry_normalized, stat='mean', na.rm=TRUE),cellStats(ry_normalized, median),cellStats(ry_normalized,cv),cellStats(ry_normalized, stat='sd'),cellStats(ry_normalized, stat='min'),cellStats(ry_normalized, stat='max'),cellStats(ry_normalized, stat='countNA'),sep=",")
	sum4 <- paste("diffLayer",ncell(diffxy),cellStats(diffxy, stat='mean', na.rm=TRUE),cellStats(diffxy, median),cellStats(diffxy, cv),cellStats(diffxy, stat='sd'),cellStats(diffxy,stat='min'),cellStats(diffxy, stat='max'),cellStats(diffxy, stat='countNA'), sep=",")
} else{
	sum2 <- paste("currentLayer",ncell(rx),cellStats(rx, stat='mean', na.rm=TRUE),cellStats(rx, median),cellStats(rx,cv),cellStats(rx, stat='sd'),cellStats(rx,stat='min'),cellStats(rx, stat='max'),cellStats(rx, stat='countNA'),sep=",")
	sum3 <- paste("predictionLayer",ncell(ry),cellStats(ry, stat='mean', na.rm=TRUE),cellStats(ry, median),cellStats(ry,cv),cellStats(ry, stat='sd'),cellStats(ry,stat='min'),cellStats(ry, stat='max'),cellStats(ry, stat='countNA'),sep=",")
	sum4 <- paste("diffLayer",ncell(diffxy),cellStats(diffxy, stat='mean', na.rm=TRUE),cellStats(diffxy, median),cellStats(diffxy, cv),cellStats(diffxy, stat='sd'),cellStats(diffxy,stat='min'),cellStats(diffxy, 'max'),cellStats(diffxy, stat='countNA'), sep=",")
}
rastersummary <- paste(sum1, sum2, sum3, sum4, sep = "\n")

#compute overall statistics
#NOTE count value -> count(rx,0) is deprecated in R version >3
#use freq instead
count_zero_rx <- 0
count_zero_ry <- 0

if(isTRUE(toString(version["major"][1]) == 2)){
	# for R version 2
	count_zero_rx <- count(rx,0)
	count_zero_ry <- count(ry,0)
} else {
	#for R version 3+
	count_zero_rx <- as.numeric(toString(freq(rx)[1,"count"]))
	count_zero_ry <- as.numeric(toString(freq(ry)[1,"count"]))
}

#compute overall cover
overall_coverX <- ((ncell(rx)-cellStats(rx, 'countNA')-count_zero_rx)/(ncell(rx)-cellStats(rx, 'countNA')))*100
overall_coverY <- ((ncell(ry)-cellStats(ry, 'countNA')-count_zero_ry)/(ncell(ry)-cellStats(rx, 'countNA')))*100

#compute overall intensity
overall_intensityX <- (cellStats(rx, sum)/(ncell(rx)-cellStats(rx, 'countNA')))
overall_intensityY <- (cellStats(ry, sum)/(ncell(ry)-cellStats(ry, 'countNA')))
diff_intensity <- overall_intensityY - overall_intensityX 
#compute nomalized intensity to have comperable values in percent
diff_intensity_normalized <- tempfile()
if(isTRUE(normalize) || isTRUE(normalizeBoth)){
	overall_intensityX_normalized <- (cellStats(rx_normalized, sum)/(ncell(rx_normalized)-cellStats(rx_normalized, 'countNA')))
	overall_intensityY_normalized <- (cellStats(ry_normalized, sum)/(ncell(ry_normalized)-cellStats(ry_normalized, 'countNA')))
	diff_intensity_normalized <- overall_intensityY_normalized - overall_intensityX_normalized
} else {
	overall_intensityX_normalized <- overall_intensityX
 	overall_intensityY_normalized <- overall_intensityY	
	diff_intensity_normalized <- overall_intensityY_normalized - overall_intensityX_normalized
}
diff_intensity_normalized <- overall_intensityY_normalized - overall_intensityX_normalized

changes1 <- tempfile()
changes2 <- tempfile()
changes3 <- tempfile()
if(isTRUE(normalize)){
	changes1 <- paste("XX","Overall_coverage","in_the_Raster_output","DIFF_coverage","Overall_intensity","absolute_","Overall_intensity","percent_","DIFF_intensity_percent",sep=",")
	changes2 <- paste("species","2012","2050","2050-2012","2012","2050","2012_%","2050_%","2050-2012_%",sep=",")
	changes3 <- paste(label,overall_coverX,overall_coverY,overall_coverY-overall_coverX,overall_intensityX,overall_intensityY,overall_intensityX_normalized,overall_intensityY_normalized,diff_intensity_normalized,sep=",")
}else if(isTRUE(normalizeBoth)){
	changes1 <- paste("XX","Overall_coverage","in_the_Raster_output","DIFF_coverage","Overall_intensity","absolute_","Overall_intensity","percent_","DIFF_intensity_absolute","DIFF_intensity_percent",sep=",")
	changes2 <- paste("species","2012","2050","2050-2012","2012","2050","2012_%","2050_%","2050-2012","2050-2012_%",sep=",")
	changes3 <- paste(label,overall_coverX,overall_coverY,overall_coverY-overall_coverX,overall_intensityX,overall_intensityY,overall_intensityX_normalized,overall_intensityY_normalized,diff_intensity,diff_intensity_normalized,sep=",")	
}else{
	changes1 <- paste("XX","Overall_coverage","in_the_Raster_output","DIFF_coverage","Overall_intensity","absolute_","Overall_intensity","percent_","DIFF_intensity_absolute","DIFF_intensity_percent", sep = ",")
	changes2 <- paste("species","2012","2050","2050-2012","2012","2050","2012_%","2050_%","2050-2012","2050-2012_%",sep=",")
	changes3 <- paste(label,overall_coverX,overall_coverY,overall_coverY-overall_coverX,overall_intensityX,overall_intensityY,overall_intensityX_normalized,overall_intensityY_normalized,diff_intensity,diff_intensity_normalized,sep=",")	
}
changes <- paste(changes1,changes2,changes3, sep = "\n")


#export of PNG overviews
png(filename=currentLayerPng,height=size, width=size, bg="white")
breakpoints <- tempfile()
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))
if(isTRUE(normalizeBoth)){
	breakpoints <- c(0,2,63,127,190,254)
}else{
	breakpoints <- c(0,1,25,50,75,100)
}
if(isTRUE(normalize)){
	plot(rx_normalized ,main='currentLayer',breaks=breakpoints,col=cols, alpha=0.7)
} else {
	plot(rx ,main='currentLayer',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=predictionLayerPng,height=size, width=size, bg="white")
#plot(rep_ry,main='predictionLayer',breaks=breakpoints,col=cols, alpha=0.7)
if(isTRUE(normalize)){
	plot(ry_normalized ,main='predictionLayer',breaks=breakpoints,col=cols, alpha=0.7)
} else {
	plot(ry ,main='predictionLayer',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()
breakpoints_diff <- tempfile()
if(isTRUE(normalizeBoth)){
	breakpoints_diff <- c(-254,-190,-127,-63,-2,2,63,127,190,254)
}else{
	breakpoints_diff <- c(-100,-70,-50,-25,-1,1,25,50,70,100)
}
cols_diff <- c(rgb(0.396,0.310,0.635),rgb(0.196,0.533,0.741),rgb(0.4,0.761,0.647),rgb(0.671,0.867,0.643),rgb(1.0,1.0,1.0),rgb(0.992,0.682,0.380),rgb(0.957,0.427,0.263),rgb(0.835,0.243,0.310),rgb(0.620,0.004,0.259))
png(filename=diffLayerPng,height=size, width=size, bg="white")
plot(diffxy,main='diffLayer',breaks=breakpoints_diff,col=cols_diff, alpha=0.7)
#plot(shapeFile, add=TRUE)
#if(isTRUE(displayPoints)){
#	points(coordinates(pts),col = "black", bg = "yellow",pch = 21)
#}
dev.off()

#store files
#-------------------------------------------------
#store input raster in a table/Arcgrid format

## NOTE: the usage of the ARCGRID format is necessary because it is only possible to pass text files between Taverna
## workflowas, any other files will be stringified and damaged
if(isTRUE(normalize)){
	writeRaster(rx_normalized, filename='x_matrix.asc', format="ascii", overwrite=TRUE)
		file.rename('x_matrix.asc','x_matrix.png')
	writeRaster(ry_normalized, filename='y_matrix.asc', format="ascii", overwrite=TRUE)
		file.rename('y_matrix.asc','y_matrix.png')
	
}else {
	writeRaster(rx, filename='x_matrix.asc', format="ascii", overwrite=TRUE)
		file.rename('x_matrix.asc','x_matrix.png')
	writeRaster(ry, filename='y_matrix.asc', format="ascii", overwrite=TRUE)
		file.rename('y_matrix.asc','y_matrix.png')
}

#store diffresult as a TIFF/IMG file and a Arcgrid file
if(isTRUE(normalizeBoth)){
	output_format <- "GTiff"
	writeRaster(diffxy, filename='diffImg.tif', format="GTiff", overwrite=TRUE)
		file.rename('diffImg.tif','diffImg.png')
} else {
	output_format <- "HFA"
	writeRaster(diffxy, filename='diffImg.img', format="HFA", overwrite=TRUE)
		file.rename('diffImg.img','diffImg.png')	
}
writeRaster(diffxy, filename='diffascii.asc', format="ascii", overwrite=TRUE)
#write.asciigrid(rep_diff, fname='diffascii.asc',na.value = novalue)
	file.rename('diffascii.asc','diffascii.png')


#write logfile with process informations
log_date <- paste("computation start: ",creation_time, sep="")
inFile1 <-tempfile()
values_in1 <- tempfile()
values_in1_norm <- tempfile()
inFile2 <-tempfile()
values_in2 <- tempfile()
values_in2_norm <- tempfile()
outFile <-tempfile()
values_out <- tempfile()
isCropped <- tempfile()
isNormalized <- tempfile()

if(driver1=="HFA"){
	inFile1 <- 'File 1 input format: ERDASImagine (.img)'
	values_in1 <- paste("value range: 0 <> 100, noData 101, current values: min ",cellStats(rx,'min'), ", max: ", cellStats(rx,'max'),sep="")
	values_in1_norm <- ""
}else{
	inFile1 <- 'File 1 input format: GeoTIFF (.tif)'
	values_in1 <- paste("value range: 0 <> 254, noData 255, current values: min ",cellStats(rx,'min'), ", max: ", cellStats(rx,'max'),sep="")

	if(isTRUE(normalize)){
		values_in1_norm <- paste("value range normalized: 0 - 100 , noData 101, current values: min ",cellStats(rx_normalized,'min'), " max: ", cellStats(rx_normalized,'max'),sep="")	
	} else {
		values_in1_norm <- ""
	}
}
if(driver2=="HFA"){
	inFile2 <- 'File 2 input format: ERDASImagine (.img)'
	values_in2 <- paste("value range: 0 <> 100, noData 101, current values: min ",cellStats(ry,'min'), ", max: ", cellStats(ry,'max'),sep="")
	values_in2_norm <- ""
}else{
	inFile2 <- 'File 2 input format: GeoTIFF (.tif)'
	values_in2 <- paste("value range: 0 <> 254, noData 255, current values: min ",cellStats(ry,'min'), ", max: ", cellStats(ry,'max'),sep="")

	if(isTRUE(normalize)){
		values_in2_norm <- paste("value range normalized: 0 - 100 , noData 101, current values: min ",cellStats(ry_normalized,'min'), ", max: ", cellStats(ry_normalized,'max'),sep="")	
	} else {
		values_in2_norm <- ""
	}
}
if (isTRUE(normalizeBoth)){
	outFile <- "output file format: GeoTIFF (.tif)"
	values_out <- paste("value range: -254 <> 254, noData: 255, current values: min ",cellStats(diffxy,'min'), ", max: ", cellStats(diffxy,'max'),sep="")
} else {
	outFile <- "output file format: ERDASImagine (.img)"
	values_out <- paste("value range: -100 <> 100, noData: 101, current values: min ",cellStats(diffxy,'min'), ", max: ", cellStats(diffxy,'max'),sep="")
}

if (isTRUE(cropped)){
	isCropped <- "input files cropped: YES"
} else {
	isCropped <- "input files cropped: NO"
}

if (isTRUE(normalize)){
	isNormalized <- "input files nomalized: YES"
} else {
	isNormalized <- "input files normalized: NO"
}
logfile <- paste("---------------",log_date,"------", inFile1, values_in1, values_in1_norm,"------", inFile2, values_in2, values_in2_norm,"------", isCropped, isNormalized,"------", outFile, values_out,"-----","point display", displayPoints, sep="\n")

R Server