BioVeL ESW DIFF - ENM Statistical Workflow with raster difference computation

Script

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

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

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

#download and unzip the shapefile
shapeDownload <- "OK"
tmpDir <- tempdir()
shapeDir <- paste(tmpDir,"tdwg_level_4.shp",sep='/')
if(!isTRUE(file.exists(shapeDir))){
   shapeDownload <- "file download"
   url <- "http://biovel.github.io/esw/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 raster values 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)	
}

#-------------------------------------------------
#Centre of Gravity or Mass calculations for spatial data 
#-------------------------------------------------
x_centre <- COGravity(rx) #or COGravity(x) 
y_centre <- COGravity(ry) #or COGravity(y)

x_COGx <- x_centre[1]
x_COGy <- x_centre[3]

y_COGx <- y_centre[1]
y_COGy <- y_centre[3]

if (x_COGx != y_COGx || x_COGy!=y_COGy){
   distance_shift <- SDMTools::distance(x_COGy, x_COGx, y_COGy, y_COGx, bearing=TRUE)
   dist_km <- distance_shift[1,5]/1000  
}else {
   dist_km <- 0
}

displayMassPoints <- FALSE
if(!is.na(pmatch('true',center_of_mass))){
   displayMassPoints <- TRUE
}

#-------------------------------------------------
#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
#-------------------------------------------------
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(ry, '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')))
} else {
   overall_intensityX_normalized <- overall_intensityX
   overall_intensityY_normalized <- overall_intensityY	
}
diff_intensity_normalized <- overall_intensityY_normalized - overall_intensityX_normalized

changes1 <- paste("XX","Overall coverage","Overall coverage","Difference overall coverage","Overall intensity (absolute) ","Overall intensity (absolute)","Difference overall intensity (absolute)","Overall intensity (percentage)","Overall intensity (percentage)","Difference overall intensity (percentage)","Mass centre longitude","Mass centre latitude","Mass centre longitude","Mass centre latitude","Shift distance", sep=",")
changes2 <- paste("species","current","prediction","prediction-current","current","prediction","prediction-current","current_%","prediction_%","prediction-current_%","current","current","prediction","prediction","km", sep=",")
changes3 <- paste(label,overall_coverX,overall_coverY,overall_coverY-overall_coverX,overall_intensityX,overall_intensityY,diff_intensity,overall_intensityX_normalized,overall_intensityY_normalized,diff_intensity_normalized,x_centre[1],x_centre[3],y_centre[1],y_centre[3],dist_km,sep=",")
changes <- paste(changes1,changes2,changes3, sep = "\n")


#-------------------------------------------------
#export of PNG overviews
#-------------------------------------------------
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
}

#print current layer
png(filename=currentLayerPng,height=size, 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))
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 ,breaks=breakpoints,col=cols, alpha=0.7, asp=getAsp(rx), maxpixels=size^2)
} else {
   plot(rx ,breaks=breakpoints,col=cols, alpha=0.7, asp=getAsp(rx), maxpixels=size^2)
}
plot(shapeFile, add=TRUE)
if(isTRUE(displayPoints)){
   points(coordinates(pts),col = "black", bg = "yellow",pch = 21)
}

if(isTRUE(displayMassPoints)){
  points(x_COGx,x_COGy,col = "black", bg = "green",pch = 21, cex = 1.5)
  title(main='currentLayer, Centre of Mass')
  #points(y_COGx,y_COGy,col = "black", bg = "blue",pch = 21)
} else {
  title(main='currentLayer')	
}
dev.off()

#print prediction layer
png(filename=predictionLayerPng,height=size, width=size, bg="white")
if(isTRUE(normalize)){
   plot(ry_normalized ,breaks=breakpoints,col=cols, alpha=0.7, asp=getAsp(rx), maxpixels=size^2)
} else {
   plot(ry ,breaks=breakpoints,col=cols, alpha=0.7, asp=getAsp(rx), maxpixels=size^2)
}
plot(shapeFile, add=TRUE)
if(isTRUE(displayPoints)){
   points(coordinates(pts),col = "black", bg = "yellow",pch = 21)
}

if(isTRUE(displayMassPoints)){
   points(y_COGx,y_COGy,col = "black", bg = "green",pch = 21, cex = 1.5)
   title(main='predictionLayer, Centre of Mass')
} else {
   title(main='predictionLayer')	
}
dev.off()

#print diff layer
png(filename=diffLayerPng,height=size, width=size, bg="white")
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))
#plot(diffxy,breaks=breakpoints_diff,col=cols_diff, alpha=0.7) #main='diffLayer',
plot(diffxy,breaks=breakpoints_diff,col=cols_diff, alpha=0.7, asp=getAsp(rx), maxpixels=size^2)
plot(shapeFile, add=TRUE)
if(isTRUE(displayPoints)){
   points(coordinates(pts),col = "black", bg = "yellow",pch = 21)
}

if(isTRUE(displayMassPoints)){
   points(x_COGx,x_COGy,col = "black", bg = "green",pch = 21,cex = 1.5) #cex = 0.5
   points(y_COGx,y_COGy,col = "black", bg = "green",pch = 21,cex = 1.5)
   arrows(x_COGx, x_COGy, y_COGx, y_COGy, col = "red", lwd=3)
   title(main='diffLayer, Shift Vector')
} else {
   title(main='diffLayer')	
}
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

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