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')
}
Comments (0)
No comments yet
Log in to make a comment