IDF.R 27.4 KB
Newer Older
Rust Henning's avatar
Rust Henning committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
##################################################
## IDF package 
## Authors: Sarah Joedicke, Christoph Ritschel
## Update: 29.01.2016  
###################################################

###############################################
############# Read Data function ##############
###############################################

##' @title Reading precipitation data 
##' @description The function \code{IDF.read} reads a file in table format and creates a \code{list} from it, writing the data frame as first entry 
##' and adds some attributes (station information, aggregation time, data source). The only data values used are: 
##' date, precipitation
##' The \code{data.frame} will have the following format:
##' | year | mon | day | hour | min | RR |
##' |------+-----+-----+------+-----+----+
##' |      |     |     |      |     |    |
##' @usage IDF.read(file, type) 
##' @param file a \code{character string} naming the file from which the data is to be read. 
##' @param type a \code{character string} defining the type of data to be read: either "stadtmessnetz" or "webwerdis", depending on if the data comes from the Stadtmessnetz Berlin
##' or WebWerdis. If type = "webwerdis", the data will be read, then sorted, formatted and missing lines added, 
##' while if type = "stadtmessnetz", the data will just be read and formatted. 
##' Both source types have a different layout in the original file.
##' @return Liste a \code{list} containing a \code{data.frame} of date and time information and precipitation values for each time step
##' @details This function is designed to prepare a data file for aggregation by functions \code{IDF.agg} or \code{IDF.magg}.
##' The time given in the data is the end time, so the precipitation was measured up to that time.  
##' @note The first entry of the list will be a \code{data.frame} named 't1'.
##' @seealso read.table, IDF.agg, IDF.magg
##' @author Sarah Joedicke \email{sarah.joedicke@@fu-berlin.de}

IDF.read <- function(file,type){
  
  if(type != "stadtmessnetz" && type != "webwerdis") {
    
    cat("Warning: wrong type declared for input file")
    stop()
  }
  
  if (type == "stadtmessnetz") {
    
    Tab_MN <- read.csv2(file)  #STADTMESSNETZ
    new_time <- strptime(Tab_MN$Zeitstempel,format="%d.%m.%Y %H:%M")   #STADTMESSNETZ Datumsvektor
  }
  
  # Da die Stadtmessnetzdaten (bisher) konstistent aussehen, wird auf das Erstellen einer neuen Tabelle mit sicher allen
  # Zeiten verzichtet, da die Minutendaten sehr gross sind. Sollte es inkonsistente Tabellen geben, sollte man diese seperat behandeln,
  # sonst wird viel Rechenzeit fuer die kompletten Tabellen verschwendet. 
  
  if (type == "webwerdis") {
    Tab <- read.table(file,header=TRUE,sep=";")   #WEBWERDIS
    Tab_kurz <- Tab[,c("Date","precipitation")]
    
    ## Tabelle in gewuenschtes Format sortieren; nur WEBWERDIS, da Stadtmessnetz sortiert und keine fehlenden Zeilen
    time <- strptime(Tab_kurz$Date,format="%Y-%m-%d T %H:%M:%S")
    Tab_sort <- Tab_kurz[order(as.character(time)),]
    time_sort <- strptime(Tab_sort$Date,format="%Y-%m-%d T %H:%M:%S")
    Tab_sort$Date <- as.character(time_sort)
    
    ## Fehlende Zeilen in Daten hinzufuegen, Niederschlagswert = NA; nur WEBWERDIS
    h_diff <- as.numeric(difftime(format(time_sort[length(time_sort)],"%Y-%m-%d T %H:%M:%S") , 
                                  format(time_sort[1],"%Y-%m-%d T %H:%M:%S"),units="hours")) #h_diff ist die gesamte Zeitdifferenz
    new_time <- seq(time_sort[1], length = h_diff+1, by = "hour")
    new_tab <- data.frame(Date=as.character(new_time), precipitation=NA)  # Tabelle mit NAs und richtigen Zeiten vordefinieren
    
    Tab_na <- (merge(Tab_sort, new_tab, "Date", all.y=TRUE))[,1:2]
  } #Daten sortieren und in gewuenschte Form bringen
  
  new_timect <- as.POSIXct(new_time)
  
  J <- as.numeric(format(new_timect,'%Y'))
  M <- as.numeric(format(new_timect,'%m'))
  d <- as.numeric(format(new_timect,'%d'))
  h <- as.numeric(format(new_timect,'%H'))
  m <- as.numeric(format(new_timect,'%M'))
  
  if (type == "webwerdis") Tab_end <- data.frame(J,M,d,h,m,Tab_na$precipitation.x) #WEBWERDIS
  if (type == "stadtmessnetz") Tab_end <- data.frame(J,M,d,h,m,Tab_MN[,2]) #STADTMESSNETZ
  
  ## Tabellen-Attribute benennen: 
  
  colnames(Tab_end) <- c("year","mon","day","hour","min","RR")
  attr(Tab_end,"accumulation time (min)") <- as.numeric(difftime(new_timect[2],new_timect[1], units="mins"))
  Liste <- list(t1=Tab_end)
  
  if (type == "webwerdis"){
    # WEBWERDIS:
    attr(Liste,"StationName") <- as.character(Tab$Stationname[1])
    attr(Liste,"StationID") <- "NA"
    attr(Liste,"Long (deg N)")  <- Tab$Longitude[1]
    attr(Liste,"Lat (deg E)") <- Tab$Latitude[1]
    attr(Liste,"Heigth (m)")   <- Tab$StationHeight[1]
    attr(Liste,"Source") <- "Web-WERDIS"
  } #Listen-Attribute benennen
  
  if (type == "stadtmessnetz"){
    # STADTMESSNETZ:
    attr(Liste,"StationName") <- colnames(Tab_MN)[2]
    attr(Liste,"StationID") <- "NA"
    attr(Liste,"Long (deg N)")  <- "NA"
    attr(Liste,"Lat (deg E)") <- "NA"
    attr(Liste,"Height (m)")   <- "NA"
    attr(Liste,"Source") <- "Stadtmessnetz"
  } #Listen-Attribute benennen
  
  cat(paste("read.data of", file , "done \n"))
  str(Liste)   # optional; so sieht man beim Einlesen, womit man es zu tun hat und ob alles geklappt hat
  
  return(Liste)
} 



############################################
########## Moving aggragation  #############
############################################

##'@title Moving aggregation of a time series
##'@description The function \code{IDF.magg} aggregates the precipitation values of a \code{data.frame} in a \code{list} by moving an aggregation
##'window over the time steps. The list should be created by funcntion \code{IDF.read} in a previous step. 
##'@usage IDF.magg(x, ks=c(2,4,8,16,32,64,128)) 
##'@param x a \code{list} containing a \code{data.frame} named 't1', preferably generated by function \code{IDF.read}.
##'@param ks a numeric \code{vector} of aggregatation levels. Each value represents an aggregation level and will create 
##'a new \code{data.frame} in the output list.
##'@details This function is designed to aggregate data which has been read by function \code{IDF.read}.
##'The time given in the data is the end time, so the precipitation was added up to that time.
##'Every line of data will be used in this process. In the end, the data frames containing aggregated precipitation 
##'will be exactly as long as the original one. In the end, a new \code{list} containing several \code{data.frames} will be created. 
##'The first one will be the initial \code{data.frame}, now named agg0. The other data frames will be named 'agg' + the aggregation level number from ks.
##'@return Liste a \code{list} containing \code{data.frames} of original and aggregated precipitation time series with time information and 
##'precipitation values.
##'@seealso IDF.read, IDF.agg.
##' @author Sarah Joedicke \email{sarah.joedicke@@fu-berlin.de}

IDF.magg <- function(x,ks=c(2,4,8,16,32,64,128)){
  NS  <- x$t1$Niederschlagssumme   #Niederschlagsvektor
  Liste <- vector("list", length(ks)+1)   # Leere Liste vordefinieren, +1 fuer Ursprungstabelle
  Liste[[1]] <- x$t1
  j <- 2
  for (y in ks){
    agg <- vector(length=nrow(x$t1))
    print(paste("Aggregation ueber Aggregationsstufe", y, "laeuft"))
    for (i in y:nrow(x$t1)){
      agg[i] <- sum(NS[(i-(y-1)):i])
    }
    agg[1:(y-1)] <- NA
    
    Tab_agg           <- x$t1[,1:5]
    Tab_agg[,6]       <- agg
    colnames(Tab_agg) <- c("Jahr","Monat","Tag","Stunde","Minute","Niederschlagssumme (mm)")
    attr(Tab_agg,"Aggregationszeit (min)") <- attr(x$t1, "Aggregationszeit (min)")*y
    Liste[[j]] <- Tab_agg
    j <- j+1
  }
  attributes(Liste) <- attributes(x)
  names(Liste) <- paste("agg", c(0,ks), sep="")
  print("Vorgang abgeschlossen")
  return(Liste)
}

############################################
########## diskret aggregation #############
############################################

##'@title Aggregation of a time series for seperated time steps 
##'@description The function \code{IDF.agg} aggregates the precipitation values of a \code{data.frame} in a \code{list} in groups.
##'The list should be created by function \code{IDF.read} in a previous step. 
##'@usage IDF.agg(x, ks=c(2,4,8,16,32,64,128)) 
##'@param x a \code{list} containing a \code{data.frame} named 't1', preferably generated by function \code{IDF.read}.
##'@param ks a numeric \code{vector} of aggregatation levels. Each value represents an aggregation level and will create 
##'a new \code{data.frame} in the output list.
##'@return Liste a \code{list} containing \code{data.frames} of original and aggregated precipitation time series with time information and 
##'precipitation values.
##'@details This function is designed to aggregate data which has been read by function \code{IDF.read}.
##'The time given in the data is the end time, so the precipitation was added up to that time.
##'In this process the data will be summed up in groups, depending on the chosen aggregation levels. In the end, 
##'the aggregated \code{data.frames} will be much shorter than the original one.
##'In the end, a new \code{list} containing several \code{data.frames} will be created. The first one will be the 
##'initial data frame, now named agg0. The other data frames will be named 'agg' + the aggregation level number from ks.
##'@seealso read.data, agg.f().
##' @author Sarah Joedicke \email{sarah.joedicke@@fu-berlin.de}

IDF.agg <- function(x,ks=c(2,4,8,16,32,64,128)){
  
  NS <- x$t1$RR
  date <- paste(as.character(x$t1$year), as.character(x$t1$mon), as.character(x$t1$day), 
                as.character(x$t1$hour), as.character(x$t1$min))
  Liste <- vector("list", length(ks)+1)   # Leere Liste vordefinieren, +1 fuer Ursprungstabelle
  Liste[[1]] <- x$t1
  j <- 2
  for (y in ks){
    agg_NS <- apply(matrix(data=NS[1:((length(NS)%/%y)*y)], nrow=(length(NS)%/%y), ncol=y, byrow=TRUE), 1, sum)
    # %/% dividiert und rundet ganzzahlig ab
    # Variable if NS%%y > 0 sind das die unbenutzen Resttage, die letzten "" Zeitschritte werden ignoriert
    if (length(NS)%%y > 0) cat("Note: the last", length(NS)%%y,"timesteps are ignored \n")
    agg_date <- matrix(data=date[1:((length(NS)%/%y)*y)], nrow=(length(NS)%/%y), ncol=y, byrow=TRUE)[,y]
    new_time2 <- (strptime(agg_date, format="%Y %m %d %H %M"))
    
    J <- as.numeric(format(new_time2,'%Y'))
    M <- as.numeric(format(new_time2,'%m'))
    d <- as.numeric(format(new_time2,'%d'))
    h <- as.numeric(format(new_time2,'%H'))
    m <- as.numeric(format(new_time2,'%M'))
    
    Tab_agg <- data.frame(J,M,d,h,m,agg_NS)
    colnames(Tab_agg) <- c("year","mon","day","hour","min","RR")
    attr(Tab_agg,"accumulation time (min)") <- as.numeric(difftime(new_time2[2],new_time2[1]),units="mins")
    Liste[[j]] <- Tab_agg
    j <- j+1
  }
  attributes(Liste) <- attributes(x)
  names(Liste) <- paste("agg", c(0,ks), sep="")
  return(Liste)
} # End of function agg
################################################################################################################


##'@title Yearly precipitation maximum  
##'@description The function \code{IDF.max} generates an \code{vector} of yearly precipitation maxima from a \code{data.frame}
##'containing time information and precipitation data at a single aggregation level generated by functions \code{IDF.agg} or \code{IDF.magg}. An option
##'for selection a single month or search the whole year exists.
##'@param df a \code{data.frame} containing date and time information and the corresponding precipitation values at a single aggregation level. 
##'@param month a \code{value} defining the month to use. A single number between 1 and 12 selects the month. The \code{character string} "all"
##'can be used if the maximum for all months is needed.
##'@return RR.max a \code{vector} containing the precipitation maxima for each year of a defined month or all months.
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}

IDF.max <- function(df,month=1) {
  
  years <- unique(df$year)
  y.len <- length(years)
  
  RR.max <- array(NA,dim=c(y.len))
  
  for(y in 1:y.len) {
    
    if(month=="all") {  index <- which(df$year==years[y])
    }else index <- which(df$year==years[y] & df$mon==month)
    if(sum(index)>0)  RR.max[y] <- max(df$RR[index],na.rm=T)
    
  }
  
  RR.max <- RR.max[!is.na(RR.max)]
  
  return(RR.max)
  
}

#######################
## Fitting Functions ##
#######################

##'@title Density function of modified generalized extreme value distribution
##'@description The function \code{dgev.d} is a modified version of the function \code{dgev} for different durations \code{d} developed by Koutsoyiannis et al. (1998).
##'@param q Vector of quantiles
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param d vector of durations
##'@param log \code{logical} option to use logarithmic parameter values, default=FALSE
##'@seealso \code{\link[evd]{dgev}}
##'@return dgev.d gives the density function
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}

dgev.d <- function(q,mu=0,sigma=1,xi=0,theta=0,eta=1,d=1,log=FALSE) {
  sigma.d <- sigma/(d+theta)^eta
  dgev(q,loc=mu*sigma.d,scale=sigma.d,shape=xi,log=log)
}


##'@title Quantile function of modified generalized extreme value distribution
##'@description The function \code{qgev.d} is a modified version of the function \code{qgev} for different durations \code{d} developed by Koutsoyiannis et al. (1998).
##'@param p Vector of probabilities
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param d vector of durations
##'@param lower.tail \code{logical} if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
##'@seealso \code{\link[evd]{qgev}}
##'@return qgev.d gives the quantile function
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
qgev.d <- function(p,mu=0,sigma=1,xi=0,theta=0,eta=1,d=1,lower.tail=TRUE) {
  
  sigma.d <- sigma/(d+theta)^eta
  qgev(p,loc=mu*sigma.d,scale=sigma.d,shape=xi,lower.tail=lower.tail)
  
}

##'@title Random generation for the modified generalized extreme value distribution
##'@description The function \code{rgev.d} is a modified version of the function \code{rgev} for different durations \code{d} developed by Koutsoyiannis et al. (1998).
##'@param n Number of observations
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param d vector of durations
##'@seealso \code{\link[evd]{rgev}}
##'@return rgev.d generates random derivates
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
rgev.d <- function(n,mu=0,sigma=1,xi=0,theta=0,eta=1,d=1) {
  ## gumbel
  sigma.d <- sigma/(d+theta)^eta
  rgev(n, loc=mu*sigma.d,scale=sigma.d,shape=xi)
}

#######################################################################
##' @title Negativ log-likelihood of modified GEV
##' @description The function \code{IDF.nll} calculates the negative log-likelihood for a given set of model parameters
##' \code{mu,sigma,xi,theta,eta}, given observations \code{x} and given durations \code{d}. Options for the usage of
##' logartihmic values \code{use.log} and a debugging function \code{DEBUG} are available.
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param x vector of observations at different durations d
##'@param d vector of durations
##'@param use.log \code{logical} value for usage of logarithmic values, default is \code{FALSE}
##'@param DEBUG \code{logical} value for usage of debugging, if \code{TRUE} the input parameters and the value of negative
##'log-likelihood are printed on console.
##'@return retruns weightes negative log-likelihood by number of observatons uesd
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}

IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,use.log=F,DEBUG=F) {
  ## mu is the mu~ from Koutsoyiannis
  
  if(use.log){
    ## ensure that critical parameters are positive
    sigma <- exp(sigma)
    theta <- exp(theta)
    eta <- exp(eta)
  }
  
  sigma.d <- sigma/((d+theta)^eta)
  if(DEBUG) debug.values <- c(mu,sigma,xi,theta,eta)
  
  
  ## Weibull und Frechet
  if(xi!=0){
    C <- 1 + xi * (x/sigma.d - mu )
    nll <- switch((sum(C<0)>0)+1,
                  sum(log(sigma.d),na.rm=T)+(1+1/xi)*sum(log(C),na.rm=T)+sum((C)^(-1/xi),na.rm=T),
                  Inf)
    #       + penalty*(sum(C[C<0]^2))
    ## Gumbel
  }else if(xi==0){# & sigma<1 & eta<1) 
    Y <- x/sigma.d-mu
    nll <- -(-sum(log(sigma.d))-sum((Y))-sum(exp(-Y)))
  }
  
  if(DEBUG){ 
    cat(debug.values,nll,"\n")
    options(digits.secs=6)
    ##    debug.values <- c(debug.values,nll,as.character(Sys.time()))
    ##    write(debug.values,file="optim.log",append=TRUE,ncolumns=length(debug.values))
    ##    cat(debug.values,nll,sum(A<0),"\n")
  }
  return(nll/length(x))
  
  } # end of function IDF.nll
######################################################################################################

##' @title Fitting function to optimize IDF model parameters
##' @description The function \code{fit.fun} fits IDF model parameters \code{mu,sigma,xi,theta,eta} to a set of given observations \code{obs}, 
##' typically a series of yearly maxima at different durations \code{d}. Options for using logarithmic parameter values and debugging
##' are given. Also the \code{optim} parameters \code{method} and \code{upper,lower} can be defined.
##' @param obs vector of yearly intensity maxima at different durations. Order: Y1D1, Y2D1,...,YnD1,Y1D2,...YnD2,Y1D3,...,YnDk
##' @param dur vector of durations with same length as \code{obs}. Order: n x D1, n x D2, ... n x Dk 
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param use.log \code{logical} value for usage of logarithmic values, default is \code{FALSE}
##'@param DEBUG \code{logical} value for usage of debugging, if \code{TRUE} the input parameters and the value of negative
##'log-likelihood are printed on console for each iteration during optimization.
##' @param method \code{character} defining the method to be used in \code{optim}, preferences are: "Nelder-Mead", "BFGS", "L-BFGS-B"e
##' @param lower \code{vector} specifying the lower boundary of parameters for "L-BFGS-B" method
##' @param upper \code{vector} specifying the upper boundary of parameters for "L-BFGS-B" method
##' @return $min value of negative log-likelihood at optimization minimum
##' @return $par vector of IDF parameters at optimization minimum
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}

fit.fun <- function(obs,dur,mu=1,sigma=1,xi=0.5,theta=1,eta=1,use.log=F,DEBUG=F,method="Nelder-Mead",upper=Inf,lower=-Inf) {

  use.log=use.log
  
  if(use.log) {
    sigma <- log(sigma)
    theta <- log(theta)
    eta <- log(eta)
    
    if(method=="L-BFG-S") {
    upper[2] <- log(upper[2])
    upper[4] <- log(upper[4])
    upper[5] <- log(upper[5])
    
    lower[2] <- log(lower[2])
    lower[4] <- log(lower[4])
    lower[5] <- log(lower[5])
    }
    
  }
  
  ## check initial value of negative log-Likelihood function
  nll <- IDF.nll(mu,sigma,xi,theta,eta,x=obs,d=dur,use.log=use.log,DEBUG=DEBUG)
  
  ## if initial value is acceptable...
  if(!is.infinite(nll)&!is.na(nll)) {
    
    
    if(method=="L-BFGS-B") {
      
      ## problem: optimization algrorithm often has difficulities concerning infinite or NA-difference values betweeen iterations
      ## solution: ignore this error message using functon tryCatch and return NULL if there was an error during optimization
      fit <- tryCatch(mle(IDF.nll,start=list(mu=mu,sigma=sigma,xi=xi,theta=theta,eta=eta),fixed=list(x=obs,d=dur,use.log=use.log,DEBUG=DEBUG),
                          control=list(trace=0,maxit=5000),
                          method=method,upper=upper,lower=lower), error=function(e) NULL)#,
      #upper=upper,lower=lower)
      
    }else{
      
      ## problem: optimization algrorithm often has difficulities concerning infinite or NA-difference values betweeen iterations
      ## solution: ignore this error message using functon tryCatch and return NULL if there was an error during optimization
      fit <- tryCatch(mle(IDF.nll,start=list(mu=mu,sigma=sigma,xi=xi,theta=theta,eta=eta),fixed=list(x=obs,d=dur,use.log=use.log,DEBUG=DEBUG),
                          control=list(trace=0,maxit=5000),
                          method=method), error=function(e) NULL)#,
      #upper=upper,lower=lower)
      
      
      
    }
    
    ## if there was no error
    if(!is.null(fit)) {
      fit.min <- fit@min
      fit.par <- fit@coef
    }else { ## else return NA
      fit.min <- NA
      fit.par <- rep(NA,5)  
    } ## end if error
    
  }else { ## else retunr NA
    
    fit.min <- NA
    fit.par <- rep(NA,5)  
    
  } ## end if initial value..
  
  if(use.log){
    fit.par[2] <- exp(fit.par[2])
    fit.par[4] <- exp(fit.par[4])
    fit.par[5] <- exp(fit.par[5])
  }
  names(fit.par) <- c("mu","sigma","xi","theta","eta")
  
  return(list("min"=fit.min,"par"=fit.par))
  
} ## end of function fit.fun
##################################################################################

##' @title Fitting IDF model parameters to observations at different durations
##' @description The function \code{IDF.fit} fits the IDF model parameters \code{mu,sigma,xi,eta,theta}
##' to a list object of observations \code{data} in a data.frame with temporal inforamtion and values of precipitation
##' at a given temporal resoultion. This precipitation time series gets aggregated at given aggregation levels
##' \code{agg.lev} and yearly maxima of intensity are caluclated for a specific month or the whole year. 
##' The starting values of the IDF model parameters can be determined by the user as well as specific options to use
##' during optimization. Logartihmic transformation, debugging, the optimization method, and an option to plot the
##' IDF curves.
##'@param data a \code{list} containing a \code{data.frame} named 't1', preferably generated by function \code{IDF.read}.
##'@param agg.lev a vector of aggregation levels used to fit the IDF curves.
##'@param month a value specifying the month to be used for estimating the IDF parameters. Type "all" for all months.
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param use.log \code{logical} value for usage of logarithmic values, default is \code{FALSE}
##'@param DEBUG \code{logical} value for usage of debugging, if \code{TRUE} the input parameters and the value of negative
##'log-likelihood are printed on console for each iteration during optimization.
##' @param method \code{character} defining the method to be used in \code{optim}, preferences are: "Nelder-Mead", "BFGS", "L-BFGS-B"e
##' @param lower \code{vector} specifying the lower boundary of parameters for "L-BFGS-B" method
##' @param upper \code{vector} specifying the upper boundary of parameters for "L-BFGS-B" method
##' @param plot \code{logical} option of creating a plot of IDF curves with estimated parameters.
##' @param station.name \code{character} naming of data input, e.g. a observational station name
##' @return $ints vector of sorted intensities for selected aggregation levels
##' @return $durs vector of sorted aggregation levels
##' @return $min minimum value of negative log-likelihood during optimization
##' @return $par vector of estimated IDF model parameters mu,sigma,xi,theta,eta at minimum value of negative log-likelihood.
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}

IDF.fit <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",mu=1,sigma=1,xi=0,theta=0,eta=1,
                    use.log=FALSE,DEBUG=FALSE,method="Nelder-Mead",upper=Inf,lower=-Inf,plot=FALSE,station.name="Berlin") {
  
  data.agg <- IDF.agg(data,ks=agg.lev)
  
  max.agg <- do.call("cbind",lapply(data.agg,FUN=IDF.max,month=month))
  d.all <- c(1,agg.lev)
  
  int <- t(max.agg)/d.all
  y.len <- dim(int)[2]
  int.vec <- as.vector(t(int))
  
  durs <- rep(d.all,each=y.len)
  
  fit <- fit.fun(obs=int.vec,dur=durs,mu=mu,sigma=sigma,xi=xi,theta=theta,eta=eta,use.log=use.log,
                 DEBUG=DEBUG,method=method,upper=upper,lower=lower)
  
  
  if(plot&& !is.na(fit$min)) {
    ds <- sort(rep(d.all,length(int.vec)/length(d.all)))
    IDF.plot(pars=fit$par,name=station.name,ints=int.vec,ds=ds)
  }
  
  if(plot && is.na(fit$min)) {
    cat("Warning: optimization did not converge and no parameters where estimated. Plot not possible. \n")
  }
  
  return(list("ints"=int.vec,"durs"=durs,"min"=fit$min,"par"=fit$par))
  
} ## End of function IDF.fit
######################################################################################################################

##' @title Plotting IDF curves
##' @description The function \code{IDF.plot} plots a set of IDF curves with given IDF model parameters \code{pars} for
##' several probability levels \code{probs} at given durations \code{dur}. The colors of the curves can be defined with
##' parameter \code{cols} (need to have same length as \code{probs}). The \code{station.name} will be printed in the legend.
##' @param pars a vector of IDF model parameters mu,sigma,xi,eta,theta
##' @param probs a vector of probabilities for which the IDF curves are calculated
##' @param dur a vector of durations at which the IDF curves are calculated
##' @param cols a vector of colors for the seperate IDF curves, needs same length as \code{probs}
##' @param name a character string defining the station name for which the IDF curves are calculated, will be printed in legend
##' @param ints \code{vector} of observational intensities (surted by durations)
##' @param ds \code{vector} of durations (same length as intensities)
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}

IDF.plot <- function(pars,probs=c(0.5,0.9,0.99),dur=c(0.5,1,2,3,6,12,24,48,72,96),cols=c(rgb(1,0,0,1),rgb(0,1,0,1),rgb(0,0,1,1)),name="Berlin-Dahlem",ints=NA,ds=NA) {
  
## initialize array for IDF values at different durations and for different probabilities
  idf.array <- array(NA,dim=c(length(dur),length(probs)))
  
  ## loop over probabilities
  for(i in 1:length(probs)) {
    
    ## calculate IDF values for given probability at all durations
    idf.array[,i] <- qgev.d(probs[i],mu=pars[1],sigma=pars[2],xi=pars[3],theta=pars[4],eta=pars[5],d=dur)
    
  } ## end of loop over probs
  
  ## initiialize plot window with limits of IDF values
  plot(NA,axes=F,xlim=c(min(dur),max(dur)),ylim=c(min(idf.array[,1]),max(idf.array[,3])),xlab="duration [h]",ylab="intensity [mm/h]",log="xy")
  axis(1,at=dur,labels=dur)
  axis(2)  
  points(ds,ints,pch=16,col=rgb(0,0,0,0.5))
  
  ## loop over probabilities
  ## plot IDF curve
  for(i in 1:length(probs)) {
    points(dur,idf.array[,i],type="l",col=cols[i],lwd=1.5)
  } 
  
  ## plot legend
  legend(x="topright",legend=c(name,"obs","IDF 0.5 quantile","IDF 0.9 quantile","IDF 0.99 quantile"),
         col=c(1,rgb(0,0,0,0.5),cols[1],cols[2],cols[3]),lty=c(NA,NA,1,1,1),pch=c(NA,16,NA,NA,NA))
  
} ## end of function IDF.plot
###################################################################################