Nothing
read.dmna <- function( file, val=1, debug=FALSE ) {
#
# read header
#
header <- dmna.header(file)
#
# get number and kind of dimensions
#
if ( ! ( "dims" %in% names(header))) {
stop (paste(file,"does not contain number of dimensions"))
}
dims=as.integer(header$dims)
if ( dims > 3 ) {
stop ("this function only supports up to three dimensions")
} else {
if (debug) print (paste("dims:",dims))
}
if ("artp" %in% names(header)) {
artp=header$artp
} else {
artp=""
}
if (debug) print (paste("artp:",artp))
# get index order and orientation
#
# index sequence gives order (slowest counting to fastest counting)
# of numbers in file e.g. "k+,j-,i+"
# index position is position of axis in list seq
# e.g. x-axis boundaries are in first column in lowb/highb
# x-index "i" is found in last position, direction is +
# -> fastest counting, increasing
# -> along data rows, lowes x left highest x right
#
if ( ! ( "sequ" %in% names(header) ) ) {
stop (paste(file,"does not contain information on index order"))
}
sequ=unlist(strsplit(header$sequ,","))
if (length(sequ) != dims) {
stop (paste("number of indices does not match number of dimensions in",file))
}
if (dims==1||dims==2||dims==3) {
# index names
inam=c("i","j","k")
# direction of each index in sequence
# take second character of sequence entry,
# assume "+" if 2nd character is missing
seqdir=substr(paste(sequ[1:dims],"+",sep=""),2,2)
seqind=substr(sequ[1:dims],1,1)
# position of each index in sequence
ipos=rep(NA,dims)
idir=rep("",dims)
for (i in 1:dims){
if (inam[i] %in% seqind){
ipos[i]=grep(inam[i],seqind)
idir[i]=seqdir[ipos[i]]
}
}
} else {
stop (paste(dims," dimensions are not supported by this version"))
}
#
# index boundaries
#
if ( ! ( "hghb" %in% names(header) && "lowb" %in% names(header))) {
stop (paste(file,"does not contain information on grid size"))
}
lowb=as.integer(read.table(header=F,text=header$lowb))
hghb=as.integer(read.table(header=F,text=header$hghb))
ilen=hghb-lowb+1
seqlen=rep(NA,dims)
seqlen[ipos]=ilen
numrec=prod(ilen)
if (debug) {
print(paste('ipos',ipos))
print(paste('idir',idir))
print(paste('ilen',ilen))
print(paste('seqdir',seqdir))
print(paste('seqlen',seqlen))
}
#
# ascii or binary ?
#
if ( "mode" %in% names(header)) {
mode=header["mode"]
} else {
mode="text"
}
if (debug) print(paste('mode:',mode))
#
# how many values per data record
#
if ( "form" %in% names(header)) {
form=header$form
form=gsub("[\t ]+", " ", form)
forms=unlist(strsplit(form," "))
nval=length(forms)
valnams=unlist(lapply(strsplit(forms,"%"), function(l) l[[1]]))
} else {
nval=1
}
if (debug) print (paste("nval:",nval))
#
# select value
#
if (nval>1) {
if (is.character(val)){
if (val %in% valnams) {
ov=which(valnams==val)
} else {
stop(paste('ERROR: value name ',val,' nof in file ',file))
}
} else {
if (val %in% 1:nval) {
ov=val
} else {
stop(paste('ERROR: value number ',val,' not in file ',file))
}
}
if (length(ov)>1) {
stop(paste('ERROR: value specification ',val,' not unique in file ',file))
}
} else {
if (val!=1){
warning(paste('value specification ',val,' not needed in file ',file))
}
ov=1
}
if (debug) print (paste("header lines: ",header$lines))
#
# read the data
#
if ( mode == "text" ){
if ( dims > 1 & nval > 1 ){
stop (paste("ERROR: number of values >1 for dimensions >1 not implemented "))
}
table=read.table(file,skip=header$lines,
header=FALSE,
comment.char = "*",quote=NULL,
as.is=TRUE,blank.lines.skip=TRUE)
# if data section contains in-line comments
# remove supposed comment column and following columns
for (i in 1:dim(table)[2]){
if (any(table[,i]=="'")){
table = table[,1:(i-1)]
break
}
}
# warp table
numbers=as.matrix(table)
values=c(t(numbers))
}else if ( mode == "binary" ){
binfile=gsub("dmna","dmnb",file)
to.read <- file(binfile, "rb")
values=readBin(to.read, "numeric", n=nval*numrec,
size = 4, endian = "little")
close(to.read)
}else{
stop (paste("ERROR: unsopported mode ",mode))
}
# data in the file are in FORTRAN order i.e. last index is counting fastest
# in the array command, the first index is the fastest counting
# => use reverse order of indexes
value=list()
for ( l in 1:nval ){
value[[l]]=array(values[l+seq(0,(nval*numrec-1),nval)],rev(seqlen))
}
#
# rearrange axes to the R-convention (equivalent to sequence "i+,j+,k+")
# i.e. un-reverse order of of indices
#
dat=list()
for ( l in 1:nval ){
dat[[l]]=aperm(value[[l]],rev(ipos))
}
#
# reverse order of values if an index was counting backwards
#
for ( l in 1:nval ){
dat[[l]]=.revdim(dat[[l]],which(idir=='-'))
}
#
# make output
#
if (dims==1) {
out=as.data.frame(do.call(cbind, dat))
names(out)=valnams
#
# if timeseries: find time column and convert to POSIXct
#
for (i in names(out)){
if (artp=="ZA" & i=="te") {
# "." in format string breaks processing !?
out$te=as.POSIXct(gsub("\\."," ",out$te), "%Y-%m-%d %H:%M:%S", tz="")
}
}
#
# try harder to output numbers not factors
#
for (i in names(out)){
# exclude time column
if (i!="te") {
if (is.factor(out[[i]])){
out[[i]] = as.character(out[[i]])
}
if (is.character(out[[i]])){
if (! any(is.na(suppressWarnings(as.numeric(out[[i]]))))) {
out[[i]] = as.numeric(out[[i]])
}
}
out[[i]] = type.convert(out[[i]], as.is = FALSE)
}
}
} else if (artp=="CZ") {
#
# artp CZ is a timeseries though described as as 2D array (why?)
#
if (debug) print ('filetype CZ: adding time column')
te = .timeaxis(header, debug)
out = cbind(as.data.frame(dat), data.frame(te=te))
} else {
out=drop(dat[[ov]])
}
return (out)
}
#
# ------------------------------------------------------------------------
#
dmna.header <- function(file) {
#
# read the file as text lines and find the divider line "*"
#
lines <- readLines(file)
divider=which(lines == "*")
if (length(divider) != 1) {
stop (paste(file,"is not in DMNA format"))
}
#
# convert the file header into named list
#
head_lines=lines[1:(divider-1)]
# remove empty lines
head_lines=head_lines[head_lines!=""]
# convert space behind line tag into tab (if not already present)
head_lines=sub("\ +","\t",head_lines)
# 1st field is name 2nd and on is content
head_names=sub("\t.*","",head_lines)
head_values=sub("[a-z0-9]+\t","",head_lines)
#remove tabs and quotes
head_values=gsub("\t"," ",head_values)
head_values=gsub("\\\"","",head_values)
#make named list
names(head_values)=head_names
head=as.list( head_values )
head$lines=divider
return(head)
}
#
# ------------------------------------------------------------------------
#
dmna.grid <- function(file) {
header <- dmna.header(file)
#
# get axes length
#
if ( ! ( "dims" %in% names(header))) {
stop (paste(file,"does not contain number of dimensions"))
}
dims=header$dims
if ( header$artp == "ZA" | dims<2 ) {
stop (paste("this function does not apply to timeseries"))
}
if ( ! ( "hghb" %in% names(header)
&& "lowb" %in% names(header))) {
stop (paste(file,"does not contain index value ranges"))
}
lowb=read.table(header=F,text=header$lowb)
hghb=read.table(header=F,text=header$hghb)
xlen=hghb[,1]-lowb[,1]+1
ylen=hghb[,2]-lowb[,2]+1
#
# get axes start
#
if ( ! ( "xmin" %in% names(header)
&& "ymin" %in% names(header)
&& "delta" %in% names(header))) {
stop (paste(file,"does not contain all information on axes"))
}
xmin=as.numeric(header$xmin)
ymin=as.numeric(header$ymin)
delta=as.numeric(header$delta)
#
# reference position
#
if ( ! ( "refx" %in% names(header)
&& "refy" %in% names(header))) {
stop (paste(file,"does not contain all information on axes"))
}
refx=as.numeric(header$refx)
refy=as.numeric(header$refy)
#
# lower left corner reference:
#
xll=refx+xmin
yll=refy+ymin
#
# return list
#
out=c(xlen,ylen,xll,yll,delta)
names(out)=c("xlen","ylen","xll","yll","delta")
return(out)
}
#
# ------------------------------------------------------------------------
#
dmna.axes <- function(file, debug=FALSE) {
header <- dmna.header(file)
#
# get axes length
#
if ( ! ( "dims" %in% names(header))) {
stop (paste(file,"does not contain number of dimensions"))
}
dims=header$dims
if ( "axes" %in% names(header)) {
axes=header$axes
}else{
if ( dims == 1) {
axes='ti'
}else if ( dims == 2 ){
axes='xy'
}else if ( dims >= 3 ){
axes='xyz'
}else{
stop (paste(file,"does not contain names of axes"))
}
}
if ( ! ( "hghb" %in% names(header) && "lowb" %in% names(header))) {
stop (paste(file,"does not contain index value ranges"))
}
lowb=read.table(header=F,text=header$lowb)
hghb=read.table(header=F,text=header$hghb)
xlen=hghb[,1]-lowb[,1]+1
if (dims>1) {
ylen=hghb[,2]-lowb[,2]+1
} else {
ylen=1
}
if (dims>2) {
zlen=hghb[,3]-lowb[,3]+1
} else {
zlen=1
}
if (axes == 'ti'){
# time in data:
if (grepl("te%", header$form)){
te = read.dmna(file)$te
} else {
te = .timeaxis(header, xlen = xlen, debug = debug)
}
return(list(te=te))
}else{
#
# get axes start
#
if ( ! ( "xmin" %in% names(header)
&& "ymin" %in% names(header)
&& "delta" %in% names(header))) {
stop (paste(file,"does not contain information on all axes"))
}
xmin=as.numeric(header$xmin)
ymin=as.numeric(header$ymin)
delta=as.numeric(header$delta)
x=xmin+delta*((1:xlen)-1)
if (dims>1) {
y=ymin+delta*((1:ylen)-1)
} else {
y=0:0
}
if (zlen==1) {
return(list(x=x,y=y))
} else {
if ( ! ( "sk" %in% names(header))) {
stop (paste(file,"does not contain level heights"))
}
sk=as.matrix(read.table(header=F,text=header$sk))
if (! ( hghb[,3] <= length(sk) && lowb[,3] >= 1 ) ) {
stop (paste(file,"level indices outside givel level heights"))
}
z=sk[lowb[,3]:hghb[,3]]
return(list(x=x,y=y,z=z))
}
}
}
#
# ------------------------------------------------------------------------
#
# define a helper function following a discussion at
# http://r.789695.n4.nabble.com/distributing-the-values-of-data-frame-to-a-vector-based-on-td837896.html
#
#' @NoRd
.revdim <- function(arr, di=NA) {
if ( is.vector(arr) ){
return(rev(arr))
} else {
ndim = dim(arr)
dims = length(ndim)
idim = 1:dims
rdim = idim %in% di
if ( sum(rdim) > 0 ) {
rind <- function(i){
if (! i %in% idim) {
return(0)
} else if (rdim[i]) {
return(rev(1:ndim[i]))
} else {
return(1:ndim[i])
}
}
do.call("[", c(list(arr), sapply(idim, rind, simplify=FALSE),drop=FALSE))
} else {
return(arr)
}
}
}
#
# ------------------------------------------------------------------------
#
#' @NoRd
.timeaxis <- function(header, xlen = -1, debug = FALSE){
if ( ! ( "dt" %in% names(header)
&& "rdat" %in% names(header))){
stop (paste(file,"does not contain enough information on time axis"))
}
dt = .parsedifftime(header$dt) # "01:00:00"
rd = gsub("([0-9]{4}-[0-9]{2}-[0-9]{2})[ T.]([0-9]{2}:[0-9]{2}:[0-9]{2}).*",
"\\1T\\2",
header$rdat) # "2000-01-01T00:00:00+0100" or "2000-01-01.00:00:00+0100" or "2000-01-01 00:00:00"
rdat = as.POSIXct(rd, format="%Y-%m-%dT%H:%M:%S")
if ( "t1" %in% names(header)){
t1 = .parsedifftime(header$t1) # "00:00:00"
}else{
t1 = difftime(0, units = 'secs')
}
if ( "t1" %in% names(header)){
t2 = .parsedifftime(header$t2) # "366.00:00:00"
}else{
t2 = difftime(0, units = 'secs')
}
if (debug) print (paste(rdat + t1, rdat + t2 - dt, dt))
te = seq(from = rdat + t1, to = rdat + t2 - dt, by = dt)
if (debug) print (paste('... ',strftime(te[1],"%F %T"),' -- ',strftime(te[length(te)],"%F %T")))
return(te)
}
#
# ------------------------------------------------------------------------
#
#' @NoRd
.parsedifftime <- function(s){
# parse difftime format [ddd.]hh:mm:ss
#
# split into days and hh:mm:ss
if (grepl(".", s, fixed=TRUE)) {
# if days are present, use them
p <- unlist(strsplit(s, ".", fixed=TRUE))
d <- p[1]
s <- p[2]
} else {
# if not, use zero days
d <- "0"
}
re1 <- as.difftime(as.numeric(d), units = "days")
re2 <- as.difftime(s, format="%H:%M:%S")
return(re1 + re2)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.